Modal Decoupling Using the Method of Weighted Residuals for the Nonlinear Elastic Dynamics of a Clamped Laminated Composite

Abstract

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 [2–4]. 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 [8–10]. 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 [2–4] 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 [2–4]. 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 [2–4]. 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 {𝜀}=𝜀⁰+𝑧𝜀¹,(2.1)

𝜀⁰=⎧⎪⎨⎪⎩𝜀𝜀^{0𝑥𝑥0𝑦𝑦}𝛾^0𝑥𝑦⎫⎪⎬⎪⎭,𝜀₁=⎧⎪⎨⎪⎩𝜅𝜅^𝑥_𝑦𝜅_𝑥𝑦𝜀⎫⎪⎬⎪⎭,(2.2a)_0𝑥𝑥=𝑢₀,_𝑥+12𝑤,^𝑥2,𝜀^0𝑦𝑦=𝑣₀,_𝑦+12𝑤,^𝑦2𝜅_𝑥,𝛾^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𝑥𝑦,_𝑥𝑦=−12𝐿(𝑤,𝑤).(2.3c) The constitutive relation between the forces and the strains is:

𝑀𝜀⁰=𝐴^∗𝐵^∗−𝐵^∗^𝑇𝐷^∗𝑁𝜀¹,(2.4) where: 𝐴^∗=𝐴⁻¹,𝐵^∗=−𝐴⁻¹𝐵,𝐷^∗=𝐷−𝐵𝐴⁻¹𝐵.(2.5a) The stiffness matrices of the laminate are defined as:

𝐶𝐴_𝑖𝑗,𝐵_𝑖𝑗,𝐷_𝑖𝑗=^ℎ2_−ℎ1𝐶^{(𝑘)𝑖𝑗}1,𝑧,𝑧²𝑑𝑧,(2.5b)_(𝑘)11=𝐶^(𝑘)22=∑^𝑘=𝑛_𝑘=1𝐸^(𝑘)1−𝜈^(𝑘)²,𝐶^(𝑘)12=𝐶^(𝑘)21=∑^𝑘=𝑛_𝑘=1𝜈^(𝑘)𝐸^(𝑘)1−𝜈^(𝑘)²,𝐶^(𝑘)66=^𝑘=𝑛_𝑘=1𝐸^(𝑘)21+𝜈^(𝑘),(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: ∇⁴𝐹=−𝑟^∗2𝐿(𝑤,𝑤).(2.9) The governing equations of motion for a laminate in a nonlinear elastic deformation are given in [21]: 𝑁𝑁_𝑥,_𝑥+𝑁_𝑥𝑦,_𝑦=0,_𝑥𝑦𝑀,_𝑥+𝑁_𝑦,_𝑦=0,_𝑥,_𝑥𝑥+2𝑀_𝑥𝑦,_𝑥𝑦+𝑀_𝑦,_𝑦𝑦+𝑁(𝑤)+𝑞=𝐼𝑤_,𝑡𝑡−𝐼₂∇²𝑤_,𝑡𝑡,(2.10) where 𝑁(𝑤)=𝑁_𝑥𝑤,_𝑥𝑥+2𝑁_𝑥𝑦𝑤,_𝑥𝑦+𝑁_𝑦𝑤,_𝑦𝑦.(2.11) Define the inertia of the laminate as 𝐼,𝐼₁,𝐼₂=^𝑁_𝑘=1^ℎ2_−ℎ1𝜌^0(𝑘)1,𝑧,𝑧²𝑑𝑧.(2.12) A homogeneous material with a symmetric layer distribution yields 𝐼₁=0. Using (2.4), (2.6) and (2.10), we obtain the equilibrium equation for the deflection of a symmetric isotropic laminate: 𝐷₁₁∇⁴𝑤+𝐼𝑤_,𝑡𝑡−𝐼₂∇²𝑤_,𝑡𝑡−𝐿(𝐹,𝑤)=𝑞,(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𝑥∈−𝑎𝑊_𝑚𝑛(𝑡)cos²𝛼_𝑚𝑥cos²𝛽_𝑛𝑦,2,𝑎2,𝑦∈−2𝑏,2𝑏,𝛼^𝑚=𝑚𝜋𝑎,𝛽^𝑛=𝑛𝜋𝑏,𝑚𝑛==22𝑘𝑗−−11,,𝑗𝑘==11,,22,,……𝑁𝑀..(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:

𝑢₀−𝑎2,𝑦,𝑡=0,𝑢⁰𝑎2,𝑦,𝑡=0,𝑣⁰𝑥,−2𝑏,𝑡=0,𝑣⁰𝑥,2𝑏,𝑡=0.(2.16b)

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 [18–20], 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

𝐸=𝐷₁₁∇⁴𝑤+𝐼𝑤_,𝑡𝑡−𝐼₂∇²𝑤_,𝑡𝑡−𝐿𝐹,𝑤−𝑞.(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_𝑚=1cos2𝛼_𝑚𝑥cos2𝛽_𝑛𝑦.(2.19) This makes (2.18) equivalent to

^{𝑏/2−𝑏/2}^{𝑎/2−𝑎/2}𝐷₁₁∇⁴𝑤+𝐼𝑤_,𝑡𝑡−𝐼₂∇²𝑤_,𝑡𝑡−𝐿𝐹,𝑤−𝑞cos2𝛼_𝑝𝑥cos2𝛽_𝑞𝑦𝑑𝑥𝑑𝑦=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:

∇⁴𝐹_𝑚𝑛=−𝑟^∗2𝛼^2𝑚𝛽^2𝑛𝑊^2𝑚𝑛cos2𝛼^𝑚co𝑥s+2𝛼c_𝑚os2𝛽_𝑛𝑦+𝑐𝑜𝑠4𝛼_𝑚𝑥+cos4𝛽_𝑛𝑦+cos4𝛼_𝑚𝑥𝑥+cos4𝛽_𝑛𝑦+2cos2𝛼_𝑚𝑥cos2𝛽_𝑛𝑦.(2.21) Solving (2.21), we obtain

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

𝐹_𝑚𝑛=𝑊^2𝑚𝑛(𝑡)𝑓₅𝑓¹cos2𝛼_𝑚𝑥+𝑓₂cos2𝛽_𝑛𝑦+𝑓₃cos4𝛼_𝑚𝑥+𝑓₄cos4𝛽_𝑛𝑦+cos2𝛼_𝑚𝑥cos4𝛽_𝑛𝑦+𝑓₆cos4𝛼_𝑚𝑥cos2𝛽_𝑛𝑦+𝑓₇cos2𝛼_𝑚𝑥cos2𝛽_𝑛𝑦(2.23a) and the coefficients 𝑓_𝑖 are given by

𝑓₁=−𝑟^∗𝛽^2𝑛32𝛼^2𝑚,𝑓₂=−𝑟^∗𝛼^2𝑚32𝛽^2𝑛,𝑓₃=−𝑟^∗𝛽^2𝑛512𝛼^2𝑚,𝑓₄=−𝑟^∗𝛼^2𝑚512𝛽^2𝑛,𝑓₅=−𝑟^∗𝛼^2𝑚𝛽^2𝑛32𝛼^2𝑚+4𝛽^2𝑛²𝑓₆,=−𝑟^∗𝛼^2𝑚𝛽^2𝑛324𝛼^2𝑚+𝛽^2𝑛²,𝑓₇=−𝑟^∗𝛼^2𝑚𝛽^2𝑛16𝛼^2𝑚+𝛽^2𝑛².(2.23b) In order to satisfy the in-plane motion constraints (2.16b) for a rigidly clamped laminate, 𝐶^1𝑚𝑛=3𝑊^𝑚𝑛(𝑡)²32𝐴¹¹𝛼^2𝑚+𝐴₁₂𝛽^2𝑛,𝐶^2𝑚𝑛=3𝑊^𝑚𝑛(𝑡)²32𝐴¹²𝛼^2𝑚+𝐴₁₁𝛽^2𝑛,(2.23c) so that for the rigidly clamped laminate

𝐹_𝑚𝑛=3𝑊^𝑚𝑛(𝑡)²64𝐴¹²𝛼^2𝑚+𝐴₁₁𝛽^2𝑛𝑥²+𝐴₁₁𝛼^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:

𝐽𝐽₀=cos2𝛼_𝑝𝑥cos²𝛽_𝑞𝑦,cos2𝛽_𝑞𝑦cos²𝛼_𝑝𝑥₁,𝐽=cos2𝛼_𝑚𝑥cos2𝛽_𝑞𝑦cos²𝛼_𝑝𝑥,cos2𝛽_𝑛𝑦cos2𝛼_𝑝𝑥cos²𝛽_𝑞𝑦₂,𝐽=cos4𝛼_𝑚𝑥cos2𝛽_𝑞𝑦cos²𝛼_𝑝𝑥,cos4𝛽_𝑛𝑦cos2𝛼_𝑝𝑥cos²𝛽_𝑞𝑦₃,𝐽=cos2𝛼_𝑚𝑥cos4𝛽_𝑛𝑦cos2𝛽_𝑞𝑦cos²𝛼_𝑝𝑥,cos4𝛼_𝑚𝑥cos2𝛽_𝑛𝑦cos2𝛼_𝑝𝑥cos²𝛽_𝑞𝑦₄,𝐽=cos4𝛼_𝑚𝑥cos2𝛽_𝑛𝑦cos2𝛽_𝑞𝑦cos²𝛼_𝑝𝑥,cos2𝛼_𝑚𝑥cos4𝛽_𝑛𝑦cos2𝛼_𝑝𝑥cos²𝛽_𝑞𝑦₅,𝐽=cos2𝛼_𝑚𝑥cos2𝛽_𝑛𝑦cos2𝛽_𝑞𝑦cos²𝛼_𝑝𝑥,cos2𝛼_𝑚𝑥cos2𝛽_𝑛𝑦cos2𝛼_𝑝𝑥cos²𝛽_𝑞𝑦₆,𝐽=sin2𝛼_𝑚𝑥sin4𝛽_𝑛𝑦sin2𝛽_𝑞𝑦sin2𝛼_𝑝𝑥,sin4𝛼_𝑚𝑥sin2𝛽_𝑛𝑦sin2𝛼_𝑝𝑥sin2𝛽_𝑞𝑦₇,=sin2𝛼_𝑚𝑥sin2𝛽_𝑛𝑦sin2𝛼_𝑝𝑥sin2𝛽_𝑞𝑦.(2.26) Each of the above pairs contains a product of a function of x and a function of 𝑦, that is,

𝐽^𝑘𝐼(𝑥,𝑦)=𝐽^𝑘𝐼𝑥𝑓₁𝛼_𝑚𝑥,𝑓₂𝛼_𝑝𝑥𝐽^𝑘𝐼𝑦𝑔₁𝛽_𝑛𝑦,𝑔₂𝛽_𝑞𝑦𝑘=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𝐼𝑥𝑓=₁𝐽𝛼_{𝑚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𝐷¹¹𝛼^2𝑚+𝛽^2𝑛²𝑟𝐼,(2.31a)_𝑚𝑛=116𝐼6𝐴¹¹+𝐴₂₁₁𝛼^4𝑚+2𝐴₁₂𝛼^2𝑚𝛽^2𝑛+𝐴₁₁𝛽^4𝑛−𝐴²¹²𝐴₁₁9𝛼^4𝑚+𝛽^4𝑛+4𝛼^4𝑚𝛽^4𝑛𝛼^2𝑚+4𝛽^2𝑛²+4𝛼^4𝑚𝛽^4𝑛𝑞4𝛼^2𝑚+𝛽^2𝑛²,(2.31b)_𝑚𝑛=𝑄^𝑚𝑛𝐼=𝐼𝐼+,𝛼_2𝑚(2.31c)+𝛽^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

𝑟^𝐿𝑚𝑛=116𝐼𝐴²¹¹−𝐴²¹²𝐴₁₁9𝛼^4𝑚+𝛽^4𝑛+4𝛼^4𝑚𝛽^4𝑛𝛼^2𝑚+4𝛽^2𝑛²+4𝛼^4𝑚𝛽^4𝑛4𝛼^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 𝐴₁₁>𝐴₁₂ 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

𝑑²𝑊𝑑𝜏²+𝑊+𝜂̃𝑊̇+𝑟𝑊̃³=𝑓cos(𝑘𝜏),(3.1) where 𝜏=𝜔𝑡, 𝑟̃=𝑟/𝜔²,𝜂̃=𝜂/𝜔²,𝑓=𝑓/𝜔²,𝑘=Ω/𝜔, 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 [26–28], 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 𝜔²<0, or a system with a negative stiffness 𝑟<0. Consequently, the results may not directly apply to the Duffing equation with 𝜔²>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 𝜔²<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

𝐿(𝑊,𝑉)=𝑉²2+𝑊2²+𝑟𝑊2⁴>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)=𝑊+𝑟𝑊³𝑊̇+𝑉𝑉̇=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

𝑊̈+𝑊+𝑟𝑊̃³=𝑓1−(𝑘𝜏)²2!+(𝑘4𝜏!)⁴−(𝑘6𝜏!)⁶+(𝑘8𝜏!)⁸+⋯+(−1)^𝑛(𝑘𝜏)^2𝑛=𝑓𝑔(𝑘𝜏),𝑛=0,1…(2𝑛)!(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

𝑊̈+𝜔²𝑊+𝑟𝑊̃³=0,(3.5a) where

𝜔²=1−𝜇1−(𝑘𝜏)²2+(𝑘4𝜏!)⁴−(𝑘6𝜏!)⁶+(𝑘8𝜏!)⁸+⋯=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 𝜔²≈1 for all values of 𝑘𝜏. When 𝜇∼1 or 𝜇≫1, 𝜔² in (3.5b) can become negative. The transition time 𝜏_𝑐, when the system passes out of the stable system at 𝜔²=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 𝜔²<0, is:

𝜏_c=max⎧⎪⎨⎪⎩12k,1k6+23+6𝜇⎫⎪⎬⎪⎭.(3.6) This threshold time 𝜏_𝑐 signifies the system’s transition into an unstable system with 𝜔²<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 130∗180∗1.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 𝑄₁₁=12cos(𝑘𝜔₁₁𝑡)N/cm². 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 𝑞₁₁/𝑄₁₁=0.5215 for the fundamental mode, Mode-1, whereas the normalized force 𝑞 and the pressure load 𝑄 are measured in [mm] and [N/cm²], respectively.

Table 1 

Physical parameters of a PWB.

Table 2 

Modal parameters of the Duffing equation.

Figure 1 

A symmetric PWB (130∗180∗1.54 mm).

Figures 2(a) and 2(b) represent the same Poincare map of Mode-1 in a quasiperiodic behavior with finite periodic orbits subject to, 𝑄₁₁=2.4cos(𝑘𝜔₁₁𝑡)N/cm², 𝑘=1, with initial condition 𝑥₀=0, and𝑣₀=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 𝑥₀=2mm and a load 𝑄₁₁=24cos(𝑘𝜔₁₁𝑡)N/cm²,𝑘=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 𝑥₀=20mm while the load remains at 𝑄₁₁=24cos(𝑘𝜔₁₁𝑡)N/cm²,𝑘=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.

(a) Map with connected periodic points

(b) Map with disconnected periodic points

Figure 2 

Mode-1 Poincare map in finite periodic orbits subject to the harmonic load 𝑄11=2.4cos(𝑘𝑤11𝑡)N/cm2,𝑘=1.

(a) 𝑥0=2mm,𝑄11=24N/cm2,𝑘=1

(b) 𝑥0=20mm,𝑄11=24N/cm2,𝑘=1.

(c) 𝑥0=0,𝑄11=432N/cm2,𝑘=1.

(d) 𝑥0=0,𝑄11=4800N/cm2,𝑘=1.

Figure 3 

Mode-1 response with respect to the harmonic load 𝑄=𝑄11cos(𝑘𝜔11𝑡)N/cm2

Figure 3(c) shows the Mode-1 Poincare map in chaos subject to a large load of 𝑄₁₁=432cos(𝑘𝜔₁₁𝑡)N/cm²,𝑘=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 𝜔²>0 and 𝜔²<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 𝑄₁₁=4800cos(𝑘𝜔₁₁𝑡)N/cm²,𝑘=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𝑏/2∫^0𝑎/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:

𝑈=𝐷2^0𝑏/2^0𝑎/2∇²𝑤+12ℎ²𝑒₁₂−2(1−𝜈)12ℎ²𝑒₂+𝜕²𝑤𝜕𝑥²𝜕²𝑤𝜕𝑦²−𝜕²𝑤𝜕𝑥𝜕𝑦²𝑑𝑥𝑑𝑦,(5.1) where the first and the second strain invariants are defined as

𝑒₁=𝜀_𝑥+𝜀_𝑦,𝑒₂=𝜀_𝑥𝜀_𝑦−14𝛾^2𝑥𝑦.(5.2) Berger’s method applied the principle of virtual work with the assumption 𝑒₂=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. [18–20] 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:

^𝑎−𝑎cosh2𝛼_𝑚𝑥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.

Nomenclature

𝐴:	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
𝑤(𝑥,𝑦,𝑡):	deflection
𝑤(𝑥,𝑦,𝑡),𝑤(𝑥,𝑦):	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
𝜇:	Factor
𝜏:	Normalized time.
𝜏𝑐:	Threshold time
𝜌0(𝑘):	Mass density of the 𝑘th lamina
𝜈(𝑘):	Poisson ratio of the 𝑘th lamina
𝜔𝑚𝑛,𝜔:	Natural frequency
Ω:	Forcing frequency
𝜙(𝑥,𝑦)∶	weighting function

Acknowledgments

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.