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

The strain field satisfies the compatibility equation:

Define an operator

equation (2.3a) can be expressed as

The constitutive relation between the forces and the strains is:

where: The stiffness matrices of the laminate are defined as:

where is the number of the layers in symmetry. For a symmetric isotropic laminate, , . In (2.4), and are the in-plane force and moment vectors, respectively. Introducing the Airy stress function defined by from (2.4), the compatibility (2.3c) can be expressed as

where a symmetric isotropic laminate has Therefore, for the symmetric isotropic laminate, (2.7) can be further reduced to: The governing equations of motion for a laminate in a nonlinear elastic deformation are given in [21]: where Define the inertia of the laminate as A homogeneous material with a symmetric layer distribution yields . Using (2.4), (2.6) and (2.10), we obtain the equilibrium equation for the deflection of a symmetric isotropic laminate: where: 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

This expression satisfies the actual boundary conditions, namely:

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 [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

Using MWR, we require that

where the weighting function has the same index as that of the approximating function in (2.15), given as:

This makes (2.18) equivalent to

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:

Solving (2.21), we obtain

where

and the coefficients are given by

In order to satisfy the in-plane motion constraints (2.16b) for a rigidly clamped laminate, so that for the rigidly clamped laminate

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 , so that

To evaluate the integral (2.20), we introduce the following symmetric pairs for integration.

(19) These pairs are defined by:

Each of the above pairs contains a product of a function of and a function of , that is,

In (2.27), not all functions and need to be present in each pair. Noticeably, the two functions and in each pair are symmetric with respect to their dependent variables, and , that is,

This means by switching the argument with , respectively, in the corresponding function, the two components in each pair is the same, equivalently;

Therefore, the evaluation of each integral pair denoted by , respectively, in the spatial coordinates has the following properties: 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

where the coefficients for an isotropic laminate in a rigidly clamped boundary condition are

For a loosely clamped laminate, we obtain an identical Duffing equation with the same and expressions as above, except that the stiffness coefficient is given by

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 which yields in (2.31b) and 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

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 , or a system with a negative stiffness . Consequently, the results may not directly apply to the Duffing equation with and a positive stiffness 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 , 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

The function is positive over the plane, except at the center . The stability of the system can be measured by the time derivative of the Lyapunov function, which vanishes at the center:

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 in the form of a series, so that the Duffing equation (3.1), without damping becomes

The series in expression (3.4) satisfies . Assuming the load is in proportion to the deflection, that is, we obtain after rearranging terms

where

Here the factor , and represent the large load, the small load and the intermediate load regimes, respectively. Equation (3.5b) implies that when , the response of the system is almost periodic since for all values of . When or , in (3.5b) can become negative. The transition time , when the system passes out of the stable system at 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 with , based on retaining only three terms in (3.5b) for , is:

This threshold time signifies the system’s transition into an unstable system with . 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 mm 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 , , and of that of the fundamental mode response for the seconh, third and fourth modes under 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 for the fundamental mode, Mode-1, whereas the normalized force and the pressure load are measured in and , 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, , , with initial condition , and. 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 and a load . 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 while the load remains at . 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**

(a) |

(b) |

(c) |

(d) |

Figure 3(c) shows the Mode-1 Poincare map in chaos subject to a large load of , 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 and , 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 , which is in contrast to the homoclinic orbit around the center of a stable Duffing system with [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 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 . 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:

where the first and the second strain invariants are defined as

Berger’s method applied the principle of virtual work with the assumption , 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, 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:

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 |

: | Stiffness coefficients of the laminate |

: | 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 |

: | 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 |

: | Mid-plane deformation in x direction |

: | Deformation in x direction |

: | Mid-plane deformation in y direction |

: | Deformation in y direction |

: | Velocity |

: | Initial velocity of |

: | Time-dependent deflection function in modal form |

: | deflection |

,: | Approximate deflection |

: | Velocity of the approximate deflection |

: | Acceleration of the approximate deflection |

: | initial deflection of |

: | Strain components |

: | Modal coefficient in |

: | Modal coefficient in |

: | Damping ratio |

: | The normalized damping ratio |

: | Factor |

: | Normalized time. |

: | Threshold time |

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