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

Authors: Xiaoling He
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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.


---

## Body

## 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}+z{ε1},(2.2a){ε0}={εxx0εyy0γxy0},{ε1}={κxκyκxy},(2.2)εxx0=u0,x+12w,x2,εyy0=v0,y+12w,y2,γxy0=u0,y+v0,x+w,xw,y,κx=-w,xx,κy=-w,yy,κxy=-2w,xy.
The strain field satisfies the compatibility equation:(2.3a)εxx0,yy+εyy0,xx-γxy0,xy=w,xy2-w,xxw,yy.
Define an operator(2.3b)L(p,q)=p,yyq,xx+p,xxq,yy-2p,xyq,xy,
equation (2.3a) can be expressed as(2.3c)εxx0,yy+εyy0,xx-γxy0,xy=-12L(w,w).
The constitutive relation between the forces and the strains is:(2.4){ε0M}=[A*B*-(B*)TD*]{Nε1},
where:
(2.5a)A*=A-1,B*=-A-1B,D*=D-BA-1B.
The stiffness matrices of the laminate are defined as:(2.5b)[Aij,Bij,Dij]=∫-h1h2Cij(k)[1,z,z2]dz,(2.5c)C11(k)=C22(k)=∑k=1k=nE(k)1-(ν(k))2,C12(k)=C21(k)=∑k=1k=nν(k)E(k)1-(ν(k))2,C66(k)=∑k=1k=nE(k)2(1+ν(k)),
where n is the number of the layers in symmetry. For a symmetric isotropic laminate, [B]=0, D*=D. In (2.4), N and M are the in-plane force and moment vectors, respectively. Introducing the Airy stress function defined by
(2.6)Nx=F,yy,Ny=F,xx,Nxy=-F,xy,
from (2.4), the compatibility (2.3c) can be expressed as(2.7)F,xxxx+p*F,xxyy+F,yyyy=-r*2L(w,w),
where a symmetric isotropic laminate has
(2.8)r*=1A22*,p*=2A12*+A66*A22*=2.
Therefore, for the symmetric isotropic laminate, (2.7) can be further reduced to:
(2.9)∇4F=-r*2L(w,w).
The governing equations of motion for a laminate in a nonlinear elastic deformation are given in [21]:
(2.10)Nx,x+Nxy,y=0,Nxy,x+Ny,y=0,Mx,xx+2Mxy,xy+My,yy+N(w)+q=Iw,tt-I2∇2w,tt,
where
(2.11)N(w)=Nxw,xx+2Nxyw,xy+Nyw,yy.
Define the inertia of the laminate as
(2.12)[I,I1,I2]=∑k=1N∫-h1h2ρ0(k)[1,z,z2]dz.
A homogeneous material with a symmetric layer distribution yields I1=0. Using (2.4), (2.6) and (2.10), we obtain the equilibrium equation for the deflection of a symmetric isotropic laminate:
(2.13)D11∇4w+Iw,tt-I2∇2w,tt-L(F,w)=q,
where:
(2.14)L(F,w)=∑i=1I∑j=1J∑n=1N∑m=1M(Fmn,xxwij,yy+Fmn,yywij,xx-2Fmn,xywij,xy).
Here Fmn are the modal functions of the Airy stress function F(x,y,t), and wij are the modal components of the deflection function w(x,y,t). For either the rigidly clamped or the loosely clamped laminates, ŵ(x,y,t) is chosen as the deflection approximating function in series as(2.15)ŵ(x,y,t)=∑n=1N∑m=1MWmn(t)*w(x,y)=∑n=1N∑m=1MWmn(t)cos2(αmx)cos2(βny),x∈[-a2,a2],y∈[-b2,b2],αm=mπa,βn=nπb,m=2k-1,k=1,2,…M.n=2j-1,j=1,2,…N.
This expression satisfies the actual boundary conditions, namely:(2.16a)ŵ(x,-b2,t)=0,ŵ(x,b2,t)=0,ŵ(a2,y,t)=0,ŵ(a2,y,t)=0,ŵ,y(x,-b2,t)=0,ŵ,y(x,b2,t)=0,ŵ,x(-a2,y,t)=0,ŵ,x(a2,y,t)=0.
In addition, the rigidly clamped laminate must satisfy the in-plane motion constraints of the mid-plane at the boundaries as:(2.16b)u0(-a2,y,t)=0,u0(a2,y,t)=0,v0(x,-b2,t)=0,v0(x,b2,t)=0.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 ŵ(x,y,t), as(2.17)Ẽ=D11∇4ŵ+Iŵ,tt-I2∇2ŵ,tt-L(F,ŵ)-q.
Using MWR, we require that(2.18)∫-b/2b/2∫-a/2a/2[Ẽϕ(x,y)]dxdy=0,
where the weighting function has the same index as that of the approximating function in (2.15), given as:(2.19)ϕ(x,y)≡ŵ,xxyy(x,y)=∑n=1N∑m=1Mcos(2αmx)cos(2βny).
This makes (2.18) equivalent to(2.20)∫-b/2b/2∫-a/2a/2{[D11∇4ŵ+Iŵ,tt-I2∇2ŵ,tt-L(F,ŵ)-q]cos(2αpx)cos(2βqy)}dxdy=0.
As a preliminary to evaluate integral (2.20), various calculations are necessary. Using the deflection approximation ŵ(x,y,t) in (2.15), the function in (2.9) becomes:(2.21)∇4Fmn=-r*2αm2βn2Wmn2[cos(2αmx)+cos(2βny)+cos(4αmx)+cos(4βny)+cos(4αmx)cos(2αmx)+cos(4βny)+2cos(2αmx)cos(2βny)].
Solving (2.21), we obtain(2.22)Fmn=12Cmn1y2+12Cmn2x2+F̃mn,
where(2.23a)F̃mn=Wmn2(t)[f1cos(2αmx)+f2cos(2βny)+f3cos(4αmx)+f4cos(4βny)+f5cos(2αmx)cos(4βny)+f6cos(4αmx)cos(2βny)+f7cos(2αmx)cos(2βny)]
and  the coefficients fi are given by(2.23b)f1=-r*βn232αm2,f2=-r*αm232βn2,f3=-r*βn2512αm2,f4=-r*αm2512βn2,f5=-r*αm2βn232(αm2+4βn2)2,f6=-r*αm2βn232(4αm2+βn2)2,f7=-r*αm2βn216(αm2+βn2)2.
In order to satisfy the in-plane motion constraints (2.16b) for a rigidly clamped laminate,
(2.23c)Cmn1=3[Wmn(t)]232(A11αm2+A12βn2),Cmn2=3[Wmn(t)]232(A12αm2+A11βn2),
so that for the rigidly clamped laminate(2.24)Fmn=3[Wmn(t)]264((A12αm2+A11βn2)x2+(A11αm2+A12βn2)y2)+F̃mn,
where F̃mn 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 Cmn1=Cmn2=0, so that(2.25)F̃mn=Fmn.
To evaluate the integral (2.20), we introduce the following symmetric pairs for integration.JI={JI1(αmx,βny)JI1(αpx,βqy),JI2(αmx,βny)JI2(αpx,βqy)}I=0,1,...7. (19) These pairs are defined by:(2.26)J0={cos(2αpx)cos2(βqy),cos(2βqy)cos2(αpx)},J1={cos(2αmx)cos(2βqy)cos2(αpx),cos(2βny)cos(2αpx)cos2(βqy)},J2={cos(4αmx)cos(2βqy)cos2(αpx),cos(4βny)cos(2αpx)cos2(βqy)},J3={cos(2αmx)cos(4βny)cos(2βqy)cos2(αpx),cos(4αmx)cos(2βny)cos(2αpx)cos2(βqy)},J4={cos(4αmx)cos(2βny)cos(2βqy)cos2(αpx),cos(2αmx)cos(4βny)cos(2αpx)cos2(βqy)},J5={cos(2αmx)cos(2βny)cos(2βqy)cos2(αpx),cos(2αmx)cos(2βny)cos(2αpx)cos2(βqy)},J6={sin(2αmx)sin(4βny)sin(2βqy)sin(2αpx),sin(4αmx)sin(2βny)sin(2αpx)sin(2βqy)},J7={sin(2αmx)sin(2βny)sin(2αpx)sin(2βqy)}.
Each of the above pairs contains a product of a function of x and a function of y, that is,(2.27)JIk(x,y)=JIxk[f1(αmx),f2(αpx)]JIyk[g1(βny),g2(βqy)]k=1,2;I=0,1,…7.
In (2.27), not all functions  fi and  gj need to be present in each pair. Noticeably, the two functions  JI1 and  JI2  in each pair are symmetric with respect to their dependent variables, xand y, that is,(2.28a)JIx1[f1(αmx),f2(αpx)]JIy1[g1(βny),g2(βqy)]=JIx2[g1(αmx),g2(αpx)]JIy2[f1(βny),f2(βqy)]I=0,1,…7.
This means by switching the argument (αmx),(αpx) with (βny),(βqy), respectively, in the corresponding function, the two components in each pair is the same, equivalently;(2.28b)JI1(f[(αmx),(αpx)],g[(βny),(βqy)])=JI2(f[(βny),(βqy)],g[(αmx),(αpx)])I=0,1,…7.
Therefore, the evaluation of each integral pair denoted by k=1,2, respectively, in the spatial coordinates has the following properties:
(2.29)∫-b/2b/2∫-a/2a/2JIk(x,y)dxdy=∫-a/2a/2JIxk(x)dx*∫-b/2b/2JIyk(y)dy,k=1,2,I=0,1,…7,∫-b/2b/2∫-a/2a/2JI1dxdy=∫-b/2b/2∫-a/2a/2JI2dxdy,∫-b/2b/2∫-a/2a/2JI1w,xxyydxdy=∫-b/2b/2∫-a/2a/2JI2w,xxyydxdy.
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.30)Ẅmn+ωmn2Wmn+rmnWmn3=qmn,
where the coefficients for an isotropic laminate in a rigidly clamped boundary condition are(2.31a)ωmn2=16D11(αm2+βn2)2Ĩ,(2.31b)rmn=116Ĩ{6[A11αm4+2A12αm2βn2+A11βn4]+A112-A122A11[9(αm4+βn4)+4αm4βn4(αm2+4βn2)2+4αm4βn4(4αm2+βn2)2]},(2.31c)qmn=QmnĨ,(2.31d)Ĩ=I+(αm2+βn2)I2.
For a loosely clamped laminate, we obtain an identical Duffing equation with the same ωmn2 and qmn expressions as above, except that the stiffness coefficient is given by(2.31e)rmnL=116Ĩ{A112-A122A11[9(αm4+βn4)+4αm4βn4(αm2+4βn2)2+4αm4βn4(4αm2+βn2)2]}.
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 A11>A12 and  Aij>0,which yields rmn>0 in  (2.31b) and rmnL>0 in (2.31e), respectively. The loosely clamped boundary condition has a reduced stiffness, that is, rmnL<rmn. Equations (2.31a), (2.31b), and (2.31e) indicate that both the natural frequency ωmn and the stiffness rmn and rmnL decrease with an increase of the inertia Ĩ, due to the increased density and thickness. On the contrary, a higher rigidity Dij and stiffness Aij 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  αmorβn. In comparison, a higher mode has a higher ωmn and rmn, or rmnL, since amn and βmn are amplified by the modal integers m and n.


## 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(3.1)d2Wdτ2+W+η̃Ẇ+r̃W3=f̃cos(kτ),
where τ=ωt, r̃=r/ω2,η̃=η/ω2,f̃=f/ω2,k=Ω/ω, 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 ω2<0, or a system with a negative stiffness r<0. Consequently, the results may not directly apply to the Duffing equation with ω2>0 and a positive stiffness r>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 V=Ẇ=dW/dt, the weak Lyapunov function of the system is(3.2)L(W,V)=V22+W22+rW42>0.
The function L(W,V) is positive over the plane, except at the center L(W,V)=L(0,0)=0. The stability of the system can be measured by the time derivative of the Lyapunov function, which vanishes at the center:(3.3)L̇(0,0)=(W+rW3)Ẇ+VV̇=0.
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 f̃cos(kτ) in the form of a series, so that the Duffing equation (3.1), without damping becomes(3.4)Ẅ+W+r̃W3=f̃(1-(kτ)22!+(kτ)44!-(kτ)66!+(kτ)88!+⋯+(-1)n(kτ)2n(2n)!)=f̃g(kτ),n=0,1…
The series g(kτ) in expression (3.4) satisfies |g(kτ)|≤1. Assuming the load is in proportion to the deflection, that is, f̃=μW, we obtain after rearranging terms(3.5a)Ẅ+ω̃2W+r̃W3=0,
where(3.5b)ω̃2=1-μ(1-(kτ)22+(kτ)44!-(kτ)66!+(kτ)88!+⋯)=1-μg(kτ).
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 ω̃2≈1 for all values of kτ. When μ~1 or μ≫1, ω̃2 in (3.5b) can become negative. The transition time τc, 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 τc when kτ<1 with μ~1, based on retaining only three terms in (3.5b) for ω̃2<0, is:(3.6)τc=max{12k,1k6+23+6μ}.
This threshold time τc 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 of130*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 Q11=12cos(kω11t)N/cm2. The computation specifies the pressure load to the laminate Qmn, instead of the normalized load qmn in (2.31c). The conversion between the pressure load and the normalized load is q11/Q11=0.5215 for the fundamental mode, Mode-1, whereas the normalized force q and the pressure load Q are measured in [mm] and [N/cm2], respectively.Table 1
Physical parameters of a PWB.
LayerMaterialSpecific weight  (g/cm3)Young's module(MPa)Poisson ratioThickness    (mm)1Cu/FR-47.59917.9*1040.2050.0302FR-41.20018.6*1040.1900.2003Cu8.31011.8*1040.4000.0354FR-41.20018.6*1040.1901.000Table 2
Modal parameters of the Duffing equation.
ModeNatural frequencyωmn,[(s)-1]Natural frequency ratioStiffnessrmn,[(mm·s)-2]M1:m=1, n=16.7839e37.1888e6M2:m=1, n=32.5381e4ω31/ω11=2.681.2548e8M3:m=3, n=14.2430e4ω31/ω13=1.684.3085e8M4.m=3, n=36.1008e4ω33/ω31=1.455.8140e8Figure 1
A symmetric PWB (130*180*1.54 mm).Figures2(a) and 2(b) represent the same Poincare map of Mode-1 in a quasiperiodic behavior with finite periodic orbits subject to, Q11=2.4cos(kω11t)N/cm2, k=1, with initial condition x0=0,  andv0=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 x0=2mm and a load Q11=24cos(kω11t)N/cm2,k=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 x0=20mm while the load remains at Q11=24cos(kω11t)N/cm2,k=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.Mode-1 Poincare map in finite periodic orbits subject to the harmonic loadQ11=2.4cos(kw11t)N/cm2,k=1.
(a)
Map with connected periodic points(b)
Map with disconnected periodic pointsMode-1 response with respect to the harmonic loadQ=Q11cos(kω11t)N/cm2
(a)
x0=2mm,Q11=24N/cm2,k=1(b)
x0=20mm,Q11=24N/cm2,k=1.(c)
x0=0,Q11=432N/cm2,k=1.(d)
x0=0,Q11=4800N/cm2,k=1.Figure3(c) shows the Mode-1 Poincare map in chaos subject to a large load of Q11=432cos(kω11t)N/cm2,k=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 ωmn2<0, which is in contrast to the homoclinic orbit around the center of a stable Duffing system with ωmn2>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 Q11=4800cos(kω11t)N/cm2,k=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 Figures3(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 ∫0b/2∫0a/2(L(F,w)w)dxdy=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:(5.1)U=D2∫0b/2∫0a/2{[(∇2w+12h2e12)]-2(1-ν)[12h2e2+∂2w∂x2∂2w∂y2-(∂2w∂x∂y)2]}dxdy,
where the first and the second strain invariants are defined as(5.2)e1=εx+εy,e2=εxεy-14γxy2.
Berger’s method applied the principle of virtual work with the assumption e2=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 ŵ,xxyy(x,y) for the weighting function, where ŵ(x,y) 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,
(5.3)ŵ(x,y,t)=∑n=1N∑m=1MWmn(t)Xm(αmx)Xn(βny),x∈[-a,a],y∈[-b,b],Xm(αmx)Xn(βny)=(cosh(αmx)cosh(αma)-cos(αmx)cos(αma))(cosh(βny)cosh(βnb)-cos(βny)cos(βnb)).
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:(5.4)∫-aacosh(2αmx)cosh(αjx)cosh(αkx)dx.
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.



---

*Source: 972930-2009-07-07.xml*