Effect of large deflection on laminated composite plates

Abstract

Keywords: Finite Element Analysis; Vibration; Laminated Composite Plate; Pre-Stressed Plate; Large Deflection

Introduction

With the development of technology, the need for more lightweight, higher resistance, and more economical materials is increasing day by day. Therefore, composite materials are widely used in engineering structures due to their lightness, flexibility, corrosion resistance, and rigidity [1]. Some of the engineering structures in the assembly stage are produced in a flat or curved form. However, some of

them are curved during assembly due to their usage places and used in the pre-stressed state. The rigidity, which varies depending on the geometry and material of the structure, also affects the dynamic behavior of the plates. For this reason, knowing the dynamic characteristics of the structures is vital in terms of predicting the behavior of the structures under dynamic loads.

*Corresponding author. *E-mail address: can.gonenli@ege.edu.tr This paper was recommended for publication in revised form by Regional Editor Pravin Katare Published by Yıldız Technical University Press, İstanbul, Turkey Copyright 2021, Yıldız Technical University. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Several studies have been presented about static and dynamic analyzes of curved plates with different geometries and materials in the literature. Wang et al. analyzed composite fiber curve laminates and then gets a damping ratio of each laminated plate based on a modal experiment using the classical plate theory [2]. Tornabene et al. examined the free vibration nature of laminated composite thick and moderately thick elliptic cones, cylinders, and plates [3]. Ye et al. developed a high-performance semi-analytical numerical model to analyze the bending responses of the angle-ply composite laminated cylindrical shells with the fiber-reinforced layers using the scaled boundary finite element method [4]. Vidal et al. studied the modeling of laminated composite and sandwich shells through a variable separation approach based on Reissner’s Variational Mixed Theorem [5]. Sahoo et al. analyzed the geometrically nonlinear deflection responses of glass/epoxy composite flat/curved shell panel structure theoretically with the help of three different displacement field kinematics and Green– Lagrange strain–displacement relation [6]. Civalek studied free vibration analysis of laminated composite panels and curved plates with functionally graded materials using Love’s shell theory and first-order shear deformation theory [7]. Aurojyoti et al. presented a polygonal finite element for nonlinear analysis of laminated plates using Reddy’s thirdorder shear deformation theory [8]. Kormanikova presented a numerical approach of mode-frequency analysis of a simply-supported laminated doubly curved shell. For the laminated shell, the first-order shear deformation theory is stated that accurately predicting shell behavior [9]. Dynamic analogues of von Karman-Donnell type shell equations for doubly curved, thin isotropic shells in rectangular platform are formulated and expressed in displacement components by Nath and Sandeep [10]. Anamagh and Bediz focused on the dynamics of doubly-curved functionally graded and laminated composite structures with arbitrary geometries and boundary conditions using an energy-based approach where the strain energy of the structure is expressed using three-dimensional elasticity equations [11]. There are other several studies which employed first-order shear deformation theory (FSDT) and higher-order shear deformation theory (HSDT) to examine various structures for many purposes [12, 13]. Structural elements encounter different loading cases under their service conditions [14]. The curved plates with/ without pre-stress are used in the production of engineering structures, depending on the desired conditions. In this context, the effect of large deflection on the natural frequencies and mode shapes of the pre-stressed curved laminated composite plates is the main goal of this study. Although there are various studies in the literature about laminated composite plates, there have not been any published papers on this particular topic and the pre-stressed curved laminated composite plates have not been studied before. The transformation of the structure into a pre-stressed curved

form by bending with an external force, and transforming it into a structure with different curvatures gives originality to the work. In order to validate the presented approach which is modeled in MATLAB, the flat and curved plate models are also created with a finite element software program, ANSYS, and natural frequency values are compared.

Mathematical Model

The main difference between thick and thin plates is that with thin plates, the transverse shear strains are negligible whilst for thick plates they are not [15]. Shear effects are considered in plate theories such as first-order shear deformation theory (FSDT) and higher-order shear deformation theory (HSDT). However, the use of these theories is not required in thin plates where shear effects are not observed. Especially in thin and ultra-thin plates, very good and consistent results can be obtained with the classical plate theory (CPT) theory. The four-node quadrilateral element is selected to create flat and curved laminated composite plate models using classical plate theory for free vibration analysis. The plate must be rotatable about any of its axes to create a curvilinear model. Therefore, the selected finite element model must also contain six generalized coordinates in each node. Figure-1 shows the laminated composite plate model. The strain energy equation for the thin plate element is given in equation (1) [15]. U=

where σ is a stress matrix, and ε is a strain matrix of the plate. Stress-strain relation can be expressed in the form of equation (2) for a thin plate. σ = Dε

Figure 1. Degrees of freedom for the quadrilateral rectangular plate.

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

1. NL k {ε}T [D]{ε}dz dA

Where NL is the number of laminae, zk, zk-1 are the top and bottom coordinates of the kth lamina, and D is the material matrix. The stress-strain relation for the laminated composite plate contains three displacements (u0, v0, w0) and can be expressed in matrix form:  ∂2 w   ∂u      2  ∂x   ε x   ∂x   ∂2 w     ∂v  { ε} =  ε y  =  −z 2   ∂y   γ   ∂y   xy   ∂u ∂v   ∂2 w  2   +   ∂x ∂y   ∂y ∂x 

Substituting equation (9) into equation (7) gives the quadratic shape function and it is given in equation (10) representing the flexural vibration of the plate.  1 (1 + ξ ξ )(1 + η η)(2 + ξ ξ + η η − ξ 2 − η 2 ) j j j j 8    b  N1Tj (ξ,η)  (1 + ξj ξ )( ηj + η)(η 2 − 1)   2   a 2 − ( ξ j +ξ )(ξ − 1)(1 + ηjη)   2  

 4 ∂2   2 2   (a) ∂ ξ   4 ∂2   N1 ( , ) B1 =  2 2  (b) ∂ η    ∂2   8  (a)(b) ∂ ξ ∂ η   

The kinetic energy equation for the thin plate element is given in equation (6) [15]. T=

Two different shape function is necessary to satisfy the stress-strain and displacement relations. The first shape function satisfies the flexural vibration, while the second one satisfies the in-plane vibration of the plate. Flexural Vibration of the Plate The flexural vibration theory is based on the bending vibration of the plate. The finite element model to be used for the bending vibration type must have a total twelve degrees of freedom corresponding to three generalized coordinates, which are w, dw/dy, and dw/dx at each node. The displacement function is given in equation (7). w = a1 + a2x + a3y + a4x2 + a5xy + a6y2 + a7x3 + a8x2y + a9xy2 + a10y3 + a11x3y + a12xy3

The flexural displacement, w, is written in terms of natural coordinates (ξ, η) instead of cartesian coordinates (x, y). Hence, dw/dy and dw/dx become:

The strain energy equation derived from equation (3) for flexural vibration is given in equation (12). U1 =

where ε1 is strain relation for flexural vibration and given in equation (13): ε1 = B1q1

Substituting equation (13) into equation (12): 1 1 1 U1 = q1T  ∫ ∫ [ B1 ]T [DB ][ B1 ]d ξ d η q1 −  2  1 −1

Element stiffness matrix [K1] can be calculated and it is given in equation (15). 1

The expanded kinetic energy equation is given in equation (16).

(ξj, ηj) are the coordinates of node j. The strain displacement matrix B1 is obtained by equation (11).

The matrix form of displacement function in terms of ξ and η coordinates is:

1 1 1 T1 = q1T ρ h∫ ∫ [N1 ]T [N B ][ B1 ]d ξ d η  q1 − 1 −1  2 

Element mass matrix [M1] is calculated from equation (16) and is given in equation (17).

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

u = N21q21 + N22q23 + N23q25 + N24q27 v = N21q22 + N22q24 + N23q26 + N24q28

In-Plane Plate Theory The in-plane vibration theory contains total of eight displacements which are uj and vj (j=1 to 4). The strain energy equation derived from equation (3) for in-plane vibration is given in equation (18). NL 1 U 2 = ∫ { ε 2}T ∑[D k ](z k − z k −1 ){ ε 2} dA 2A k =1

where q2i (i=1 to 8) represents the element displacement vector. The shape function for satisfying the stress-strain relation is given in equation (26), where ξj, ηj are the coordinates of node j in the rectangular plate model. 1 N 2i = (1 + ξ i ξ )(1 + η i η) 4

where ε2 is strain relation for in-plane vibration and given in equation (19):

The in-plane element stiffness matrix [K2] is obtained with using matrix B from equation (21), and in-plane material matrix DL.

Substituting equation (19) into equation (18): 1 1 1 U 2 = q2T  ∫ ∫ [ B2 ]T [DL ][ B2 ]det J d ξ dη  q2  2  −1 −1

1 1 1 T1 = q2T  ρh∫ ∫ [N 2T N 2 det J d ξ dη  q2 − 1 −1  2 

where J is the Jacobian matrix.  ∂x ∂y   ∂ ξ ∂ξ   J=  ∂x ∂y   η η ∂ ∂ 

[K]21 and [K]12, obtained from the substitution of equation (5) and equation (3), give the couple effect between flexural and in-plane strain relation. These matrices can be calculated through equation (30). 1

From the interpolation equations, we have:  ∂u   ∂ξ     ∂u   η ∂  Gq2 =    ∂v  ∂ξ    ∂v   ∂η   

where A2 is given by:  J 22 − J12 0 1  A2 = 0 0 − J 21 det J   − J 21 J11 J 22 

The element mass matrix [M2] can also be found from the kinetic energy.

The element displacement components u and v are given in equation (25).

Material Matrix of an Anisotropic Material The material matrix is given in equation (32) for an anisotropic material, which has different material properties in each direction [15]:  Ex′ Ex′ v xy  [D *] =  Ex′  sym  where,

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Ex is the modulus of elasticity in the x-direction, Ey is the modulus of elasticity in the y-direction, νxy and νyx are the Poisson ratios. The laminated composite has fibers, which can be located at different angles in each lamina. The material matrix is calculated with considering the fiber angles, and the transformation matrix, which includes the angles, is given in equation (34).  cos 2 β sin2 β 0.5 sin 2 β    [R*] =  sin2 β cos2 β −0.5 sin 2 β  − sin 2 β sin 2 β cos 2 β  

The material matrix for each lamina is calculated through the equation (35). [Dk] = [R*]T [D*][R*]

Finite Element Model As mentioned before, there might be a couple effect between flexural and in-plane strain relation. However, in non-axisymmetric materials, the value of couple effect matrices are equal to zero [15]. In this study, both of the axisymmetric ([0 45 –45 0]), and the non-axisymmetric ([0 0 0 0], [90 0 0 90] and [0 90 0 90]) orientations are examined. The element stiffness and the mass matrices are obtained with flexural and in-plane matrices are combined. The relation is given in equation (36). [In plane, K2 & M2]8×8 + [Flexural, K1 & M1]12×12 = K 21  8 × 8 K 12 × 12 20 × 20  12

This model has six degrees of freedom at each node, but the stiffness and mass matrices have (20×20) size. Although these matrices contain the effect of θz, the last value has to be added formally into matrices for the drilling effect. This relation is given in equation (37). 0  8 × 8 + K e , Me =  12 × 12 20 × 20  0 8 × 8

The value of θz has taken absolute value that 1/1000 of the minimum value in (20×20) stiffness and mass matrices [17]. The generalized displacement vector of node j for local reference coordinates is now expressed as: qj = [uj vj wj θx θy θz ] j

Curvilinear Model The curvature of the deflection curve of a plate under loading at any point depends on the magnitude of the bending moment at that point under the assumption that the material of the plate remains linearly elastic. In many such applications, the deformation can be analyzed by using a thin plate (or beam) theories [18]. The curvature form of the plate is shown in Figure-2a and Figure-2b. Q represents the uniformly distributed load. Ang et al. [19] and Ozturk [20] obtained the deflection equation, and it is given in equation (39). z(x ) =

2. EI  3

12 2 5  x7 6 3 4  − 7 + lx − 5 l x + 2l x  +  3  P  5  x11 40 − + lx10 − l 2 x 9 + 10l 3 x 8 8  2 EI   11 9 16 80 − l 4 x 7 + l 5 x 6  +  7 3

Eef represents an effective modulus of elasticity, which is given in equation (40). Eef =

P represents force, I represents the moment of inertia for the beam, and l represents the projected length of the beam after the force has been applied. P is obtained from equation (41). P=Qb

b represents the width of the plate. The projected beam length l is still an unknown. To obtain the curved length, the projected length l can be calculated from the knowledge of the plate length a using the following formula [20]. l

A deflection curvature is needed to create a curvilinear model for the pre-stressed laminated composite plate. Therefore, rotation angles have to be found to create a

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Figure 2. The curvilinear form of the thin plate. (a) perspective view, (b) side view and the exact length (l) of the plate, (c) two side fixed curved model model.

curvilinear pattern. These rotation angles are obtained by taking the first derivative of equation (39) respect to x. Afterward, these values are placed into a transformation matrix, which is given in Table-1.

It should be noted that the curvilinear plate is modeled by the assembly of a rotated finite element with respect to the y-axis by using the transformation matrix. Rotation of the coordinate system is given in Figure-3. The stiffness and mass matrices can be obtained through the transformation matrix is given in equation (43). Kr = TT × Ke × T Mr = TT × Me × T

where Ke and Me represent stiffness and mass matrices of the reference plate, then rotated matrices Kr and Mr can be obtained with the multiplication of transformation and reference matrices. By assembling the rotated element matrices, global K and M matrices are obtained. Initial Stress The curved plate has initial stress because of the vertical load. After the creation of the curved form, the stress effect must be added to the model. The work done by the in-plane

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Figure 3. Rotation of the coordinate system. (a) the rotation of the axes of translation, (b) the rotation of the axes of rotation.

loads when they move through the bending displacements is given in equation (44) [21]. 2

2. A  dx 

where Nx is an initial axial force acting in the local x direction, and it is given in Figure-4. As seen in Figure-4, axial force Nx is obtained as:

where Fy is a total distributed vertical load. These forces can be obtained from statics. If equation (44) is written in the matrix form: 1 V = qSq 2

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Table 2. Comparison results between the present model and the ANSYS model λ1 C1

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Stiffness matrix Ke, mass matrix Me, and initial stress matrix Se of each plate element are used to form global stiffness, mass, and initial stress matrices. The dynamic response of a plate for the total system can be formulated via equation of motion (47).

In this study, large deflected pre-stressed laminated composite thin plates that are fixed on both sides are examined, as shown in Figure-2c. The use of the four-node quadrilateral plate element approach is convenient since the radius of the curvature of the pre-stressed curved plate is large. The local coordinates of the plate element are transformed into global coordinates using its angle of rotation.

In the scope of the study, four different configurations are examined. For simplicity, the configurations are denoted by C1 ([0 0 0 0]), C2 ([90 0 0 90]), C3 ([0 45 -45 0]) and C4 ([0 90 0 90]). C2 is called cross-ply laminate, C3 is called angleply laminate. In addition, the plate is discretized into 16x16 finite elements. While determining the number of finite elements, the mesh density graph is used for considering the first natural frequency parameters, which is shown in Figure-5. The material properties are: Ex=45 GPa, Ey=12 GPa, νxy= 0.3, ρ = 2080 kg/m3 [22]. The dimension parameters are: a = 0.8 m, b = 0.8 m and h = 5 mm. The glass-fiber epoxy consists of four laminae, therefore each lamina has a 1.25 mm thickness. The effects of vertical non-dimensional load parameter (β) on the natural frequency parameter (λ) and mode shapes for the first five modes are investigated. The non-dimensional load parameter β is considered between the no-load condition and the load parameter value that corresponding to deflection in the free end up to %25 of the plate length value of each composite. λ and β are given in equation (49), where ω is natural frequency, A is the crosssection area of the plate.

Figure 7. Effect of load parameter on the first five natural frequency parameters for C1.

Figure 9. Effect of load parameter on the first five natural frequency parameters for C3.

Figure 8. Effect of load parameter on the first five natural frequency parameters for C2.

Figure 10. Effect of load parameter on the first five natural frequency parameters for C4.

Natural frequencies can be obtained via eigenvalue equation as, ([K]+[S]) – ω2[M] = 0

Results

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

After the plate is deflected by the vertical distributed load, the large deflected plate is fixed at the loading point. Thus, the pre-stressed curved plate is obtained, and free

vibration analysis is performed for all calculated geometric forms. The results of the present study are compared with ANSYS to verify the reliability and validity of the present model. A flow chart of the ANSYS software process is given in Figure 6. Table-2 gives the comparative results between the present model and ANSYS for the first five natural frequencies for both flat and curved plates. Table-2 shows that the natural frequency parameter results of the present model and the error rates with respect

Figure 11. First mode shapes of each configuration under different load parameter value. (a1, a2) C1, (b1, b2) C2, (c1, c2) C3, (d1, d2) C4.

Figure 12. Second mode shapes of each configuration under different load parameter value. (a1, a2) C1, (b1, b2) C2, (c1, c2) C3, (d1, d2) C4.

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Figure 13. Third mode shapes of each configuration under different load parameter value. (a1, a2) C1, (b1, b2, b3, b4) C2, (c1, c2) C3, (d1, d2, d3) C4.

to ANSYS results are quite consistent for each laminated composite configuration. For each composite material, the effect of the load parameter on the natural frequency parameter is examined for the first five natural frequencies. Figure-7 to 10 show the effect of the load parameter on the natural frequency parameter for C1, C2, C3, and C4 composites, respectively. The first two natural frequency parameters behave similarly in all four composites. In the third, fourth, and fifth natural frequency parameters, the natural frequency parameter curves behave differently for all four composites. The curve between the natural frequency parameter and the load parameter differs for each natural frequency on each configuration of the laminated composite. With the increase of the load parameter, the curvature of each

natural frequency parameter shows similarity but progresses by showing different behaviors. Different curves or changes in curvature for each configuration of laminated composites reveal different mode shapes. The first two natural frequency parameter curves are the same for all stacking orders. However, the effect of stacking order on dynamic characteristics becomes more apparent with the third parametric natural frequency. This situation can be observed much more easily especially with the increase of curvature. The first, second, and third mode shapes of each configuration, which are affected by the load parameter value between zero to the max, are shown in Figure-11 to 13, respectively. Different mode shape is not observed in any configuration for the first two natural frequency parameter values.

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Figure 14. Fourth mode shapes of each configuration under different load parameter value. (a1, a2, a3) C1, (b1, b2, b3, b4, b5) C2, (c1, c2, c3) C3, (d1, d2, d3) C4. However, two extra mode shapes occur in the C2 configuration, and one extra mode shape occurs in the C5 configuration for the third natural frequency parameter value. C1 and C3 have the same mode shapes throughout all load values. The fourth mode shapes of each configuration are shown in Figure-14. While three different mode shapes occur for each configuration, two extra mode shapes occur in the C2

configuration due to the curve form difference, as also can be seen in Figure-8. This situation may mean that C2 is a more unstable structure according to other examined composites. The fifth mode shapes of each configuration are shown in Figure-15. Four different mode shapes occur in C1 and C3, three different mode shapes occur in C2, and two different mode shapes occur in C4. When the load parameter is equal to

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

Figure 15. Fifth mode shapes of each configuration under different load parameter value. (a1, a2, a3, a4) C1, (b1, b2, b3) C2, (c1, c2, c3, c4) C3, (d1, d2).

zero, the same mode shapes are seen for all configurations in the first three natural frequencies. However, with the load parameter increase, different numbers of mode shapes can be observed in all configurations. This proves that plate curvature has different effects on natural frequencies due to the change of the load parameter. Within the scope of the study, a total of thirteen different mode shapes are seen in C1, and C3 configurations for a maximum loading condition corresponding to a deflection of %25 of the plate length, while a total of 16 different mode shapes are seen for C2, and 12 different mode shapes are seen for C4. When the load parameter value is equal to zero, the same mode shape occurs in all configurations. With the

increase of the load parameter, more mode shapes may be observed especially after the first two mode shapes.

Conclusion

In this study, the effect of the large deflection on the natural frequency of the pre-stressed laminated composite curved thin plate is investigated. In order to see the effect of the large deflected model on the natural frequency parameter, four different configurations of the laminated composite plate are examined. According to the numerical and graphical results, the following conclusions are drawn:

Sigma J Eng Nat Sci, Vol. 40, No. 4, pp. 772–786, December, 2022

• The theoretical approach of the laminated composite curved plate is performed in different configurations, and the approach gives substantially accurate results with ANSYS. • It is seen that the large deflection equation, which has been used for beams before, can also be used in composite thin plates. • When the load parameter value is equal to zero, the same mode shape occurs in the first three natural frequencies for all examined stacking order configurations. • As the load parameter value increases, the natural frequency parameter also increases. However, this happens in different forms of curvature for each natural frequency due to the effective modulus of elasticity. • The main reason for the change in the parametric natural frequency and mode shapes is that the effective modulus of elasticity and stiffness differ due to the stacking order of the composite plates. • Each change of slope in curvature in β - λ graphs causes a different mode shape. • In different laminated composite configurations, different mode shapes may occur in different order and numbers, especially in higher frequencies. • The fiber orientation is the decisive factor in the fact that plates show differences at different points as the load parameter value increases. • As the plate curvature increases, different mode shapes can be seen for different stacking order configurations, regardless of the plate geometry. This clearly shows the variability of the effect of plate curvature on vibration and natural frequency due to the load parameter.

Share and Cite