We develop an approach to reduce the governing equation of motion for the nonlinear vibration of a clamped laminated composite to the Duffing equation in a decoupled modal form. The method of weighted residuals enables such a reduction for laminates with clamped boundary conditions. Both rigidly clamped and loosely clamped boundary conditions are analyzed using this method. The reduction method conserves the total energy of the system. The decoupled modal form Duffing equation has constant modal parameters in terms of the laminated composite material's properties and geometries. The numerical computations illustrate the individual modal response with an emphasis of the transitional phenomena to chaos caused by the large load.

1. Introduction

The nonlinear dynamic response of a thin laminate occurs with a large curvature deformation. A common threshold for nonlinear phenomena is a deflection greater than half of the plate thickness. Large deformation associated with high-order strain can be described by a second-order von Karman strain field or a strain field of higher order. Nonlinear vibrations often arise when structures and membranes are subject to high pressure or intense thermal fields. This is typical of microstructures subject to aerodynamic forces, such as those as analyzed by Suhir [1]. Vibration can also be due to steady state thermal field variations, as studied by Pal [24]. In [5, 6], we discussed thermal vibrations of laminates caused by the transient thermal fluctuations.

One approach to analyzing the dynamic response of plates and laminates in nonlinear elastic deformation is to use a Galerkin-type method, or the energy method, to reduce the governing partial differential equation to an ordinary differential equation for the time-dependent variable. Such reduction procedures are complicated by the nonlinear strain compatibility equation. This makes it very difficult to obtain a decoupled modal form ordinary differential equation for the vibration of laminates; see Bloom and Coffin [7]. These reduction procedures often lead to a differential-integral equation with modal coupling, as detailed in [810]. Some authors have obtained a decoupled modal form ordinary differential equation for the laminate nonlinear vibration by ignoring the coupled modal terms, as done by Berger [11]. Berger’s approximation assumes that the second strain invariant of the mid-plane can be neglected when using the energy approach. Berger’s assumption compromises the total energy of the system. Nevertheless, Berger’s approximation has been widely applied to many nonlinear buckling and vibration problems. For example, Basuli [12] and Pal [24] studied the thermal buckling and thermal vibration of plates using Berger’s approximation. A nonlinear equation of the Duffing type is obtained for the plate buckling due to the static in-plane temperature variation. Pal also studied the thermal vibrations of both circular and rectangular plates due to spatial temperature variations for both simply supported and clamped boundary conditions. Pal’s reduction yielded a decoupled modal form Duffing equation by restricting the thermal field to be in the form of a polynomial [24]. In contrast, Boley and Weiner [13] developed an alternative approach to successively solve for the modal coefficients by using an iterative procedure. Their method transforms the modal coupling functions, which reduces to coupled modal coefficients in the Duffing equation, into successively modal-dependent coefficients that can be solved in a sequential manner for each mode. However, the method is applicable only to certain boundary conditions due to the validity of the approximations for the sequential solution.

Thermal vibration of functionally graded materials has been studied using a method that bears similar characteristics to that for the plate thermal vibrations. For example, Hao et al. [14] used a two-mode approximation for the nonlinear thermal vibration of bilayer materials in a simply supported boundary condition, with the restriction that the spatial thermal field is in a power law with respect to the thickness. A Galerkin’s approach led to two coupled modal equations for the plate vibration. The weighted residual terms associated with a thermal field, expressed either as a power law or as a polynomial, can be reduced to constants, similar to the procedure used by Pal [24]. It is only because of such restrictive thermal field variations that Pal and Hao were able to avoid coupled terms related to the thermal field, while the coupling of the two modal equations remains due to the use of a third order strain field with a two mode approximation. This two-mode approximation approach was also applied to the vibration of a cross-ply laminate due to the in-plane load by Zhang et al. [15]. The two mode coupling makes computation less intensive, compared to the coupling of infinite modes. However, our studies in [5, 6], indicate that modal coupling for the simply supported laminate can be eliminated by using a Galerkin method in conjunction with a second order von Karman strain field.

For the laminate vibration in clamped boundary conditions, two different forms of the deflection approximating function have been used previously, for both linear and nonlinear deformations. Sundara Raja Iyengar and Naqvi [9] and Young [16] used a combination of hyperbolic and harmonic functions, whereas Yamaki [17] and Berger [11] both used a product of second-degree harmonic functions. Pal [3] assumed a deflection function which is a variation of the function used by Yamaki [17]. Both types of approximating functions satisfy the clamped boundary conditions. However, the two approximating functions gave rise to Duffing equations in different modal forms. Sundara Raja Iyengar and Naqvi [9] and Young [16] obtained Duffing’s equation with complicated coefficients in terms of integrals of coupled Fourier series. These coefficients result from the orthogonality of the approximating function that cannot decouple all the members in the weighted residual. To avoid such complications, Yamaki [17] assumed a one-term truncation such that the reduction in Galerkin’s method led to a decoupled modal form Duffing equation with constant coefficients. The method of Yamaki [17] compromises the total energy of the system by eliminating coupled modal terms, essentially the same as Berger’s assumption does in terms of the energy conservation. These earlier studies suggest that the reduction from a partial differential equation in a multimode representation of deflection will produce modal coupling coefficients in the Duffing equation associated with clamped boundary conditions; only approximations that ignore the total energy conservation can lead to the decoupled form Duffing equation. We note that the Duffing equation in various forms is obtained as a direct consequence of using the second-order strain field for the nonlinear deformation.

Motivated by finding a unified approach to obtaining a decoupled modal form equation while preserving the total energy, we have found that the method of weighted residuals (MWR), as described by Crandall [18], Finlayson [19] and Stakgold [20], can be used for the laminate vibration in clamped boundary conditions. For the simply supported laminate, our study confirmed that Galerkin’s method is useful unless we have an arbitrary thermal field. This implies that for an arbitrary thermal field, we must use MWR as well for the simply supported boundary condition [5]. The MWR is also referred to as the generalized Galerkin method, as discussed by Finlayson in [19], where the weighting function need not be the approximating function. That is, the MWR allows us to use a weighting function that is different from that used in the standard Galerkin method.

In this paper, the MWR is applied to study the deflection of symmetric thin isotropic laminates for both rigidly clamped and loosely clamped boundary conditions. This method leads to a decoupled modal form Duffing equation for each mode, with constant coefficients expressed in terms of the material properties and laminate geometries. However, our reduction approach requires a particular choice of the approximating function, as we will demonstrate. In the following, we first obtain the Duffing equation in a decoupled modal form. Our key interest in the Duffing system here lies in the transitional phenomena associated with the loading and initial conditions, instead of parametric analyses as we did earlier [5, 6]. To this end, we use a Lyapunov function to justify the weak chaotic attractors of the stable Duffing equation, and identify the stability condition in the time domain for the transition from a stable to an unstable Duffing system due to a large load. Such transitional behaviors are illustrated by the fundamental mode responses from numerical computations. We verified this phenomenon by computing the Lyapunov exponents, which indicated that the transition has sensitive dependence on the initial and loading conditions. Finally, the paper closes with discussions and conclusions.

2. Equation of Motion

The nonlinear von-Karman strain field is defined as 𝜀{𝜀}=0𝜀+𝑧1,(2.1)

𝜀0=𝜀0𝑥𝑥𝜀0𝑦𝑦𝛾0𝑥𝑦,𝜀1=𝜅𝑥𝜅𝑦𝜅𝑥𝑦𝜀,(2.2a)0𝑥𝑥=𝑢0,𝑥+12𝑤,𝑥2,𝜀0𝑦𝑦=𝑣0,𝑦+12𝑤,𝑦2,𝛾0𝑥𝑦=𝑢0,𝑦+𝑣0,𝑥+𝑤,𝑥𝑤,𝑦,𝜅𝑥=𝑤,𝑥𝑥,𝜅𝑦=𝑤,𝑦𝑦,𝜅𝑥𝑦=2𝑤,𝑥𝑦.(2.2) The strain field satisfies the compatibility equation:

𝜀0𝑥𝑥,𝑦𝑦+𝜀0𝑦𝑦,𝑥𝑥𝛾0𝑥𝑦,𝑥𝑦=𝑤,2𝑥𝑦𝑤,𝑥𝑥𝑤,𝑦𝑦.(2.3a) Define an operator

𝐿(𝑝,𝑞)=𝑝,𝑦𝑦𝑞,𝑥𝑥+𝑝,𝑥𝑥𝑞,𝑦𝑦2𝑝,𝑥𝑦𝑞,𝑥𝑦,(2.3b) equation (2.3a) can be expressed as

𝜀0𝑥𝑥,𝑦𝑦+𝜀0𝑦𝑦,𝑥𝑥𝛾0𝑥𝑦,𝑥𝑦1=2𝐿(𝑤,𝑤).(2.3c) The constitutive relation between the forces and the strains is:

𝜀0𝑀=𝐴𝐵𝐵𝑇𝐷𝑁𝜀1,(2.4) where: 𝐴=𝐴1,𝐵=𝐴1𝐵,𝐷=𝐷𝐵𝐴1𝐵.(2.5a) The stiffness matrices of the laminate are defined as:

𝐴𝑖𝑗,𝐵𝑖𝑗,𝐷𝑖𝑗=21𝐶(𝑘)𝑖𝑗1,𝑧,𝑧2𝐶𝑑𝑧,(2.5b)(𝑘)11=𝐶(𝑘)22=𝑘=𝑛𝑘=1𝐸(𝑘)𝜈1(𝑘)2,𝐶(𝑘)12=𝐶(𝑘)21=𝑘=𝑛𝑘=1𝜈(𝑘)𝐸(𝑘)𝜈1(𝑘)2,𝐶(𝑘)66=𝑘=𝑛𝑘=1𝐸(𝑘)21+𝜈(𝑘),(2.5c) where 𝑛 is the number of the layers in symmetry. For a symmetric isotropic laminate, [𝐵]=0, 𝐷=𝐷. In (2.4), 𝑁 and 𝑀 are the in-plane force and moment vectors, respectively. Introducing the Airy stress function defined by 𝑁𝑥=𝐹,𝑦𝑦,𝑁𝑦=𝐹,𝑥𝑥,𝑁𝑥𝑦=𝐹,𝑥𝑦,(2.6) from (2.4), the compatibility (2.3c) can be expressed as

𝐹,𝑥𝑥𝑥𝑥+𝑝𝐹,𝑥𝑥𝑦𝑦+𝐹,𝑦𝑦𝑦𝑦𝑟=2𝐿(𝑤,𝑤),(2.7) where a symmetric isotropic laminate has 𝑟=1𝐴22,𝑝=2𝐴12+𝐴66𝐴22=2.(2.8) Therefore, for the symmetric isotropic laminate, (2.7) can be further reduced to: 4𝑟𝐹=2𝐿(𝑤,𝑤).(2.9) The governing equations of motion for a laminate in a nonlinear elastic deformation are given in [21]: 𝑁𝑥,𝑥+𝑁𝑥𝑦,𝑦𝑁=0,𝑥𝑦,𝑥+𝑁𝑦,𝑦𝑀=0,𝑥,𝑥𝑥+2𝑀𝑥𝑦,𝑥𝑦+𝑀𝑦,𝑦𝑦+𝑁(𝑤)+𝑞=𝐼𝑤,𝑡𝑡𝐼22𝑤,𝑡𝑡,(2.10) where 𝑁(𝑤)=𝑁𝑥𝑤,𝑥𝑥+2𝑁𝑥𝑦𝑤,𝑥𝑦+𝑁𝑦𝑤,𝑦𝑦.(2.11) Define the inertia of the laminate as 𝐼,𝐼1,𝐼2=𝑁𝑘=121𝜌0(𝑘)1,𝑧,𝑧2𝑑𝑧.(2.12) A homogeneous material with a symmetric layer distribution yields 𝐼1=0. Using (2.4), (2.6) and (2.10), we obtain the equilibrium equation for the deflection of a symmetric isotropic laminate: 𝐷114𝑤+𝐼𝑤,𝑡𝑡𝐼22𝑤,𝑡𝑡𝐿(𝐹,𝑤)=𝑞,(2.13) where: 𝐿(𝐹,𝑤)=𝐼𝐽𝑖=1𝑁𝑗=1𝑀𝑛=1𝑚=1𝐹𝑚𝑛,𝑥𝑥𝑤𝑖𝑗,𝑦𝑦+𝐹𝑚𝑛,𝑦𝑦𝑤𝑖𝑗,𝑥𝑥2𝐹𝑚𝑛,𝑥𝑦𝑤𝑖𝑗,𝑥𝑦.(2.14) Here 𝐹𝑚𝑛 are the modal functions of the Airy stress function 𝐹(𝑥,𝑦,𝑡), and 𝑤𝑖𝑗 are the modal components of the deflection function 𝑤(𝑥,𝑦,𝑡). For either the rigidly clamped or the loosely clamped laminates, 𝑤(𝑥,𝑦,𝑡) is chosen as the deflection approximating function in series as

𝑤(𝑥,𝑦,𝑡)=𝑁𝑀𝑛=1𝑚=1𝑊𝑚𝑛(𝑡)𝑤(𝑥,𝑦)=𝑁𝑀𝑛=1𝑚=1𝑊𝑚𝑛(𝑡)cos2𝛼𝑚𝑥cos2𝛽𝑛𝑦,𝑎𝑥2,𝑎2𝑏,𝑦2,𝑏2,𝛼𝑚=𝑚𝜋𝑎,𝛽𝑛=𝑛𝜋𝑏,𝑚=2𝑘1,𝑘=1,2,𝑀.𝑛=2𝑗1,𝑗=1,2,𝑁.(2.15) This expression satisfies the actual boundary conditions, namely:

𝑤𝑏𝑥,2𝑤𝑏,𝑡=0,𝑥,2𝑤𝑎,𝑡=0,2𝑤𝑎,𝑦,𝑡=0,2,𝑦,𝑡=0,𝑤,𝑦𝑏𝑥,2,𝑡=0,𝑤,𝑦𝑏𝑥,2,𝑡=0,𝑤,𝑥𝑎2,𝑦,𝑡=0,𝑤,𝑥𝑎2,𝑦,𝑡=0.(2.16a) In addition, the rigidly clamped laminate must satisfy the in-plane motion constraints of the mid-plane at the boundaries as:


Note that the approximating functions in (2.15) are the same as those previously used by Yamaki [17]; however, our procedure for reduction is not the same as in [17] using Galerkin’s method. Instead, MWR, as described in [1820], will be used. In particular, a weighting function that is a derivative of the set of the approximating function is used for modal decoupling.

Define the residual of the equilibrium equation, (2.13), due to the approximating function 𝑤(𝑥,𝑦,𝑡), as

𝐸=𝐷114𝑤𝑤+𝐼,𝑡𝑡𝐼22𝑤,𝑡𝑡𝑤𝐿𝐹,𝑞.(2.17) Using MWR, we require that

𝑏/2𝑏/2𝑎/2𝑎/2𝐸𝜙(𝑥,𝑦)𝑑𝑥𝑑𝑦=0,(2.18) where the weighting function has the same index as that of the approximating function in (2.15), given as:

𝑤𝜙(𝑥,𝑦),𝑥𝑥𝑦𝑦(𝑥,𝑦)=𝑁𝑀𝑛=1𝑚=1cos2𝛼𝑚𝑥cos2𝛽𝑛𝑦.(2.19) This makes (2.18) equivalent to

𝑏/2𝑏/2𝑎/2𝑎/2𝐷114𝑤𝑤+𝐼,𝑡𝑡𝐼22𝑤,𝑡𝑡𝑤𝐿𝐹,𝑞cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦𝑑𝑥𝑑𝑦=0.(2.20) As a preliminary to evaluate integral (2.20), various calculations are necessary. Using the deflection approximation 𝑤(𝑥,𝑦,𝑡) in (2.15), the function in (2.9) becomes:

4𝐹𝑚𝑛𝑟=2𝛼2𝑚𝛽2𝑛𝑊2𝑚𝑛cos2𝛼𝑚𝑥+cos2𝛽𝑛𝑦+𝑐𝑜𝑠4𝛼𝑚𝑥+cos4𝛽𝑛𝑦+cos4𝛼𝑚𝑥cos2𝛼𝑚𝑥+cos4𝛽𝑛𝑦+2cos2𝛼𝑚𝑥cos2𝛽𝑛𝑦.(2.21) Solving (2.21), we obtain

𝐹𝑚𝑛=12𝐶1𝑚𝑛𝑦2+12𝐶2𝑚𝑛𝑥2+𝐹𝑚𝑛,(2.22) where

𝐹𝑚𝑛=𝑊2𝑚𝑛𝑓(𝑡)1cos2𝛼𝑚𝑥+𝑓2cos2𝛽𝑛𝑦+𝑓3cos4𝛼𝑚𝑥+𝑓4cos4𝛽𝑛𝑦+𝑓5cos2𝛼𝑚𝑥cos4𝛽𝑛𝑦+𝑓6cos4𝛼𝑚𝑥cos2𝛽𝑛𝑦+𝑓7cos2𝛼𝑚𝑥cos2𝛽𝑛𝑦(2.23a) and the coefficients 𝑓𝑖 are given by

𝑓1𝑟=𝛽2𝑛32𝛼2𝑚,𝑓2𝑟=𝛼2𝑚32𝛽2𝑛,𝑓3𝑟=𝛽2𝑛512𝛼2𝑚,𝑓4𝑟=𝛼2𝑚512𝛽2𝑛,𝑓5𝑟=𝛼2𝑚𝛽2𝑛𝛼322𝑚+4𝛽2𝑛2,𝑓6𝑟=𝛼2𝑚𝛽2𝑛324𝛼2𝑚+𝛽2𝑛2,𝑓7𝑟=𝛼2𝑚𝛽2𝑛𝛼162𝑚+𝛽2𝑛2.(2.23b) In order to satisfy the in-plane motion constraints (2.16b) for a rigidly clamped laminate, 𝐶1𝑚𝑛=3𝑊𝑚𝑛(𝑡)2𝐴3211𝛼2𝑚+𝐴12𝛽2𝑛,𝐶2𝑚𝑛=3𝑊𝑚𝑛(𝑡)2𝐴3212𝛼2𝑚+𝐴11𝛽2𝑛,(2.23c) so that for the rigidly clamped laminate

𝐹𝑚𝑛=3𝑊𝑚𝑛(𝑡)2𝐴6412𝛼2𝑚+𝐴11𝛽2𝑛𝑥2+𝐴11𝛼2𝑚+𝐴12𝛽2𝑛𝑦2+𝐹𝑚𝑛,(2.24) where 𝐹𝑚𝑛 is as expressed in (2.23a). For the loosely clamped laminate, which permits the in-plane motion at the boundaries, it is sufficient to take 𝐶1𝑚𝑛=𝐶2𝑚𝑛=0, so that

𝐹𝑚𝑛=𝐹𝑚𝑛.(2.25) To evaluate the integral (2.20), we introduce the following symmetric pairs for integration.

𝐽𝐼={𝐽1𝐼(𝛼𝑚𝑥,𝛽𝑛𝑦)𝐽1𝐼(𝛼𝑝𝑥,𝛽𝑞𝑦),𝐽2𝐼(𝛼𝑚𝑥,𝛽𝑛𝑦)𝐽2𝐼(𝛼𝑝𝑥,𝛽𝑞𝑦)}𝐼=0,1,...7. (19) These pairs are defined by:

𝐽0=cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦,cos2𝛽𝑞𝑦cos2𝛼𝑝𝑥,𝐽1=cos2𝛼𝑚𝑥cos2𝛽𝑞𝑦cos2𝛼𝑝𝑥,cos2𝛽𝑛𝑦cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦,𝐽2=cos4𝛼𝑚𝑥cos2𝛽𝑞𝑦cos2𝛼𝑝𝑥,cos4𝛽𝑛𝑦cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦,𝐽3=cos2𝛼𝑚𝑥cos4𝛽𝑛𝑦cos2𝛽𝑞𝑦cos2𝛼𝑝𝑥,cos4𝛼𝑚𝑥cos2𝛽𝑛𝑦cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦,𝐽4=cos4𝛼𝑚𝑥cos2𝛽𝑛𝑦cos2𝛽𝑞𝑦cos2𝛼𝑝𝑥,cos2𝛼𝑚𝑥cos4𝛽𝑛𝑦cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦,𝐽5=cos2𝛼𝑚𝑥cos2𝛽𝑛𝑦cos2𝛽𝑞𝑦cos2𝛼𝑝𝑥,cos2𝛼𝑚𝑥cos2𝛽𝑛𝑦cos2𝛼𝑝𝑥cos2𝛽𝑞𝑦,𝐽6=sin2𝛼𝑚𝑥sin4𝛽𝑛𝑦sin2𝛽𝑞𝑦sin2𝛼𝑝𝑥,sin4𝛼𝑚𝑥sin2𝛽𝑛𝑦sin2𝛼𝑝𝑥sin2𝛽𝑞𝑦,𝐽7=sin2𝛼𝑚𝑥sin2𝛽𝑛𝑦sin2𝛼𝑝𝑥sin2𝛽𝑞𝑦.(2.26) Each of the above pairs contains a product of a function of x and a function of 𝑦, that is,

𝐽𝑘𝐼(𝑥,𝑦)=𝐽𝑘𝐼𝑥𝑓1𝛼𝑚𝑥,𝑓2𝛼𝑝𝑥𝐽𝑘𝐼𝑦𝑔1𝛽𝑛𝑦,𝑔2𝛽𝑞𝑦𝑘=1,2;𝐼=0,1,7.(2.27) In (2.27), not all functions 𝑓𝑖 and 𝑔𝑗 need to be present in each pair. Noticeably, the two functions 𝐽1𝐼 and 𝐽2𝐼 in each pair are symmetric with respect to their dependent variables, 𝑥and 𝑦, that is,

𝐽1𝐼𝑥𝑓1𝛼𝑚𝑥,𝑓2𝛼𝑝𝑥𝐽1𝐼𝑦𝑔1𝛽𝑛𝑦,𝑔2𝛽𝑞𝑦=𝐽2𝐼𝑥𝑔1𝛼𝑚𝑥,𝑔2𝛼𝑝𝑥𝐽2𝐼𝑦𝑓1𝛽𝑛𝑦,𝑓2𝛽𝑞𝑦𝐼=0,1,7.(2.28a) This means by switching the argument (𝛼𝑚𝑥),(𝛼𝑝𝑥) with (𝛽𝑛𝑦),(𝛽𝑞𝑦), respectively, in the corresponding function, the two components in each pair is the same, equivalently;

𝐽1𝐼𝑓𝛼𝑚𝑥,𝛼𝑝𝑥𝛽,𝑔𝑛𝑦,𝛽𝑞𝑦=𝐽2𝐼𝑓𝛽𝑛𝑦,𝛽𝑞𝑦𝛼,𝑔𝑚𝑥,𝛼𝑝𝑥𝐼=0,1,7.(2.28b) Therefore, the evaluation of each integral pair denoted by 𝑘=1,2, respectively, in the spatial coordinates has the following properties: 𝑏/2𝑏/2𝑎/2𝑎/2𝐽𝑘𝐼(𝑥,𝑦)𝑑𝑥𝑑𝑦=𝑎/2𝑎/2𝐽𝑘𝐼𝑥(𝑥)𝑑𝑥𝑏/2𝑏/2𝐽𝑘𝐼𝑦(𝑦)𝑑𝑦,𝑘=1,2,𝐼=0,1,7,𝑏/2𝑏/2𝑎/2𝑎/2𝐽1𝐼𝑑𝑥𝑑𝑦=𝑏/2𝑏/2𝑎/2𝑎/2𝐽2𝐼𝑑𝑥𝑑𝑦,𝑏/2𝑏/2𝑎/2𝑎/2𝐽1𝐼𝑤,𝑥𝑥𝑦𝑦𝑑𝑥𝑑𝑦=𝑏/2𝑏/2𝑎/2𝑎/2𝐽2𝐼𝑤,𝑥𝑥𝑦𝑦𝑑𝑥𝑑𝑦.(2.29) By using the orthogonal properties between each member in (2.17) and the weighting function and the symmetry properties in (2.29), all of the modal functions in (2.20) are decoupled, giving rise to constant coefficients in the Duffing equation. Therefore, MWR reduces (2.20) to the decoupled modal form Duffing equation for the transient response of each mode as

̈𝑊𝑚𝑛+𝜔2𝑚𝑛𝑊𝑚𝑛+𝑟𝑚𝑛𝑊3𝑚𝑛=𝑞𝑚𝑛,(2.30) where the coefficients for an isotropic laminate in a rigidly clamped boundary condition are

𝜔2𝑚𝑛=16𝐷11𝛼2𝑚+𝛽2𝑛2𝐼,𝑟(2.31a)𝑚𝑛=1𝐼6𝐴1611𝛼4𝑚+2𝐴12𝛼2𝑚𝛽2𝑛+𝐴11𝛽4𝑛+𝐴211𝐴212𝐴119𝛼4𝑚+𝛽4𝑛+4𝛼4𝑚𝛽4𝑛𝛼2𝑚+4𝛽2𝑛2+4𝛼4𝑚𝛽4𝑛4𝛼2𝑚+𝛽2𝑛2,𝑞(2.31b)𝑚𝑛=𝑄𝑚𝑛𝐼,𝛼(2.31c)𝐼=𝐼+2𝑚+𝛽2𝑛𝐼2.(2.31d) For a loosely clamped laminate, we obtain an identical Duffing equation with the same 𝜔2𝑚𝑛 and 𝑞𝑚𝑛 expressions as above, except that the stiffness coefficient is given by

𝑟𝐿𝑚𝑛=1𝐼𝐴16211𝐴212𝐴119𝛼4𝑚+𝛽4𝑛+4𝛼4𝑚𝛽4𝑛𝛼2𝑚+4𝛽2𝑛2+4𝛼4𝑚𝛽4𝑛4𝛼2𝑚+𝛽2𝑛2.(2.31e) The Duffing equation (2.30) represents a hard spring system for both the rigidly clamped and the loosely clamped boundaries, due to the relation between the laminate stiffness coefficients of 𝐴11>𝐴12 and 𝐴𝑖𝑗>0,which yields 𝑟𝑚𝑛>0 in (2.31b) and 𝑟𝐿𝑚𝑛>0 in (2.31e), respectively. The loosely clamped boundary condition has a reduced stiffness, that is, 𝑟𝐿𝑚𝑛<𝑟𝑚𝑛. Equations (2.31a), (2.31b), and (2.31e) indicate that both the natural frequency 𝜔𝑚𝑛 and the stiffness 𝑟𝑚𝑛 and 𝑟𝐿𝑚𝑛 decrease with an increase of the inertia 𝐼, due to the increased density and thickness. On the contrary, a higher rigidity 𝐷𝑖𝑗 and stiffness 𝐴𝑖𝑗 of the laminate, owing to an increased thickness, increases the natural frequency and stiffness. In addition, an increase in the geometric dimensions of the laminate plane decreases the natural frequency and stiffness as a result of reduced modal parameter 𝛼𝑚𝑜𝑟𝛽𝑛. In comparison, a higher mode has a higher 𝜔𝑚𝑛 and 𝑟𝑚𝑛, or 𝑟𝐿𝑚𝑛, since 𝑎𝑚𝑛 and 𝛽𝑚𝑛 are amplified by the modal integers 𝑚 and 𝑛.

3. The Chaotic Attractor and Stability Transitions

The general characteristics of the Duffing equation have been studied extensively for various behaviors that are relevant for laminate response analysis. The Duffing equation with damping in a normalized form is

𝑑2𝑊𝑑𝜏2̇+𝑊+̃𝜂𝑊+̃𝑟𝑊3=𝑓cos(𝑘𝜏),(3.1) where 𝜏=𝜔𝑡, ̃𝑟=𝑟/𝜔2,̃𝜂=𝜂/𝜔2,𝑓=𝑓/𝜔2,𝑘=Ω/𝜔, and 𝜂 is the damping coefficient. Closed form solutions in terms of elliptical functions were found in [1, 22, 23] for the Duffing equation without damping, in the free response and the response subject to a constant load. Marinca and Herisanu [24] and Tang [25] demonstrated that the periodic solutions can be obtained by numerical computation for a nonlinear system subject to harmonic excitation. Approximate solutions were found for the finite periodic orbits of the Duffing equation by perturbation of the averaged system, as described by Nayfeh and Balachandran, Guckenheimer and Holmes in [2628], and by the method of harmonic balancing by Hayashi [29]. It was shown in [29, 30] that the nonlinear resonance frequency and amplitude depend on the load. Studies in [28, 31, 32] demonstrate that chaotic attractors of the Duffing system have a sensitive dependence on the initial and loading conditions. However, these studies were concerned with a Duffing system with an unstable separatrix where 𝜔2<0, or a system with a negative stiffness 𝑟<0. Consequently, the results may not directly apply to the Duffing equation with 𝜔2>0 and a positive stiffness 𝑟>0 and without damping, as in (2.30).

For our Duffing equation (2.30), it can be shown that the Mel’nikov function for this system is identically zero due to the homoclinic orbit around the center. This means that the Mel’nikov function is not adaptable to establish any bifurcation condition with respect to the loading for the onset of chaos, which is in contrast to the case of 𝜔2<0, as studied in [33]. To describe the behavior of our system subject to the large load, we found that the weak Lyapunov function of the Duffing equation can serve as an indicator of a weak attractor, with a tendency toward chaos when the system is subjected to perturbations. Denoting ̇𝑉=𝑊=𝑑𝑊/𝑑𝑡, the weak Lyapunov function of the system is

𝑉𝐿(𝑊,𝑉)=22+𝑊22+𝑟𝑊42>0.(3.2) The function 𝐿(𝑊,𝑉) is positive over the plane, except at the center 𝐿(𝑊,𝑉)=𝐿(0,0)=0. The stability of the system can be measured by the time derivative of the Lyapunov function, which vanishes at the center:

̇𝐿(0,0)=𝑊+𝑟𝑊3̇̇𝑊+𝑉𝑉=0.(3.3) This means that the center is a weak attractor and is Lyapunov stable [34]. This system becomes chaotic when subjected to certain levels of perturbations, such as a large dynamic load as we demonstrate hereinafter. The time elapsed can be estimated before the system changes from a stable system to an unstable system. To this end, we express the harmonic forcing 𝑓cos(𝑘𝜏) in the form of a series, so that the Duffing equation (3.1), without damping becomes

̈𝑊+𝑊+̃𝑟𝑊3=𝑓1(𝑘𝜏)2+2!(𝑘𝜏)44!(𝑘𝜏)6+6!(𝑘𝜏)88!++(1)𝑛(𝑘𝜏)2𝑛=(2𝑛)!𝑓𝑔(𝑘𝜏),𝑛=0,1(3.4) The series 𝑔(𝑘𝜏) in expression (3.4) satisfies |𝑔(𝑘𝜏)|1. Assuming the load is in proportion to the deflection, that is, 𝑓=𝜇𝑊, we obtain after rearranging terms

̈𝑊+𝜔2𝑊+̃𝑟𝑊3=0,(3.5a) where

𝜔2=1𝜇1(𝑘𝜏)22+(𝑘𝜏)44!(𝑘𝜏)6+6!(𝑘𝜏)88!+=1𝜇𝑔(𝑘𝜏).(3.5b) Here the factor 𝜇1, 𝜇1 and 𝜇1 represent the large load, the small load and the intermediate load regimes, respectively. Equation (3.5b) implies that when 𝜇1, the response of the system is almost periodic since 𝜔21 for all values of 𝑘𝜏. When 𝜇1 or 𝜇1, 𝜔2 in (3.5b) can become negative. The transition time 𝜏𝑐, when the system passes out of the stable system at 𝜔2=0 can be obtained numerically as infinite number of roots of (3.5b), which corresponds to infinite transitions from a stable to an unstable system. An approximation of the first passage time 𝜏𝑐 when 𝑘𝜏<1 with 𝜇1, based on retaining only three terms in (3.5b) for 𝜔2<0, is:

𝜏c=max12k,1k6+263+𝜇.(3.6) This threshold time 𝜏𝑐 signifies the system’s transition into an unstable system with 𝜔2<0. This unstable system becomes chaotic if the load is very large, or the initial deflection and the initial velocity are large, or a combination of these two situations, as established in [27, 28].

4. Numerical Computation

We will use the fourth-order Runge-Kutta method to illustrate the fundamental mode response of the Duffing equation for the chaotic and periodic behaviors in relation to the loading and initial conditions. The computation also generates the Lyapunov exponents to verify the transition to chaos. The case study examines a thin laminate of 1301801.54mm in symmetric layers as shown in Figure 1. The physical properties of the laminate are listed in Table 1. The parameters for the second, third and fourth modes of the Duffing equation are shown in Table 2. Table 2 indicates that a higher mode has higher stiffness and higher natural frequency, which give rise to reduced deflections when subjected to the same loading compared to the lower mode responses. For example, our computations show that the deflection is in the ratio of 1/3, 1/4, and 1/6 of that of the fundamental mode response for the seconh, third and fourth modes under 𝑄11=12cos(𝑘𝜔11𝑡)N/cm2. The computation specifies the pressure load to the laminate 𝑄𝑚𝑛, instead of the normalized load 𝑞𝑚𝑛 in (2.31c). The conversion between the pressure load and the normalized load is 𝑞11/𝑄11=0.5215 for the fundamental mode, Mode-1, whereas the normalized force 𝑞 and the pressure load 𝑄 are measured in [mm] and [N/cm2], respectively.

Figures 2(a) and 2(b) represent the same Poincare map of Mode-1 in a quasiperiodic behavior with finite periodic orbits subject to, 𝑄11=2.4cos(𝑘𝜔11𝑡)N/cm2, 𝑘=1, with initial condition 𝑥0=0, and𝑣0=0. A periodic pattern is visible by the consecutive traces of each period, as shown in Figure 2(a). The periodic behavior remains with a moderate increase in the load and the initial conditions, as shown in Figure 3(a) where a finite period orbit emerges subject to 𝑥0=2mm and a load 𝑄11=24cos(𝑘𝜔11𝑡)N/cm2,𝑘=1. A transition to chaos arises from a higher load or large initial conditions. Figure 3(b) shows the transition to chaos when the initial condition increases to 𝑥0=20mm while the load remains at 𝑄11=24cos(𝑘𝜔11𝑡)N/cm2,𝑘=1. The result confirms the center of the phase map as a weak attractor, as predicted by the Lyapunov function as analyzed above. Around the center, there is a limit cycle. The attractor has sensitive dependence on initial conditions, namely, a chaotic attractor.

Figure 3(c) shows the Mode-1 Poincare map in chaos subject to a large load of 𝑄11=432cos(𝑘𝜔11𝑡)N/cm2,𝑘=1, whereas the initial 100 transient points are eliminated. The map contains both chaotic and periodic orbits, whereas the heteroclinic manifold for the unstable system has two attractors, and the stable system has one attractor at the center. The result indicates a coexistence of chaotic and periodic orbits [31]. The coexistence of the two types of orbits persists, regardless of the transient time eliminated from the map, which suggests that the large load causes alternating transitions of the system between 𝜔2>0 and 𝜔2<0, as we analyzed as we analyzed the roots of (3.5b) in the afore discussed section. Generally speaking, a two-attractor map symbolizes the heteroclinic orbit of an unstable Duffing system with 𝜔2𝑚𝑛<0, which is in contrast to the homoclinic orbit around the center of a stable Duffing system with 𝜔2𝑚𝑛>0 [27]. This means that a larger load effectively transforms a stable system to an unstable system. Figure 3(d) shows the response subject to a further increasing load of 𝑄11=4800cos(𝑘𝜔11𝑡)N/cm2,𝑘=1 without the initial 100 transient periods. Evidently, the periodic orbits with two attractors are in chaos, identical to that shown in Figure 3(c).

Other than the Poincare map, the Lyapunov exponent provides an alternative indicator for chaos. The Lyapunov exponents for the Duffing system were computed using the algorithms developed for low-dimensional systems by Dieci [35], which converges for 10 000 to 30 000 iterations. It is found that the largest Lyapunov exponent is negative for the case in Figure 3(a), which means that it is a periodic response. The largest Lyapunov exponent for the 2D Duffing system of the differential equation converges to 0.00283, 0.36466, and 1.29960 for the responses in Figures 3(b), 3(c), and 3(d), which indicate that each response is chaotic. In addition, the sum of the Lyapunov exponents is negative for each case, which suggests that the chaos is bounded [36], as indicated by the Poincare map. Noticeable also is that a larger deflection prevails at chaos compared with that of a periodic response.

These results point out that a large amplitude harmonic forcing can transform a stable Duffing system into an unstable system. In addition, chaos shows a sensitive dependence on the initial and loading conditions, as observed from the comparisons of Figures 3(a), 3(b), and 3(c). Our results suggest that the transitional phenomenon of the Duffing system without damping, that has sensitive dependence on the initial conditions, is similar to that of a system with damping as studied earlier by Ueda, Holmes, and Guckenheimer [28, 31, 32].

5. Discussion

The MWR method used in the present study, for obtaining the Duffing equation, eliminated the modal coupling resulting from a use of the standard Galerkin method. The modal decoupling avoids the computational difficulties associated with coupled modal coefficients in integral forms, such as those obtained by Sundara Raja Iyengar and Naqvi [8, 9]. A distinction between the present method and the method used by Berger [11] and Yamaki [17] is that their methods do not conserve the total energy of the system. Specifically, in Yamaki’s analysis using the Galerkin method, the coupling term is eliminated by imposing the condition 0𝑏/20𝑎/2(𝐿(𝐹,𝑤)𝑤)𝑑𝑥𝑑𝑦=0. Similarly, in Berger’s approach using the energy method for the single mode response, the total energy is altered by neglecting the second strain invariant in the strain energy expression:

𝐷𝑈=20𝑏/20𝑎/22𝑤+122𝑒122(1𝜈)122𝑒2+𝜕2𝑤𝜕𝑥2𝜕2𝑤𝜕𝑦2𝜕2𝑤𝜕𝑥𝜕𝑦2𝑑𝑥𝑑𝑦,(5.1) where the first and the second strain invariants are defined as

𝑒1=𝜀𝑥+𝜀𝑦,𝑒2=𝜀𝑥𝜀𝑦14𝛾2𝑥𝑦.(5.2) Berger’s method applied the principle of virtual work with the assumption 𝑒2=0, in order to obtain the decouple modal form Duffing equation. In contrast to Berger and by Yamaki, the Duffing equation obtained using MWR retains the total energy. Consequently, the results from our method are more accurate, in the sense of Galerkin type averaging.

It has been discussed by Crandall et al. [1820] that in using the MWR there is flexibility in choice of the approximating and weighting functions. For example, Shen and Pesheck used a weighting function that is different from their approximating function, in discussing ordinary and partial differential equations [37, 38]. Our choice of 𝑤,𝑥𝑥𝑦𝑦(𝑥,𝑦) for the weighting function, where 𝑤(𝑥,𝑦) is our approximating function, enables modal decoupling in all the weighted residual terms in (2.20), while retaining all these terms. This leads directly to the Duffing equation with constant coefficients. The modal form Duffing equation can readily be applied to each individual mode response analysis, such as resonance, periodic, and chaotic behaviors with respect to various parameters, as we studied in [39] for the clamped laminate with thermal effect. An accurate account for each modal response provides a basis for truncations in a multimode superposition and for identification of internal resonance among different modes. This is of particular interest when the response is chaotic, since an increased deflection at chaos of higher modes becomes much more significant in the total response of multimodes.

Our decoupling method would not be valid for the approximating function for the clamped laminate as that used by Sundara Raja Iyengar and Naqvi [9] and Young [16], that is, 𝑋𝑚𝛼𝑚𝑥𝑋𝑛𝛽𝑛𝑦=𝛼cosh𝑚𝑥𝛼cosh𝑚𝑎𝛼cos𝑚𝑥𝛼cos𝑚𝑎𝛽cosh𝑛𝑦𝛽cosh𝑛𝑏𝛽cos𝑛𝑦𝛽cos𝑛𝑏.(5.3) It can be verified that the reduction using the standard Galerkin method with this approximating function, and the weighting function is identically equal to the approximating function, leads to multimode integral terms that cannot be reduced to constants. For example, the following integral function remains:

𝑎𝑎cosh2𝛼𝑚𝑥𝛼cosh𝑗𝑥𝛼cosh𝑘𝑥𝑑𝑥.(5.4) It is apparent that orthogonality of this approximating function determines that any other weighting function will also give rise to integral coefficients. Therefore, neither the standard Galerkin method nor MWR will lead to a decoupled modal form Duffing equation using the function in (5.3).

We note that our reduction procedure resulted in a decoupled modal form Duffing equation. This differs from the reduction for the normal forms of a nonlinear vibration system, as studied by Shaw and Pierre [40] as well as by Pesheck et al. [38]. Specifically, our objective is to obtain an equation for the time-dependent deflection functions for the laminate dynamic response, whereas the normal form reduction schemes focus on finding the spatial functions governing the normal manifolds for a specific equation of motion. This suggests that the search for the normal forms for each decoupled modal form Duffing equation, as obtained from our study, can be conducted using the procedures outlined by Pesheck and Shaw et al. An interesting observation from these two types of drastically different studies is that Pesheck essentially used the method of weighted residuals in obtaining the normal forms. This indicates the usefulness of the MWR in solving various problems.

6. Conclusions

A decoupled modal form Duffing equation with constant coefficients is obtained for the nonlinear vibration of a thin symmetric isotropic laminate for both the rigidly clamped and loosely clamped boundary conditions. This equation results from a reduction of the governing equation of motion by using the method of weighted residuals, with a total conservation of energy in the sense of Galerkin’s averaging. The benefit, other than the modal decoupling, is an improved accuracy and computation efficiency. The decoupled modal form Duffing equation provides a direct formulation for the laminates’ modal response. Our analytical and computation studies identified the transitional phenomena of the Duffing equation from a stable system to an unstable system subject to a large load. The weak Lyapunov function justifies the weak attractors of the stable system as a path to chaos. The results also illustrate that chaos has a sensitive dependence on the loading and the initial conditions. By generalizing the Galerkin approach to the MWR, the present study found a particular procedure for modal decoupling of the Duffing equation for the dynamics of laminates with clamped boundary conditions.


𝐴:Laminate rigidity matrix of the laminate
𝐴:Inverse of A
𝐴𝑖𝑗:Elements of 𝐴
𝐵:Coupling rigidity matrix of the laminate
𝐶:Stiffness matrix of the laminate
𝐶(𝑘)𝑖𝑗,𝐶(𝑘)𝑖𝑗,𝑖,𝑗=1,2,6:Stiffness coefficients of the laminate
𝐶1𝑚𝑛,𝐶2𝑚𝑛:Complementary function coefficients of the Airy stress function
𝐷:Flexural rigidity matrix of the laminate
𝐷𝑖𝑗: Matrix components of 𝐷
𝐸:Residual of the equation of motion
𝐸(𝑘):𝑘th lamina Young's Modulus
𝑓: Normalized load
𝑓𝑖:Airy’s stress function coefficients
𝐹(𝑥,𝑦,𝑡),𝐹(𝑥,𝑦,𝑡):Airy's stress function
𝐼,𝐼1,𝐼2,𝐼:Inertia of the laminate
𝐽𝑘𝐼(𝑥,𝑦),𝐽𝑘𝐼𝑥((𝛼𝑚𝑥),(𝛼𝑖𝑥)),𝐽𝑘𝐼𝑦((𝛽𝑛𝑦)(𝛽𝑗𝑦)):integrant functions.
𝑀:Applied moment vector
𝑁:In-plane force vector
𝑁(𝑤):In-plane force as a function of deflection
𝑄𝑚𝑛:Laminate pressure load
Fime step in numerical integration
𝑘Forcing frequency index
𝑝:Stiffness coefficient
𝑞𝑚𝑛:Duffing equation forcing load
𝑟:Stiffness coefficient
𝑟𝑚𝑛:Stiffness in Duffing’s equation
̃𝑟:Normalized stiffness in Duffing’s equation
𝑢0(𝑥,𝑦):Mid-plane deformation in x direction
𝑢(𝑥,𝑦):Deformation in x direction
𝑣0(𝑥,𝑦):Mid-plane deformation in y direction
𝑣(𝑥,𝑦):Deformation in y direction
𝑣:Velocity ̇𝑤
̇𝑤0:Initial velocity of |𝑊(𝑡)|
|𝑊𝑚𝑛(𝑡)|:Time-dependent deflection function in modal form
𝑤(𝑥,𝑦,𝑡),𝑤(𝑥,𝑦):Approximate deflection
̇𝑤:Velocity of the approximate deflection
̈𝑤:Acceleration of the approximate deflection
𝑥0:initial deflection of |𝑊(𝑡)|
𝜀𝑖,𝜀𝑖𝑗,𝛾𝑖𝑗:Strain components
𝛼𝑚,𝛼𝑖,𝛼𝑝:Modal coefficient in 𝑥
𝛽𝑛,𝛽𝑗,𝛽𝑞:Modal coefficient in 𝑦
𝜂:Damping ratio
̃𝜂:The normalized damping ratio
𝜏:Normalized time.
𝜏𝑐:Threshold time
𝜌0(𝑘):Mass density of the 𝑘th lamina
𝜈(𝑘):Poisson ratio of the 𝑘th lamina
𝜔𝑚𝑛,𝜔:Natural frequency
Ω:Forcing frequency
𝜙(𝑥,𝑦)weighting function


The author would like to express her appreciation to Professor Stallybrass, Professor Dieci, and Professor Weiss in the School of Mathematics at Georgia Institute of Technology, for helpful discussions.