Nonlinear continuum mechanics of thick hyperelastic sandwich beams using various shear deformable beam theories

In this study, the time-dependent mechanics of multilayered thick hyperelastic beams are investigated for the first time using five different types of shear deformation models for modelling the beam (i.e. the Euler–Bernoulli, Timoshenko, third-order, trigonometric and exponential shear deformable models), together with the von Kármán geometrical nonlinearity and Mooney–Rivlin hyperelastic strain energy density. The laminated hyperelastic beam is assumed to be resting on a nonlinear foundation and undergoing a time-dependent external force. The coupled highly nonlinear hyperelastic equations of motion are obtained by considering the longitudinal, transverse and rotation motions and are solved using a dynamic equilibrium technique. Both the linear and nonlinear time-dependent mechanics of the structure are analysed for clamped–clamped and pinned–pinned boundaries, and the impact of considering the shear effect using different shear deformation theories is discussed in detail. The influence of layering, each layer’s thickness, hyperelastic material positioning and many other parameters on the nonlinear frequency response is analysed, and it is shown that the resonance position, maximum amplitude, coupled motion and natural frequencies vary significantly for various hyperelastic and layer properties. The results of this study should be useful when studying layered soft structures, such as multilayer plastic packaging and laminated tubes, as well as modelling layered soft tissues.

In many engineering applications, structures are fabricated with three layers (e.g. in laminated packaging); the outer layers are called the skins (upper skin and lower skin) and the middle layer is the core. If the purpose of the design is to increase the strength against bending loads, the core is usually made of a softer material with a higher thickness and the outer skins are thinner and stiffer. If the structure is used in an uncertain environment, the outer skins are made of softer materials to propagate the force throughout the structure and avoid cracks in a brittle core. In some other engineering structures, researchers gradually change the material properties to fabricate functionally graded structures which have been studied in Refs. [3][4][5][6].
Over the past few years, researchers have analysed different layered elastic (versus hyperelastic in this paper) structures to understand their mechanical behaviour under diverse conditions. A brief review of previous studies in this field is given in this section; however, interested readers seeking a more detailed literature review are referred to Refs. [7][8][9][10].
For bending analysis of layered elastic beams [11,12], researchers have used different beam theories and analysed layered beam structures. By neglecting the shear effects, Chai and Yap [13] examined the bending behaviour of layered composite structures using classic beam theory (CBT). By comparing the results with those of a finite element modelling (FEM), it was shown that the coupling terms between the layers play a significant role in predicting the mechanical behaviour accurately. Chen et al [14] added the shear influence to the model using the first-order shear deformable beam theory; it was shown that the bending deformation of the Timoshenko beam model is higher than the CBT and the difference is higher for any lower thicknessto-length ratio. Özütok and Madenci [15] used a higher-order shear deformation theory together with FEM for timed-independent analysis of layered composite beams and verified the results by comparing them with those available in the literature.
For linear vibration analysis of layered elastic beams [16][17][18], Emam and Nayfeh [19] investigated the free vibration of layered composites using the Euler-Bernoulli beam theory. The axial motion was considered in the formulation and was substituted in the transverse equation of motion. Natural frequencies were obtained for different axial loads and boundary conditions; it was reported that for the studied six-layer beam model, the layer positioning has a significant effect in varying the fundamental frequency for different axial loads. A first-order shear deformable beam theory was used by Banerjee and Sobey [20] to examine the free oscillation of laminated structures under axial loading. A closed-form analytical solution was presented showing accuracy in predicting the natural frequencies and mode shapes of layered beams. Damanpack and Khalili [21] used higher-order beam theory to model the longitudinal-transverse coupled motion of layered structures. Hamilton's principle was used to reach the equations of motion showing that the natural frequencies obtained from this model agree with the experimental results presented in Ref. [22].
For nonlinear vibration analysis of elastic beams, Zhang et al. [23] studied the nonlinear forced oscillation of laminated beams in a humid thermal condition. Only the transverse motion was considered in the formulation, and a nonlinear energy sink was added to the model. The equations of motion were solved using a fourth-order Rung-Kutta (RK4) method and a two-term harmonic balance (HB) scheme. It was shown that the nonlinear studied model has a hardening behaviour and adding a nonlinear energy sink can reduce the maximum amplitude of the transverse vibration significantly. Shen et al. [24] examined the nonlinear free oscillation behaviour of layered graphene-strengthened beams; Halpin-Tsai model was employed to predict the physical properties of the combination of the matrix and fibres. A higher-order shear deformable theory, together with the von-Kármán geometrical nonlinearity, was used to model the beam. By employing a two-step perturbation technique, it was shown that the layering and positioning of fibres has a significant effect in varying the nonlinear frequency ratio and the natural frequencies. Three-layered and bi-layered shear deformable elastic beams have been studied by Ghayesh et al. [25]; it was shown that the resonance region change based on the layer sorting and both hardening and softening behaviour can be observed while considering geometrical imperfections. Other types of layered structures have also been examined by researchers, and interested authors are referred to Refs. [26,27] for more information about layered-plate and layered-shell structures focused investigations.
The previously mentioned studies focus on elastic structures that undergo small strains while facing different types of loadings; even the nonlinear ones are based on geometric nonlinearities and small strains. Soft structures cannot be modelled as a linear elastic material and require more accurate hyperelastic assumptions and formulations; in other words, a material nonlinearity should be taken into account, as done in this paper.
Multilayered hyperelastic structures have been used widely in many products. For instance, multilayer plastic packaging is used for packaging different products, such as food packaging which has been a trending topic over the last few years in waste management, recycling and optimised packaging [28][29][30][31]. Furthermore, human body organs act as hyperelastic components and many parts are made of different layers, such as the artery, which is made of three hyperelastic layers (intima, media and adventitia) [32].
Previously, axially moving hyperelastic structures and porous-hyperelastic structures have been investigated by Khaniki et al. [33,34]; however, until now, there have been no studies on the time-dependent mechanics of layered hyperelastic structures. In this study, a comprehensive analysis of the mechanics of multilayered thick hyperelastic structures is presented for the first time as having different boundary conditions via different shear theories. Large strain modelling, together with large-amplitude displacements, is modelled using the Mooney-Rivlin strain energy model (for material nonlinearity) and the von Kármán geometrical nonlinearity. Different types of shear deformable beam theories (Euler-Bernoulli, Timoshenko, third-order, trigonometric and exponential shear deformable models) are used to model the layered structures and the influence of layering, material positioning and the thickness of the core and skins is discussed in detail.

Layered hyperelastic shear deformable beam formulation
For a thick-layered shear deformable beam (Fig. 1), considering the plane motion, the axial and transverse deformations are written as [35] u where u and w are the axial and transverse displacements in the x and z directions ( Fig. 1), respectively, subscript x indicates derivations with respect to x, φ is the rotation and χ(z) is the higher-order shear function, which the definition depends on the shear deformable theory used. In this study, five different definitions for χ(z) are used, which are defined as [36][37][38][39] Case 1 (Euler−Bernoulli beam theory) : Case 2 (Timoshenko beam theory) : Case 3 (Third−order beam theory) : Case 4 (Exponential beam theory) : Case 5 (Trigonometric beam theory) : where Case 1 indicates that the shear deformation effect is neglected, following the Euler-Bernoulli beam assumption (classic beam theory) when the rotary inertia is also neglected; the advantage of using this model is that the model is simplified by neglecting the rotational degree of freedom by only considering axial and transverse motions. Case 2 defines a linearly varying shear deformation through the thickness of the beam, which represents the Timoshenko beam theory. In this model, the transverse normal plates do not necessarily remain perpendicular (to the mid-surface) after deformation leading to an extra degree of freedom for the rotary inertia. Case 3 shows the shear deformation using a higher-order model following the third-order shear deformable beam theory [36,37]. In this model, the transverse shear stress in the Cases 4 and 5 defines the shear deformation of the beam using an exponential [38] and trigonometric [39] functions, which follow the beam theories of exponential shear deformable and trigonometric, respectively. For all the five cases, the strains are written as where ε x x is the axial strain, γ xz is the shear strain and ε zz is the transverse strain that is not equal to zero and should satisfy the incompressibility condition of hyperelastic structures. The deformation gradient (F) is given in Refs. [33,40] which can be used for defining the left Cauchy-Green strain tensor (C) for planar motion as where the higher-order strain terms are due to the soft behaviour of hyperelastic materials, which undergo large strains and cannot be ignored. The strain term ε zz can be obtained as a function of ε zz and ε zz by considering det(C) = 0 (the incompressibility condition). For linear elastic models, by neglecting the higher-order strain terms, the Cauchy-Green strain tensor will be [27] which is equal to the elastic material strain tensor model. By using the Mooney-Rivlin hyperelastic strain energy density model [41,42], the variation of the total potential energy of the system is written as where for which, by assuming the hyperelastic beam has three layers ( Fig. 1) with layers attached perfectly (with no delamination), the stiffness parameters (A ii , B ii , D ii and E ii ) are defined as which, for the case considering the Timoshenko beam model (χ(z) = z), uses a shear correction factor of κ = 5/6 [43]. Accordingly, by using Eqs. (19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31)(32), the potential energy of a three-layered hyperelastic beam can be rewritten, which is not presented here for the sake of brevity. It should be mentioned that for having more layers in the laminated hyperelastic beam model, Eqs. (19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31)(32) should be rewritten on the right side, based on the number of layers. The kinetic energy of the structure is written as which, by defining the moment of inertia in terms of area (I i ) for a three-layered hyperelastic thick beam, is written as and by substituting Eqs. (34)(35)(36)(37)(38)(39) into Eq. (33), one can reach the kinetic energy of layered thick hyperelastic beams, as Moreover, the beam is assumed to be lying on a nonlinear elastic foundation ( Fig. 1) and externally actuated by a periodic load, which gives the external work on the beam as where K L , K N L and K S are the linear, nonlinear and shear coefficients of the foundation, respectively, and F is the external force acting periodically with a frequency of ω. Using Hamilton's principle, the coupled longitudinal, transverse and rotation equations of motion are obtained as It can therefore be seen that the equations of motion are highly nonlinear and coupled, which emphasises the importance of having the axial motion while analysing the mechanical behaviour of the structure. The linear coupling sources between the axial and transverse motions are the stiffness terms B 11 and B 22 and inertia terms I 1 and I 2 , for which, for a single-layer hyperelastic beam, due to the homogeneity in the thickness direction, these terms will be equal to zero and the linear coupling terms vanish. In the next section, the equations of motion are nondimensionalised, discretised and solved.

Solution procedure
For the first step of solving the equations of motion, the new nondimensional parameters are defined as and by using these terms, the equations of motion are nondimensionalised as For the sake of brevity and ease notation, the superscript * is neglected in the formulation. Using the Galerkin scheme, each motion component can be written in a series expansion of orthogonal spatial functions as where U n , n , W n are the set of trial functions that should satisfy the boundary conditions of longitudinal, rotation and transverse motions at the middle of the thickness. (Studies on layered structures based on various thickness assumptions can be found in Refs. [44][45][46][47][48][49][50].) For pinned-pinned and fixed-fixed boundary conditions, these trial functions are written as where where M i j are the mass tensors of the discretised equations and the components are defined as and KL i j are the linear stiffness tensors, defined as KN1, KN2 and KN3 are the nonlinear stiffness terms in the coupled equation where and the coefficients of the nonlinear stiffnesses (KN i j ) are defined in "Appendix A" for the sake of brevity. By employing a dynamic equilibrium technique, the time-dependent terms of the degrees of freedom are written in a series expansion of exponential functions as and the general dynamic complex equilibrium equations are written as By solving the discretised coupled 3N equations of motion presented in Eq. (79) using a continuation technique [51][52][53], the nonlinear forced time-dependent mechanical response of the layered hyperelastic structure is obtained.

Results and discussion for linear and nonlinear analyses
The coupled motion of shear deformable layered hyperelastic structures is formulated, discretised and solved in the previous sections. In this section, the linear and nonlinear mechanics of the system are analysed and discussed in two main subsections.

Linear analysis for different shear models
In this subsection, the natural frequencies of the layered hyperelastic structures are analysed and the influence of layering, shear deformation theory and the foundation is discussed in detail.

Model verification
In the first step, in order to verify the current model for analysing layered beam models, the natural frequencies of a three-layered elastic beam structure are obtained for clamped-clamped and pinned-pinned boundaries. Layering is achieved by having steel as the top layer, aluminium in the middle and silicon carbide at the bottom (the elastic material properties of the layers can be found in Ref. [54]), with geometry properties at a length of L = 2.4 m, total thickness of h = 0.03 m (0.01 m for each layer) and width of b = 0.12 m. The first four transverse natural frequencies and the first axial and rotational natural frequencies are obtained using the third-order shear deformable beam theory and compared with those obtained by ANSYS Workbench [55] in Table 1, and the mode shapes are shown in Fig. 2; it can be seen that the current modelling shows great accuracy in computing the natural frequency terms.

Layer sorting, slenderness ratio and shear model effect
Properly layering the structure can have a significant effect in changing the natural frequencies of a hyperelastic beam. To make this possible, three different soft materials are assumed with known Mooney-Rivlin hyperelastic coefficients, as shown in Table 2. For the first step, by assuming a single-layer homogeneous hyperelastic beam with the length of 0.2 mm and nondimensional foundation of K L = 1 and K S = 1, the slenderness ratio (r = L/ h) is varied and the fundamental natural frequency is obtained for different shear deformation theories in Table 3. It can be seen that for the single-layer model, the fundamental natural frequency of all the five shear models is in good agreement, especially for the vulcanised rubber beam. By increasing the number of layers to two and three, the fundamental natural frequencies are given in Tables 4 and 5, respectively. It can be seen that by having more than a single layer, the homogeneity in the thickness direction is lost, the coupling terms are considerable and, therefore, the difference between different shear models is slightly more than for the homogeneous single-layer model. Furthermore, it can be seen that for all the layered models and boundary conditions, by increasing the slenderness ratio (L/ h), the shear effect loses its importance and the results for different shear models become unique.

Foundation effect
To show the influence of a foundation support on the natural frequencies of layered hyperelastic beams, the structure is assumed to be three-layered  Tables 6, 7 and  8

for layering as [thermoplastic-silicone-vulcanised rubber], [thermoplastic-vulcanised rubber-silicone] and
[silicone-thermoplastic-vulcanised rubber], respectively, using the third-order shear deformable beam model. It can be seen that the layered hyperelastic beam models are more sensitive to variations of the shear term (K s ) compared with the linear term (K L ), and adding the foundation increases the fundamental frequency parameter significantly.

Nonlinear analysis
The nonlinear dynamics of layered hyperelastic thick beam structures is investigated in this section for various layer, shear and material position models.     Table 7 Fundamental natural frequencies of three-layered hyperelastic beams [thermoplastic-vulcanised rubber-silicone] with different slenderness ratios and shear deformation theories  Table 8 Fundamental natural frequencies of three-layered hyperelastic beams [silicone-thermoplastic-vulcanised rubber] with different slenderness ratios and shear deformation theories

Comparison of different shear models
Five different shear deformable beam theories were used in previous sections for modelling and formulation of the layered hyperelastic beam structure. In this section, the amplitude response of the sandwich structure for different shear models is examined. Three layers are assumed, and the thickness of the layers is assumed as h = [h/4, h/2, h/4] with L/ h = 10, K L = K N L = 1, K s = 0.1, with the material sorting as [siliconevulcanised rubber-thermoplastic]. In the first step, the influence of considering the shear effect as negligible (CBT), linear (Timoshenko) or higher order (third order) is examined, and the axial, transverse and rotation amplitudes are obtained, as shown in Figs. 3, 4 and 5 for clamped-clamped and Figs. 6, 7 and 8 for pinnedpinned boundary conditions. It can be seen that for both boundary conditions, the Timoshenko beam model gives the highest maximum amplitude for the dominant axial, transverse and rotation deformation. Furthermore, increasing the shear effect from the linear model to higher-order moves the resonance peak to lower frequencies.
The influence of using different types of higher-order shear deformable theories on the amplitude response is investigated for the similar three-layered hyperelastic beams with the given properties. Figures 9, 10, 11, 12, 13 and 14 show the axial, transverse and rotation amplitudes of clamped-clamped and pinned-pinned thick beam models, respectively. The results for the different higher-order shear models are fairly in good agreement, especially for the pinned-pinned boundary condition. Comparing the amplitude response of the dominant terms of axial, transverse and rotation motions, it can be seen that the highest amplitude is for the third-order beam model and the lowest belongs to the exponential beam model.

Influence of layering
Layering a hyperelastic structure for different purposes can change the dynamic behaviour of the structure significantly. In Sect. 4.1.1., this influence on the linear vibration response is analysed, and in this section, the nonlinear vibration response of the structure based on the number of layers is examined. To this end, for a homogeneous single-layer structure, the dominant amplitude responses are shown in Fig. 15, for thermoplastic, vulcanised rubber and silicone hyperelastic thick beams using the third-order shear deformable beam theory model. It can be seen that thermoplastic and silicone show a stiffness hardening effect, while the vulcanised rubber beam model shows a drop in the amplitude response of the axial motion while increasing the frequency and continues to go up afterwards. By increasing the layers to two by having h 1 = h 2 = h/2, the dominant longitudinal, transverse and rotation amplitude responses of the structure are shown for different hyperelastic materials in Fig. 16 Fig. 17. It can be seen that the material distribution [mat2-mat1-mat3] has the lowest maximum axial amplitude and the highest transverse and rotation amplitudes; this shows that for the other two material models, the coupling between the axial and transverse motion is noticeably higher and should be considered.

Influence of layer thickness
Changing the thickness of the layers can change the nonlinear vibration response of the structure significantly.

Influence of material positioning
To show the importance of the position of the hyperelastic material in a layered structure, a three-layered beam is assumed with the thickness of h = [h/4, h/2, h/4] which shows that the middle layer (the core) is thicker than the outer layers. In this section, by varying the materials in the layers, the nonlinear forced vibration response of the system is studied. To do this, the position of the material is varied and the dominant amplitudefrequency responses in axial, transverse and rotation motions of the structure are shown in Figs. 21 and 22 for pinned-pinned and clamped-clamped conditions, respectively. It can be seen that for both boundary conditions, the maximum amplitudes sweep to lower frequencies when the material positioning is [mat1-mat2-mat3]; the clamped-clamped and pinned-pinned beam models show that the dominant rotation coordinate reaches its highest when the material positioning is [mat1-mat2-mat3].

Influence of layer positioning
Another important goal of this study is to show the influence of layering and positioning the layers in manufacturing sandwich hyperelastic structures. To demonstrate this effect, the position of the layers is varied with

Summary and conclusions
In this study, a detailed layerwise analysis of the time-dependent mechanics of hyperelastic beams structures is given using five different shear deformation theories. The hyperelasticity was modelled following the Mooney-Rivlin hyperelastic strain energy density model, and the beam was assumed to lie on a foundation including linear, nonlinear and shear layers. Equations of motion were derived using Hamilton's principle together with different shear deformable beam theories. Solving the equations of motion using a dynamic equilibrium technique, it was shown that: • The classic beam theory predicts the fundamental frequency term of hyperelastic-layered beam models higher than the studied shear deformable beam models. • By increasing the length-to-thickness ratio, the influence of the type of boundary conditions decreases.
• Similar to elastic beams, the hyperelastic-layered beam model shows higher sensitivity to shear effects when the length-to-thickness ratio is lower. • The higher-order shear deformable beam models (namely third-order, exponential and trigonometric) give very close results for natural frequency terms, while the nonlinear behaviours are slightly different. (e) (f) Fig. 13 Influence of different higher-order shear effect models on the transverse amplitude response of the sandwich hyperelastic pinned-pinned beam a C 11 , b)C 31 , c C 12 , d C 32 , e C 13 and f C 33 • For linear natural frequencies, it was shown that having a foundation can increase the fundamental natural frequency of hyperelastic-layered beams significantly and this effect is higher for the shear layer compared with the linear layer. • For nonlinear forced time-dependent analysis of layered hyperelastic beams, the Timoshenko beam model gives the highest maximum amplitude for the dominant longitudinal, transverse and rotation deformations, compared with the other shear deformable models. Similarly, for higher-order shear models, the highest amplitudes are for the third-order beam model and the lowest one belongs to the exponential beam model. • For the two-layered hyperelastic beam model examined in this study, it is shown that the highest transverse amplitude belongs to the material combination of silicone and vulcanised rubber. Similarly, for the threelayered model, the material distribution silicone-thermoplastic-vulcanised rubber has the lowest maximum axial amplitude and the highest transverse and rotation amplitudes. • The influence of each layer's thickness is examined showing that changing the thickness of the layers can change the nonlinear vibration response of the structure significantly. • Hyperelastic material sorting has been shown to be a promising designing parameter in changing the dynamic behaviour of the structure for both linear and nonlinear analysis. For the studied model, both boundary conditions, the maximum amplitudes sweep to lower frequencies when the material positioning is [mat1-mat2-mat3].
Overall, since layered hyperelastic structures are replacing single-layer hyperelastic structures in many designs and applications such as food packaging, bottles, tubes and pipes, it is important to comprehend and predict the changes in their mechanical behaviour before employing them, especially their coupling motion behaviour. The results of this study are a step forward in realising the changes in the mechanical response of laminated hyperelastic beams when different materials and layers are used.