Nonlinear vibrations of thin hyperelastic plates

Abstracts Static deflection as well as free and forced nonlinear vibration of thin square plates made of hyperelastic materials are investigated. Two types of materials, namely rubber and soft biological tissues, are considered. The involved physical (material) nonlinearities are described through Neo-Hookean, Mooney–Rivlin, and Ogden hyperelastic laws; geometrical nonlinearities are modeled by the Novozhilov nonlinear shell theory. Dynamic local models are first built in the vicinity of a static configuration of interest that has been previously calculated. This gives rise to the approximation of the plate's behavior in the form of a system of ordinary differential equations with quadratic and cubic nonlinear terms in displacement. Numerical results are compared and validated in the static case via a commercial FE software package: they are found to be accurate for deflections reaching 100 times the thickness of the plate. The frequency shift between low- and large-amplitude vibrations weakens with an increased initial deflection.


Introduction
Thin-walled structures made of hyperelastic materials, such as rubbers and biomaterials, are common in mechanical and biomechanical engineering applications. Such structures are frequently subjected to significant static and dynamic loadings potentially yielding large deflections and deformations. Associated previous works reported in the literature employ, in a vast majority, a simplifying assumption that considers as a priori known the shape of the deformed structure; for instance, a cylindrical shell remains cylindrical after deformation. A general review on these works is available in [1]. The present investigation discards this limitation.
One possible way to account for complicated deformed shapes is to use the finite element (FE) approach. For instance, Einstein [2] developed a finite element strategy to analyze the dynamic behavior of hyperelastic and viscoelastic membranes with amplitudes of the order of their thickness. With the proposed method, the dynamic response of a hemisphere made of a Mooney-Rivlin material under impact pressure load is studied. A relatively similar problem concerning the dynamic inflation of a Mooney-Rivlin spherical membrane is explored in [3] also through the FE method.
Another technique lies in the approximation of the deformed configuration as a truncated series of continuous functions satisfying the essential geometric boundary conditions. This approach is implemented for circular membranes in [4,5]. The sensitivity of the respective vibratory behavior to the pre-stretch configuration of the membranes is carried out in these scientific contributions. The incompressible hyperelastic material is either described by a Neo-Hookean model [5] or by Mooney-Rivlin, Yeoh, Ogden, and Arruda-Boyce models in [4]. In this latter work, the frequencyamplitude relationships for various hyperelastic models are found to be analogous. The authors come to the conclusion that the higher the stretch, the closer the nonlinear forced response to its linear counterpart.
In [1], the static deflection and vibration around a deformed configuration of a rectangular Neo-Hookean plate are investigated accounting for both geometrical (the nonlinearity of straindisplacements relations) and physical nonlinearities (the material nonlinearity of stress-strain relations). Is proposed a method which systematically builds approximate local models (LMM) in the form of polynomial expansions of the non-polynomial strain energy densities: this greatly simplifies the description of the plate's dynamic behavior. The present paper extends this investigation to various hyperelastic laws by using a more sophisticated nonlinear plate theory. The two targeted types of material (rubber and biomaterial) are described through (i) Neo-Hookean, (ii) Mooney-Rivlin, and (iii) Ogden hyperelastic laws as well as (iv) a physically linear material law. The comparison of static force-deflection curves and frequency responses obtained with these models is carried out and the contribution of in-plane nonlinearities in the sought displacement is estimated. The numerical solutions are verified through direct comparison with results from a commercial FE software package.

Kinematics and energy equations
The geometrical nonlinearity in the kinematics of the investigated flexible plate is described with the help of the Novozhilov nonlinear plate theory [6], which stands as a limit case of the Novozhilov nonlinear shell theory. It is governed by the following strain-displacement relationships: where " 1 , " 2 , and " 12 are the components of the Green-Lagrange strain tensor for thin plates. The Novozhilov nonlinear plate theory, commonly recognized as the finest classical plate theory, reduces to the well-known von Kármán nonlinear plate theory [6] when the in-plane terms in square brackets in (1) are neglected.
Although it neglects shear deformations and rotary inertia, it is accurate for thin plates. Also, it is known that in-plane nonlinearities play a major role in large deflections, and they are accordingly retained in the present work to improve accuracy. Lagrange equations are utilized to derive the dynamics of the plate, that is: where L D T … is Lagrange's functional, T is the kinetic energy of the plate, … is the potential elastic deformation energy, and Q n are the generalized forces. The potential and kinetic energies are expressed as follows [6]: where W is the strain energy density (SED), V is the volume of the plate, S is the surface of the middle plane of the plate, is the mass-density of the plate material, h is the thickness of the plate, and u; v; w are the displacements along the axes of the rectangular coordinate system x; y; z, respectively. The dot stands for differentiation with respect to time. The displacements are expanded into corresponding truncated series involving the appropriate generalized coordinates q n : In Eq. (5), quantities W i ; U i ; V i are the admissible functions that satisfy the homogeneous boundary conditions (i.e. the geometric constraints) of the problem. The linear modes of vibration, which form a complete set properly capturing the dynamics of a structure, are eligible admissible functions and are selected in the present work. The total number of degrees of freedom is

Strain energy density
Usually, nonlinear elasticity of rubbers and soft biomaterials is described by hyperelastic laws and in most cases, such materials are assumed to be incompressible [7,8]. Three hyperelastic laws together with their associated SED W are considered: Neo-Hookean: Mooney-Rivlin: Ogden: Also, the strain energy density for a linear material reads: The following notation is used: I 1 is the first invariant of the right Cauchy-Green deformation tensor C; E is Young's modulus of the plate's material; I 2 stands for the second invariant of the right Cauchy-Green deformation tensor; 1 ; 2 ; 3 are the principal stretches of the plate; i ;˛i denote the material parameters. Since SED (9) is polynomial in strains, spatial and temporal components of the solution can be uncoupled and the Lagrange equations reduce to simple ordinary differential equations with quadratic and cubic nonlinearities [6]: k nij`qi q j q`D Q n ; n D 1; : : : ; N where n is the natural frequency of mode n and n is the corresponding damping ratio; k ni ; k nij ; k nij`a re known coefficients that result from integration in space.

Cauchy-Green tensor invariants and principal stretches
In order to derive the expressions of the invariants of the right Cauchy-Green strain tensor C in terms of displacements, the Green-Lagrange strain tensor: is used, where the expressions of " 1 , " 2 , and " 12 (but not " 3 ) are given in (1). The right Cauchy-Green deformation tensor C is then defined as [8]: and its three invariants are: The third invariant J is used to reflect the incompressibility condition through J D 1 [7] and the principal stretches are the square roots of the eigenvalues of C [8]:

Strain energy density local expansion
Expressions (6)-(8) together with " 3 from (18) are not polynomials in strains, which essentially complicates the investigation of the plate's behavior. The analysis is thus simplified by introducing a transformation of SEDs (6)-(8) in order to derive approximate governing equations in the form of ordinary differential equations with nonlinearities of order not higher than three. The corresponding local model is reliable only in the vicinity of a configuration of interest around which the SED is expanded into a series in the generalized coordinates truncated at order 4.
To reach highly deformed configurations, successive local models have to be constructed. Further details on the local models method (LMM) can be found in [1].

Exact low-dimensional models with both material and geometric nonlinearities
A distinct additional approach is implemented to measure the accuracy of the LMM. It involves the numerical solution of Equa-tion (2) in its static version, that is: • V @W @q n dV D Q n ; n D 1; : : : ; N with SEDs (6)- (8) and expression (18) previously substituted into it. Since the latter is a smooth function in the generalized coordinates, integration and differentiation operators commute, that is: System (19) of nonlinear algebraic equations can be solved numerically through the Newton-Raphson iterative technique; this approach is named the "exact solution" hereinafter.
6. Numerical example: static and dynamic bending of a rubber plate

Hyperelastic models parameter identification
Experimental data for 8% sulfur rubber obtained by Treloar [9] for uniaxial and equibiaxial tensions, as well as for pure shear, are exploited. This common rubber is chosen because experimental data for uniaxial and multiaxial loads is available in the literature: a good approximation for the strain energy density is then possible for this material. Corresponding details on the approximation of the experimental data are provided in Appendix A. From the given experimental data, Young's modulus is E D 1 247 060:2 Pa, which is in almost perfect agreement with Ogden's one, E D 1 242 992:9 Pa [7]. Other relevant parameters are listed in Table 1.

Model Parameters
Neo-Hookean E D 1 247 060:2 Pa  2, and 3 display the strain-stress relationships for the three hyperelastic laws together with the experimental points. The nominal stress S 1 versus engineering strain " is shown-see Appendix A for details. It is worthy to note that the Neo-Hookean and Mooney-Rivlin laws are very close and both laws satisfactorily approximate the rubber behavior at strains not higher than 30% (that is 0:3 in Figures 1 to 3). However, the approximation of the rubber behavior offered by the Ogden law is much more accurate on the full range of experimental strains.

Problem description
A simply supported square rubber plate illustrated in Fig. 4 is considered. It is defined on the following domain: with a D 0:1 m, b D 0:1 m, h D 0:0005 m. The material characteristics are listed in Table 1 and the rubber mass density is D 1100 kg m 3 .  The plate is simply supported with immovable edges yielding the following boundary conditions [6]: where @S denotes the boundary of the plate's middle surface. The bending moment per unit length M [6] reads: and quantities n and are the outer normal and tangent directions to @S , respectively, as shown in w.x; y; t/ D X n;m2N u.x; y; t/ D X n;m2N The problem is made non-dimensional by introducing a nondimensional time D 1 t , where 1 is the circular frequency of the first natural mode of the deformed (pressure loaded) plate. Also, a change of notation has been introduced; the two-subscript generalized coordinates (time functions) w 2nC1;2mC1 , u 2n;2mC1 , and v 2nC1;2m are divided by the plate thickness h and replaced by the single-subscript generalized coordinates q i with corresponding increasing numbering.

Static analysis
First, the problem of static plate bending under uniform pressure is explored. As illustrated in Figure 5, a quick convergence analysis shows that a 12 degree-of-freedom (DOF) model stands as a convincing compromise between prediction capabilities and computational cost for the force-deflection curve of the Neo-Hookean plate. The generalized coordinates plugged into the expansion of the displacements are given in Table 2. Responses are calculated with the LMM, starting from the configuration with w 1;1 D 20h captured by the less sophisticated model involving geometrical nonlinearity only. Accordingly, the curves start at a central deformation of 20h in Figure 5. The maximal discrepancy between the 12 DOF and 27 DOF models in the deflection range OE0; 100h is 2.3%. Similar results are found for other material models.
Size Participating eigenmodes 3 DOF w 1;1 , u 2;1 , v 1;2 12 DOF w i;j i; j D 1; 3 u i;j , v j;i i D 2; 4, j D 1; 3 27 DOF w i;j i; j D 1; 3; 5 u i;j , v j;i i D 2; 4; 6, j D 1; 3; 5 Table 2. Generalized coordinates utilized in the models Force-deflection relationships for different materials with 12 DOF are depicted in Figure 6. The displacement sensitivity to physical nonlinearities is highlighted by plotting the solution for linear elastic material and geometric nonlinearities only. Neo-Hookean and Mooney-Rivlin results almost coincide, while the Ogden curve slightly deviates for large deflections. Also, it is worthy to note that the LMM results agree very well with the exact solution.
Material nonlinearity for deflections smaller than 20h can be neglected as shown in Figure 6. For deflections up to 20h the model with only geometrical nonlinearity is a very good approximation. As a consequence, the configuration computed with the linear elastic material is used as an initial guess for LMM numerical iterations for hyperelastic materials [1]. An additional validation of the results is carried out. The static deflection of the hyperelastic plate is also explored with the commercial finite element solution ANSYS [11] for the Neo-Hookean and Ogden materials. The central deflection with respect to the applied external force is displayed in Figure 7: FEA and LMM predictions are in good agreement.
The contribution of the in-plane nonlinearities in expressions (1) is now estimated. As an illustrative example, the pressure-deflection curves for the Neo-Hookean law with and without in-plane nonlinearities are shown in Figure 8. The maximal difference is 2%. Other materials feature discrepancies of the same order. Despite of their relatively small influence, in-plane nonlinearities are incorporated in all numerical investigations except in Figure 8.

Nonlinear vibration analysis
At small strains, the effect of physical nonlinearities is negligible ( Figure 6), and the vibrations of the initially flat plate can be explored by only retaining geometrical nonlinearities in the model [6]. Here, we focus on the investigation of a more challenging configuration where both nonlinearities are participating, i.e. the vibrations around a pre-loaded state. The deflection with principal generalized coordinate w 1;1 D 80h is chosen as an initial deformed configuration. Comparison with the exact static solution shows that the local model around this deformed configuration is accurate for deflections up to 10h, which therefore stands as an upper limit in the dynamic analysis. From previous analysis ( Figure 5), it follows that the 12-DOF model is sufficiently accurate and can be selected in order to reduce computation costs. The harmonic balance method [12] is implemented to find the sought periodic solutions to system (10) through a Fourier expansion of the generalized coordinates in time: A nj cos.j /; n D 1; : : : ; N where is the non-dimensional frequency, normalized with respect to 1 . Expression (27) considers cosine terms only since the external force and damping in (10) are ignored in the free vibration analysis. Coefficients A nj are determined from the system of nonlinear algebraic equations resulting from balancing the coefficients associated to the same harmonics in equations (10) where expression (27) has been previously substituted. A conver-  the deformed plate are given in Table 3 for an initial deflection w 1;1 D 80h. Figure 9 displays the backbone curves of the plate's free vibrations, with frequencies close to the first eigenfrequency of the pre-loaded plate and normalized with respect to the natural frequency 1 of the corresponding deflected plate. Comparison with the exact static solution shows that the Neo-Hookean and Mooney-Rivlin local models are accurate for deflections up to 10h, while the Ogden model can be considered for deflections not higher than 8h only. Again, the Neo-Hookean and Mooney-Rivlin models yield identical results. This is true for both the eigenfrequency of the deformed plate and the backbone curve. The Ogden backbone curve exhibits a slightly weaker nonlinearity. For all materials, the nonlinear nature of the deformed plate, as opposed to its flat counterpart [6], is very weak: the frequency of large-amplitude vibrations is almost identical to the natural frequency. In order to show the sensitivity of the backbone curves to the initial deflection, they are constructed around various initial configurations in Figure 10. The Neo-Hookean material is selected since other materials exhibit similar features. Two effects are observed: the first one consists of an increased range of vibra- Ogden [ ] tory amplitudes where the backbone curve displays a softening behavior. This effect is well-known and might be attributed to the contribution of the quadratic terms in expression (10) with respect to the initial deflection amplitude [6]. The second effect is the attenuation of the nonlinear nature of the model with an increased initial deflection. A similar behavior has been previously reported for bended plates [1] and stretched membranes [7] and is related to the large in-plane stretching associated to the initial deflection of the loaded plate. The forced vibrations are also explored with the AUTO package [13]. External forcing derives from a time-dependent periodic pressure having an harmonic component P d D 4:53 Pa and a constant mean-value corresponding to the deformed configuration w 1;1 D 80h, around which the local model is built. A modal damping ratio n D D 0:001 is adopted. The frequency response for the principal bending coordinate w 1;1 is shown in Figure 11. The vibratory response around a highly deformed configuration is almost identical to the linear response, in contrast to the response around a moderately deflected plate for which a significant softening behavior turning to hardening for vibration amplitude around 5h is found as shown in Figure 12.

Hyperelastic models parameters
The experimental data for uniaxial tension available in [14] are exploited. The present study is intended to model only one key feature of biological materials, i.e. a sharp increase in stiffness after a given strain threshold is reached. To this end, the experimental data corresponding to the tunica adventitia of a human aorta in [14] are selected. The procedure described in Appendix A is implemented with the difference that only the uniaxial test data are available. Corresponding material parameters are listed in Table 4. The associated stress-strain relationships are shown in Figures 13 to 15. The difference between the Neo-Hookean and Mooney-Rivlin laws is clearly distinguishable for this type of material. However, both laws are in agreement with the experimental points for strains less than 8% and both are unable to reproduce the increase in stiffness, as opposed to the Ogden law.

Static analysis
The problem detailed in subsection 6.2 is solved for the hyperelastic materials listed in Table 4. Attention is paid to the static bending of the plate (Figure 4) with the help of the 12 DOF model ( Table 2). The exact pressure-deflection response is depicted in Figure 16. The Neo-Hookean and Mooney-Rivlin results are almost identical, while the Ogden model starts to depart from these curves for amplitudes about 80h, as displayed in Figure 16, where it becomes much stiffer. Strains at the point where the Ogden and the Neo-Hookean curves cross each other in Figure 16 are also analyzed through the strain intensity, given by [15]:  (Figures 13 and 14).

Vibration analysis
We investigate the free nonlinear vibrations around the deformed response with principal generalized coordinate w 1;1 D 70h. The comparison with the exact static solution shows that the Neo-Hookean and Mooney-Rivlin local models are accurate for deflections up to 10h whereas the Ogden model is limited to deflections not larger 3h only. The corresponding backbone curves are obtained through the harmonic balance method (27) and are depicted in Figure 17 for vibrations around the first eigenfrequency of the pre-loaded plate. The first eigenfrequency of the deformed Model Natural Frequency Neo-Hookean 183 rad/s Mooney-Rivlin 177.88 rad/s Ogden 181.96 rad/s plate is listed in Table 5 for the three hyperelastic laws. The mass density is D 1380 kg m 3 , which is the density of Polyethylene terephthalate, frequently used for artificial arteries [16]. We notice that the Neo-Hookean and Mooney-Rivlin results are again in good agreement for the full investigated range of deflections. The Ogden plate's eigenfrequencies are also similar for central deflections up to 80h but then become higher. Similarly to the rubber plate example, the behavior of a deformed plate made of biomaterial is weakly nonlinear for all materials and the nonlinearity is of softening type. However, certain features differ. Although static results for the Neo-Hookean and Mooney-Rivlin models are almost identical, the Mooney-Rivlin backbone curve shows a softer behavior. Still, the Ogden curve is the softest one and is presented for vibration amplitude up to 3h in Figure 17.

Conclusions
Vibration of plates made of rubber and biological materials are explored with a dedicated local models method. Local models allow for the analysis of static bending as well as physically and geometrically large-amplitude nonlinear vibrations of an initially distorted plate. Static results are validated through a systematic comparison with available exact solutions as well as against commercial finite element software results. It is found that the local models method provides accurate predictions for a wide range of deflections. However, the sharp increase in stiffness peculiar to biological materials limits the range of achievable vibration amplitudes. In most cases, it is found that the Mooney-Rivlin and Neo-Hookean materials exhibit similar static and vibratory behaviors. Corresponding constitutive laws properly capture the behavior of the actual material at moderate strains (30% for rubber and 8% for biomaterial). The best approximation is provided by the Ogden model. The latter correctly reproduces the behavior at high strains, including the well-known sharp increase in stiffness. However, the Ogden model has significant drawbacks. Due to the complicated form of the strain energy density and its formulation in terms of principal strains (unlike the Neo-Hookean and Mooney-Rivlin models which allow for a formulation in terms of strain invariants), it is computationally much more expensive. It is also shown that the pre-loaded plate exhibits very weak dynamic nonlinearity, i.e. the frequencies of the oscillations around the deformed configuration are almost identical to the associated eigenfrequencies.
The sensitivity of the free vibration backbone curves to the initial deflection is also discussed. It is shown that the higher the initial static deflection, the higher the range of amplitudes at which the backbone curve displays a softening behavior. The nonlinear nature of the system of interest (the frequency shift between low-and large-amplitude vibrations) weakens with an increased initial deflection.