Viscoelasticity in Large Deformation Analysis of Hyperelastic Structures

In this paper, an annular/circular plate made of hyperelastic material and considering the viscoelastic property was investigated based on a novel nonlinear elasticity theory. A new approach for hyperelastic materials in conjunction with the Kelvin–Voigt scheme is employed to obtain the structure’s large deformation under uniform transverse loading. The constitutive equations were extracted using the energy method. The derived partial differential time-dependent equations have been solved via the semi-analytical polynomial method (SAPM). The obtained results have been validated by ABAQUS software and the available paper. In consequence, a good agreement between the results was observed. Finally, several affecting parameters on the analysis have been attended to and studied, such as the nonlinear elasticity analysis, the boundary conditions, loading, and the material’s viscosity. It can be possible to obtain the needed time for achieving the final deformation of the structure based on the applied analysis in this research.


Introduction
Various materials can be found in nature that demonstrate different mechanical behaviors under applied loads, such as metals displaying linear elastic behavior until yielding. At the same time, irreversible deformation is then observed up to failure. On the other hand, brittle materials exhibit a slight elastic deformation linearly and then fail without irreversible deformation. However, it is notable that all materials do not exhibit linear elasticity. Hyperelastic materials will exhibit elastic deformations excessively before failure without irreversible deformation. These materials exhibit an extremely nonlinear stressstrain behavior that ascends monotonically up to fracture. It is well-known that linear elastic materials can be defined with two material constants: modulus of elasticity and Poisson's ratio. On the contrary, hyperelastic materials can be defined by a strain-energy density function. The stored strain energy remains constant when a sample is subjected to constant strain. Therefore, hyperelastic materials are modeled in the context of strain energy potentials such as Arruda-Boyce [1], Gent [2], Neo-Hookean [3], and Yeoh [4].
Some recent studies on the mechanical responses of hyperelastic structures are summarized below. Erchiqui et al. [5] implemented a dynamic finite element method to model the visco-hyperelastic behaviors of a thin, isotropic, and incompressible thermoplastic membrane based on the Ogden and Mooney-Rivlin models. Kocatürk and Akbaş [6] examined the geometrically nonlinear static response of a hyperelastic simply supported beam under a non-follower load by the finite element and Newton-Raphson iteration methods. Li et al. [7] perused the dynamic behavior of the visco-hyperelastic dielectric elastomer structures based on the Gent hyperelastic model. Alibakhshi et al. [8] analyzed the nonlinear vibration of dielectric elastomer microbeam resonators based on a hyperelastic Cosserat  [9] studied the thermomechanical analysis of a hyperelastic thick cylindrical pressure vessel by giving analytical and numerical solutions. Asgari and Hashemi [10] developed an efficient visco-hyperelastic constitutive model for a hollow cylinder elastomer under dynamic and impact loadings. Pascon [11,12] implemented a finite element formulation based on a two-dimensional beam element to investigate various viscoelastic functionally graded materials and beams made of functionally graded hyperelastic material. The constitutive equations are extracted based on the neo-Hookean model. Gharooni and Ghannad [13] studied the nonlinear analysis of functionally graded tapered hyperelastic cylindrical pressure vessels subjected to non-uniform pressure load.
Hosseini and Rahimi [14] conducted the nonlinear bending analysis of a neo-Hookean hyperelastic plate based on the Mindlin plate theory. Xu et al. [15] introduced the plate element formulation by quadratic interpolation. They carried out static and dynamic analyses of incompressible hyperelastic silicone plates. Dadgar-Rad and Firouzi [16] presented a nonlinear finite element formulation for the viscoelastic deformation of hyperelastic structures under several loadings and boundary conditions. Ansari et al. [17] developed a numerical approach to survey the deformations of hyperelastic Mindlin rectangular plates in compressible and nearly incompressible regimes based on the Neo-Hookean model. Tashiro et al. [18] analyzed blood clots using a nonlinear viscoelastic and hyperelastic model. They employed the visco-hyperelastic finite element method to estimate the mechanical behavior of blood clots. Shariyat and Abadi [19] studied the nonlinear dynamic and impact responses of incompressible neo-Hookean hyperelastic plates with stiff elastic reinforcing particles. Runge-Kutta time integration and penalty methods are used in the solution. Karimi et al. [20] conducted the nonlinear dynamic analysis of an embedded neo-Hookean hyperelastic membrane under a uniformly distributed hydrostatic pressure. Alibakhshi and Heidari [21] carried out the nonlinear vibration of a dielectric elastomer balloon via the Gent hyperelastic model. A time integration-based solver is utilized to solve the resulting equations. Alibakhshi et al. [22] elaborated on a suitable detection mechanism for scanning the surface profile of a micro-sample by AFM, considering that the probe is made of a hyperelastic material. Falope et al. [23] proposed a finite element-based theoretical model to study the bending behavior of a hyperelastic solid.
Hosseini et al. [24] perused the nonlinear static behavior of functionally graded hyperelastic plates based on FSDT. The potential energy function is formulated according to the neo-Hookean model and the Cauchy-Green tensor. Coda et al. [25] elaborated a finite element formulation to analyze the laminated and functionally graded hyperelastic one-dimensional structure with transverse shear stress distribution. Dastjerdi et al. [26] introduced a comprehensive theoretical method to examine the mechanical responses of hyperelastic structures. The proposed method can analyze the geometrically and physically nonlinear hyperelastic materials. Zhao et al. [27] investigated the nonlinear dynamic responses of the visco-hyperelastic spherical shells under uniform radial loads. Additionally, Zhao et al. [28] perused dynamic loads and structural damping influences for incompressible hyperelastic spherical shells. They considered the Yeoh strain energy function in their study. Bacciocchi and Tarantino [29] examined the finite bending of hyperelastic beams based on the compressible Mooney-Rivlin model. More information about the nonlinear dynamics of hyperelastic structures can be found in the recent review article [30]. Additionally, many studies have been carried out on analyzing viscoelastic structures [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45]. Li et al. [46] presented a perturbation approach for the lateral vibration analysis of viscoelastic microstructures. The clamped microbeam has been subjected to external harmonic excitation. Recently, some researchers developed a structure-preserving approach to solve the macroscopic/microscopic coupling dynamic problems with large deformation [47]. Additionally, a large deformation analysis of a nano-sized structure has been attended by Yan et al. [48].
As stated before, many recent studies have been performed on the mechanical analysis of structures made of hyperelastic materials based on various strain energy functions, especially the neo-Hookean model. However, there is no study on the static response of an annular circular visco-hyperelastic plate. The new approach for hyperelastic materials simulates the annular circular visco-hyperelastic plate in conjunction with the Kelvin-Voigt scheme. The general formulations are derived according to the first-order shear deformation theory (FSDT). The constitutive equations and boundary conditions are extracted by the energy method. Then, the semi-analytical polynomial method is utilized for solving the partial differential time-dependent equations. Finally, the effects of various parameters on the bending analysis of the annular circular visco-hyperelastic plate are investigated and discussed in detail.

Geometry of the Structure
First, the geometry of the problem is discussed. The geometry of the analyzed structure will play a significant role in obtaining the governing equations. The more comprehensive the geometry of the problem, the more different structures can be modeled only by considering a single formulation. Analysis of the structures made of hyperelastic materials will make the simulation more complicated. Therefore, an attempt has been made to focus more on the effect of viscosity on hyperelastic materials. Hence, an annular-circular structure made of visco-hyperelastic material has been assumed. The geometry of the analyzed problem is shown in Figure 1, along with the considered coordinate system, which is the cylindrical coordinate system. As stated before, many recent studies have been performed on the mechanical analysis of structures made of hyperelastic materials based on various strain energy functions, especially the neo-Hookean model. However, there is no study on the static response of an annular circular visco-hyperelastic plate. The new approach for hyperelastic materials simulates the annular circular visco-hyperelastic plate in conjunction with the Kelvin-Voigt scheme. The general formulations are derived according to the first-order shear deformation theory (FSDT). The constitutive equations and boundary conditions are extracted by the energy method. Then, the semi-analytical polynomial method is utilized for solving the partial differential time-dependent equations. Finally, the effects of various parameters on the bending analysis of the annular circular visco-hyperelastic plate are investigated and discussed in detail.

Geometry of the Structure
First, the geometry of the problem is discussed. The geometry of the analyzed structure will play a significant role in obtaining the governing equations. The more comprehensive the geometry of the problem, the more different structures can be modeled only by considering a single formulation. Analysis of the structures made of hyperelastic materials will make the simulation more complicated. Therefore, an attempt has been made to focus more on the effect of viscosity on hyperelastic materials. Hence, an annular-circular structure made of visco-hyperelastic material has been assumed. The geometry of the analyzed problem is shown in Figure 1, along with the considered coordinate system, which is the cylindrical coordinate system. As can be seen, the inner radius of the structure is , and its outer radius is . Additionally, the thickness of the structure is constant and equal to ℎ. The structure is placed on an elastic foundation with two components, the Winkler ( ) and Pasternak ( ). The external load on the structure is considered a uniform distributed load in the direction equal to .

Viscoelastic Property
In this research, an attempt has been made to perform two nonlinear elastic properties simultaneously with the changes in the deformation of the material with respect to the strain rate. In other words, large strains against the strain rate are considered. The structure's material is considered in such a way that the deformation created in it due to the applied load is different during the duration of the load application. On the other hand, the final deformation of the structure does not happen as soon as the load is applied, and it takes place over time. The Kelvin-Voigt model is used to consider the viscoelastic properties of the material [43]. Additionally, other models such as Zener and mixing linear As can be seen, the inner radius of the structure is R i , and its outer radius is R o . Additionally, the thickness of the structure is constant and equal to h. The structure is placed on an elastic foundation with two components, the Winkler (k w ) and Pasternak (k p ). The external load on the structure is considered a uniform distributed load in the z direction equal to q z .

Viscoelastic Property
In this research, an attempt has been made to perform two nonlinear elastic properties simultaneously with the changes in the deformation of the material with respect to the strain rate. In other words, large strains against the strain rate are considered. The structure's material is considered in such a way that the deformation created in it due to the applied load is different during the duration of the load application. On the other hand, the final deformation of the structure does not happen as soon as the load is applied, and it takes place over time. The Kelvin-Voigt model is used to consider the viscoelastic properties of the material [43]. Additionally, other models such as Zener and mixing linear viscoelasticity and nonlinear elasticity (as has been done in this paper) have been used by other researchers [49,50]. In this type of time-dependent modeling for stress-strain, a spring and a damper are placed in parallel, which can be seen in Figure 2. As shown in Figure 2, the viscous property of the material (time-dependent strain) is simulated by a damper. Additionally, the material's elasticity is modeled by a spring (E represents Young's elasticity modulus). In linear elastic materials, the value of E is constant and changes linearly. However, this value is not fixed or variable in this research. Therefore, the model (as explained in the next part) uses linear viscoelasticity according to the Kelvin-Voigt scheme. However, the principle of elasticity has been considered nonlinearly. viscoelasticity and nonlinear elasticity (as has been done in this paper) have been used by other researchers [49,50]. In this type of time-dependent modeling for stress-strain, a spring and a damper are placed in parallel, which can be seen in Figure 2. As shown in Figure 2, the viscous property of the material (time-dependent strain) is simulated by a damper. Additionally, the material's elasticity is modeled by a spring ( represents Young's elasticity modulus). In linear elastic materials, the value of is constant and changes linearly. However, this value is not fixed or variable in this research. Therefore, the model (as explained in the next part) uses linear viscoelasticity according to the Kelvin-Voigt scheme. However, the principle of elasticity has been considered nonlinearly. Stress and strain will be written according to the following equations based on the Kelvin-Voigt simulation.
It can be seen in Figure 2 that is the elasticity modulus, and is the viscosity of the material. The amount of is linear here according to Equation (1). However, it is not constant in this study due to the nonlinear elasticity analysis. This issue will be discussed in detail in the next part.

Nonlinear Elastic Material and Governing Equations
In a linear elastic material, it can be observed that the changes in stress and strain are linear. However, in a structure with nonlinear elastic property, the stress increases nonlinearly with the increase of strain. In this research, the material with nonlinear elastic properties is assumed. Therefore, calculations based on nonlinear elasticity theory are considered. Various mathematical simulations and theories have been presented regarding the analysis of nonlinear elastic or hyperelastic structures [2,[9][10][11][12]. A new method has been presented previously for analyzing the mechanical behavior of nonlinear elastic structures [26]. According to this theory, the stress-strain diagram of a nonlinear elastic material obtained from the experiment is approximated by a polynomial function ( ( ) = ∑ =1 ). According to the polynomial degree, a higher accuracy can be obtained by choosing more . Consequently, the expressed viscous stress can be formulated as . The calculations related to the strain field and the governing equations can be obtained according to the following equations by considering the viscoelastic property. Stress and strain will be written according to the following equations based on the Kelvin-Voigt simulation.
It can be seen in Figure 2 that E is the elasticity modulus, and g is the viscosity of the material. The amount of σ s is linear here according to Equation (1). However, it is not constant in this study due to the nonlinear elasticity analysis. This issue will be discussed in detail in the next part.

Nonlinear Elastic Material and Governing Equations
In a linear elastic material, it can be observed that the changes in stress and strain are linear. However, in a structure with nonlinear elastic property, the stress increases nonlinearly with the increase of strain. In this research, the material with nonlinear elastic properties is assumed. Therefore, calculations based on nonlinear elasticity theory are considered. Various mathematical simulations and theories have been presented regarding the analysis of nonlinear elastic or hyperelastic structures [2,[9][10][11][12]. A new method has been presented previously for analyzing the mechanical behavior of nonlinear elastic structures [26]. According to this theory, the stress-strain diagram of a nonlinear elastic material obtained from the experiment is approximated by a polynomial function (σ(ε) = According to the polynomial degree, a higher accuracy can be obtained by choosing more n. Consequently, the expressed viscous stress can be formulated as σ s (ε) = The calculations related to the strain field and the governing equations can be obtained according to the following equations by considering the viscoelastic property.
The strain field is assumed based on FSDT as In the above equations, u 0 and w 0 are transport displacements, and ϕ is the rotation function around the θ axes. Now, using the energy method [43][44][45] , the mathematical description of boundary conditions and governing equations can be introduced as the following equations.
δϕ : ∂M rr ∂r Considering the uniform transverse load on the structure, the problem is symmetrical, and there will be no changes in the θ direction. The only independent variables of the problem are r and t, where t represents the time duration. The mathematical definition of the different types of boundary conditions is introduced according to the following relations: Clamped (C), Simply supported (S), and Free (F) at the edges of r = R i , R o .

Solution Method
This research uses the semi-analytical method based on polynomials (SAPM) to solve the governing equations [43,44]. According to this method, the displacement functions (u 0 , w 0 , and ϕ) are formulated as comprehensive polynomials based on the number of nodes (N and M) in each direction of the independent variable of the problem (r and z).
Now, a system of algebraic equations can be obtained by inserting the functions (Equations (10)- (12)) in the mathematical definition of the boundary conditions at the edges and constitutive equations. Finally, by solving the obtained algebraic equations, the unknown functions in the displacement field and, accordingly, other unknowns, including stresses and strains, will be obtained for the visco-hyperelastic annular-circular sheet under uniform transverse load. First, the accuracy of the used solution method is checked. Figure 3 shows the background changes in terms of time for an annular/circular sheet with the following specifications for different values of the number of nodes in the direction of r and t. Now, a system of algebraic equations can be obtained by inserting the functions (Equations (10)- (12)) in the mathematical definition of the boundary conditions at the edges and constitutive equations. Finally, by solving the obtained algebraic equations, the unknown functions in the displacement field and, accordingly, other unknowns, including stresses and strains, will be obtained for the visco-hyperelastic annular-circular sheet under uniform transverse load.

The Solving Method (SAPM)
First, the accuracy of the used solution method is checked. Figure 3 shows the background changes in terms of time for an annular/circular sheet with the following specifications for different values of the number of nodes in the direction of and .  As can be seen, the accuracy of the results increases versus the increase of nodes in the calculations. A sudden jump is observed between the results of = = 3 and = = 5, so the used solution method has a very high convergence. It can be seen that the results of seven and nine nodes are very close, and therefore, the results obtained from seven or nine nodes can be used with reasonable confidence. Increasing the number of nodes in the calculations will significantly increase the time to obtain the results. Therefore, the optimized number of nodes in the network for solving the problem is an important point that should be considered. In the first form, by choosing only nine nodes in each direction ( and ), you can get the results with the desired accuracy for the subsequent analysis.  Table 1 are provided to validate the results. Figure 4 shows the maximum deflection between the results of (a) the present study and (b) the obtained results As can be seen, the accuracy of the results increases versus the increase of nodes in the calculations. A sudden jump is observed between the results of N = M = 3 and N = M = 5, so the used solution method has a very high convergence. It can be seen that the results of seven and nine nodes are very close, and therefore, the results obtained from seven or nine nodes can be used with reasonable confidence. Increasing the number of nodes in the calculations will significantly increase the time to obtain the results. Therefore, the optimized number of nodes in the network for solving the problem is an important point that should be considered. In the first form, by choosing only nine nodes in each direction (r and t), you can get the results with the desired accuracy for the subsequent analysis.  Table 1 are provided to validate the results. Figure 4 shows the maximum deflection between the results of (a) the present study and (b) the obtained results of ABAQUS software for an annular/circular plate. It is demonstrated that the results agree well. Finding relevant papers in visco-hyperelasticity to fit with the present work is hard. For example, some papers can be found. However, their modeled geometry might be different from the present study. Another assessment was done, and the results can be observed in Figure 5 for the visco-hyperelastic structure. As the tested model [51] is a thick circular sheet, we considered a quasi-three-dimensional analysis [52]. Therefore, the captured results have been monitored into a single plot in conjunction with the experimental results [51]. The compared results in Figure 5 show that an acceptable difference is available between the obtained results. The obtained results in this paper are near to the ABAQUS results, according to the originally depicted figure in [51]. Eventually, the applied theory and solving method are reliable. Additionally, according to Table 1, the same conclusion can be made for different boundary conditions. of ABAQUS software for an annular/circular plate. It is demonstrated that the results agree well. Finding relevant papers in visco-hyperelasticity to fit with the present work is hard. For example, some papers can be found. However, their modeled geometry might be different from the present study. Another assessment was done, and the results can be observed in Figure 5 for the visco-hyperelastic structure. As the tested model [51] is a thick circular sheet, we considered a quasi-three-dimensional analysis [52]. Therefore, the captured results have been monitored into a single plot in conjunction with the experimental results [51]. The compared results in Figure 5 show that an acceptable difference is available between the obtained results. The obtained results in this paper are near to the ABAQUS results, according to the originally depicted figure in [51]. Eventually, the applied theory and solving method are reliable. Additionally, according to Table 1, the same conclusion can be made for different boundary conditions.    of ABAQUS software for an annular/circular plate. It is demonstrated that the results agree well. Finding relevant papers in visco-hyperelasticity to fit with the present work is hard. For example, some papers can be found. However, their modeled geometry might be different from the present study. Another assessment was done, and the results can be observed in Figure 5 for the visco-hyperelastic structure. As the tested model [51] is a thick circular sheet, we considered a quasi-three-dimensional analysis [52]. Therefore, the captured results have been monitored into a single plot in conjunction with the experimental results [51]. The compared results in Figure 5 show that an acceptable difference is available between the obtained results. The obtained results in this paper are near to the ABAQUS results, according to the originally depicted figure in [51]. Eventually, the applied theory and solving method are reliable. Additionally, according to Table 1, the same conclusion can be made for different boundary conditions.

Comparisons
(a) (b)  Comparison between the results of this study and [51] for the visco-hyperelastic analysis. Figure 5. Comparison between the results of this study and [51] for the visco-hyperelastic analysis.

Numerical Results
To investigate different boundary conditions on nonlinear-elastic analysis, Figure 6 depicts a circular sheet with the same characteristics as Figure 3. As can be seen, the maximum deflection increases with the increase of transverse load on the structure. The increase in the maximum deflection is completely nonlinear and will be accompanied by a decreasing slope. In other words, as the transverse load on the sheet increases, it is observed that the slope of the changes decreases. The importance of this matter is that if the linear analysis were considered, the obtained results would increase linearly. These linear results would be acceptable only for low load values. However, because the nonlinear analysis of the sheet made of material with nonlinear elastic properties has been studied in this research, the simulation results can be used with appropriate accuracy. Another point (according to Figure 6, which is drawn for two boundary conditions, CC and SS, at the edges of R i and R o ) is that the decreasing trend of the results at the beginning of loading is more for the CC boundary conditions than SS. However, with increasing load on the plate, the changes are almost the same, with a slight difference for the two boundary conditions, CC and SS. In other words, by increasing the load on the structure, the effects related to the boundary conditions are reduced, and the behavior of the sheet against the applied load will be less affected by the applied boundary conditions.

Numerical Results
To investigate different boundary conditions on nonlinear-elastic analysis, Figure 6 depicts a circular sheet with the same characteristics as Figure 3. As can be seen, the maximum deflection increases with the increase of transverse load on the structure. The increase in the maximum deflection is completely nonlinear and will be accompanied by a decreasing slope. In other words, as the transverse load on the sheet increases, it is observed that the slope of the changes decreases. The importance of this matter is that if the linear analysis were considered, the obtained results would increase linearly. These linear results would be acceptable only for low load values. However, because the nonlinear analysis of the sheet made of material with nonlinear elastic properties has been studied in this research, the simulation results can be used with appropriate accuracy. Another point (according to Figure 6, which is drawn for two boundary conditions, and , at the edges of and ) is that the decreasing trend of the results at the beginning of loading is more for the boundary conditions than SS. However, with increasing load on the plate, the changes are almost the same, with a slight difference for the two boundary conditions, and . In other words, by increasing the load on the structure, the effects related to the boundary conditions are reduced, and the behavior of the sheet against the applied load will be less affected by the applied boundary conditions. The viscoelastic property in nonlinear elasticity analysis is a significant issue considered in this research. Therefore, Figures 7 and 8 show the effect of viscosity (changes in deformation over time) on the results of the maximum deflection of the circular sheet (specifications of Figure 3). Viscosity = 0 means that, when the load is applied on the sheet, the maximum deflection will occur, and the structure will reach the absolute limit of deformation, which can be seen in Figure 7. Of course, it can be observed that the mentioned result will be for = 1 . This calculation error is due to the limitation in choosing the number of nodes in time when solving the governing equations. It can be observed The viscoelastic property in nonlinear elasticity analysis is a significant issue considered in this research. Therefore, Figures 7 and 8 show the effect of viscosity (changes in deformation over time) on the results of the maximum deflection of the circular sheet (specifications of Figure 3). Viscosity g = 0 means that, when the load is applied on the sheet, the maximum deflection will occur, and the structure will reach the absolute limit of deformation, which can be seen in Figure 7. Of course, it can be observed that the mentioned result will be for t = 1s. This calculation error is due to the limitation in choosing the number of nodes in time when solving the governing equations. It can be observed that with the rise of the viscosity value (g), the needed time for the sheet to reach its final deformation will increase. In other words, with the increase in viscosity, the slope variations will decrease over time. Of course, the effect of viscosity on the results is not linear. In other words, the distance between the results from g = 0 to g = 5 is more remarkable than from g = 5 to g = 10. As the viscosity of the structure material increases, the intensity of its impact on the results decreases. Figure 8 is one of the results of the curves in Figure 7 (g = 5), which here show the changes in the deflection in two directions: r and t, simultaneously. According to Figure 8, the maximum deflection on the sheet will occur in its middle (r = R i , R o ), according to the CC boundary conditions. tions will decrease over time. Of course, the effect of viscosity on the results is not linear. In other words, the distance between the results from = 0 to = 5 is more remarkable than from = 5 to = 10. As the viscosity of the structure material increases, the intensity of its impact on the results decreases. Figure 8 is one of the results of the curves in Figure 7 ( = 5), which here show the changes in the deflection in two directions: r and , simultaneously. According to Figure 8, the maximum deflection on the sheet will occur in its middle ( = , ), according to the boundary conditions.  In Figures 7 and 8, the effect of viscosity on the nonlinear elastic structure is investigated. As seen, with the increase in viscosity, the time for the structure to reach its final shape increases. Now, this topic, how much time is required for the final deflection (considering a specific viscosity), is further investigated. Figure 9a is drawn for viscosity = 5, and Figure 9b is drawn for = 10. Here, the criterion for drawing two shapes is to reach the absolute limit of the deformation and be equal to a specific number. In other words, the final deflection for both Figure 9a,b is equal. It can be seen that the final deflection in Figure 9a is for = 25 . However, this result in Figure 9b is equal to = 50 . This analysis is crucial, because it is observed that, with the doubling of the viscosity of the tions will decrease over time. Of course, the effect of viscosity on the results is not linear. In other words, the distance between the results from = 0 to = 5 is more remarkable than from = 5 to = 10. As the viscosity of the structure material increases, the intensity of its impact on the results decreases. Figure 8 is one of the results of the curves in Figure 7 ( = 5), which here show the changes in the deflection in two directions: r and , simultaneously. According to Figure 8, the maximum deflection on the sheet will occur in its middle ( = , ), according to the boundary conditions.  In Figures 7 and 8, the effect of viscosity on the nonlinear elastic structure is investigated. As seen, with the increase in viscosity, the time for the structure to reach its final shape increases. Now, this topic, how much time is required for the final deflection (considering a specific viscosity), is further investigated. Figure 9a is drawn for viscosity = 5, and Figure 9b is drawn for = 10. Here, the criterion for drawing two shapes is to reach the absolute limit of the deformation and be equal to a specific number. In other words, the final deflection for both Figure 9a,b is equal. It can be seen that the final deflection in Figure 9a is for = 25 . However, this result in Figure 9b is equal to = 50 . This analysis is crucial, because it is observed that, with the doubling of the viscosity of the In Figures 7 and 8, the effect of viscosity on the nonlinear elastic structure is investigated. As seen, with the increase in viscosity, the time for the structure to reach its final shape increases. Now, this topic, how much time is required for the final deflection (considering a specific viscosity), is further investigated. Figure 9a is drawn for viscosity g = 5, and Figure 9b is drawn for g = 10. Here, the criterion for drawing two shapes is to reach the absolute limit of the deformation and be equal to a specific number. In other words, the final deflection for both Figure 9a,b is equal. It can be seen that the final deflection in Figure 9a is for t = 25s. However, this result in Figure 9b is equal to t = 50s. This analysis is crucial, because it is observed that, with the doubling of the viscosity of the structure, the time it takes to reach the final shape change also doubles. Although the variation rate is nonlinear, the deformation of the structure at the final time has a linear and direct relationship with the value of the viscosity of the structure. Of course, this conclusion has been reached for this particular problem. Such a result may not be obtained by considering other assumptions, including environmental factors such as temperature, humidity, or electric fields. Investigating the impact of the mentioned factors can also be considered by other researchers who work in this field. structure, the time it takes to reach the final shape change also doubles. Although the variation rate is nonlinear, the deformation of the structure at the final time has a linear and direct relationship with the value of the viscosity of the structure. Of course, this conclusion has been reached for this particular problem. Such a result may not be obtained by considering other assumptions, including environmental factors such as temperature, humidity, or electric fields. Investigating the impact of the mentioned factors can also be considered by other researchers who work in this field.

Conclusions
The present study investigated the viscoelastic analysis of nonlinear elastic (hyperelastic) materials. The governing partial differential equations and mathematical definitions of boundary conditions were derived based on a new method and were solved by the semi-analytical method based on polynomials (SAPM). The viscoelastic property of the structural material is assumed by Kelvin-Voigt modeling. In general, the significant results obtained from the research can be categorized as follows: