Geometrically Nonlinear Analysis of Laminated Composite Plates subjected to Uniform Distributed Load Using a New Hypothesis:  the finite element method (FEM) Approach

D. P. Bhaskar ^a,*, A. G. Thakur ^b

^{a, b} Department of Mechanical Engineering, Sanjivani College of Engineering, Koprgaon,423 603, Savitribai Phule Pune University, Pune, India

1.     Introduction

In all laminated composite plate (LCP) structures, a certain amount of nonlinear behaviour is shown. Once exposed to large displacements and rotations, the plate structure shows geometric nonlinearity. In such cases, linear FE models do not sufficiently estimate the response of structures . Therefore, the expansion of effective and precise nonlinear FE models becomes vital. However,  many researchers have concentrated primarily on linear analysis of composite plates. This study discusses the development of FE models for linear and geometrically nonlinear bending of laminated composite sheets using a new kinematic function for the higher-order shear deformation theory (HSDT).

In the past few decades, several hypotheses have been proposed by well-known researchers. The novel theory of plate by Kirchhoff [1] is recognized as the classical theory of plate (CPT). However, this theory is found to be insufficient in providing realistic responses to deflections since, it does not take into account the effect of transverse shear deformation. In further developments in plate theories, Mindlin [2] developed the first order shear deformation theory (FSDT), wherein the SCF is multiplied to shear modulus to improve the variation of shear strains and stresses through the plate thickness. In addition, the determination of SCF is tedious as it depends on loading conditions, BCs, geometrical and material parameters. Correspondingly, FSDT does not fulfil traction free BCs at the top and bottom faces of the plate. Furthermore, various higher order shear deformation theories (HSDTs) have been proposed to overcome the limitations of CLPT and FSDT and one of the innovative HSDT was developed by Reddy[3] and three-dimensional (3D) elasticity theory for bidirectional bending of laminated composite and sandwich by Pagano [4] and Zenkour [5]. It is understood that 3D elasticity solutions were computationally tough for geometrically nonlinear analysis and led to the development of approximate theories for the analysis of plate structures. The analytical solution of governing equations in the case of geometrically nonlinear analysis of plates becomes difficult after the inclusion of nonlinear sense of Von Karman. Moreover, very few such analyses are available with SS-BCs. FEM by Panda and Natarajan[6] was used for analysis of laminated composite plates in which stress-resultants were computed instead of shear stresses. The two new displacement-based quadrilateral plate elements were suggested by Sung and Kim [7] in FEM studies. The Gauss quadrature scheme of integration was used for calculations of the stiffness matrix. Ren and Hinton [8] extended the Reddy’s HSDT, developed two finite elements (FE), and analyzed laminated composite plates for bending studies. Pandya and Kant [9] presented isoparametric finite elements (FEs) for estimation of inter-laminae stress through the plate thickness using  C⁰ continuity in the displacements. Savithri and Varadan [10] presented a Galerkin method for the solution of displacement based on the higher order theory for the geometric nonlinear analysis of laminated plates. The  HT-FE model is suggested by Qin [11] and studied nonlinear analysis for simplification of coupling detachments between in-plane and out-of-plane displacements. The PVD based least square technique and Galerkin approach are offered by Reddy [12] for analysis by considering deformations and stress resultants as unknown field variables.

Additionally, Reddy [13] presented the effect of E₁/E₂ and a/b (aspect ratio) on the central deflection, in plane and transverse stresses through the thickness of the plate. These remarks [12-13] were used in the present study in the selection of material properties. Goswami and Becker [14] proposed a new non-polynomial HSDT by introducing a new inverse parable and secant kinematic function and achieved remarkable results. Sayyad and Ghugal  [15] proposed bidirectional-bending analysis of plates analytically by Navier’s technique using the trigonometric shear and normal deformation theories. Sayyad and Ghugal [16-18] presented an equivalent single layer trigonometric shear deformation theory for static flexure laminated composite plates. In this study, transverse shear deformation and transverse normal strain effects were considered.

In addition, Fereidoon et al. [19-20] used the power law function [20] to variable modulus of elasticity for analysis of the annular sector wherein extension was made to the Kantorovich method (EKM) and Kirchhoff theory. The methods of polynomial and harmonic differential quadrature (P & HDQ) were used for bending study of functionally graded plates by the varying modulus of elasticity [20]. Ghalebahman [21] carried out research on interlaminar stresses to cross ply laminated composite plates subjected to uniform axial strain by extracting layerwise displacement field by successive iteration. Mantari et al. [22-23] developed a new shear deformation theory [22] and new trigonometric shear deformation theory [23] for the sandwich and composite plates  by suggesting the displacement field depending on the parameter "m". In addition to this Arya et al. [24] presented a zigzag model which used a sine term to represent the nonlinear displacement field. A cosine term was used to represent transverse shear stress and strain. Sayyad and Ghugal [25] presented numerous methods for the analysis of laminated composite and sandwich plates  that could benefit researchers in their work. Kant and Komminen [26] presented a theory that accounted for parabolic distribution of the transverse shear strains through the thickness of the laminate and higher-order terms in Green’s strain vector in the sense of von Karman. A simple C⁰  FE formulation was presented with a total Lagrangian approach and a 9-nodes Lagrangian quadrilateral element was chosen with 9 DoFs/node. For the study of buckling and free vibration of isotropic and FG sandwich beams, Rakočević & Popović [27] analyzed the SS rectangular LCP using a novel computational process based on the layerwise theory of Reddy and compared it to the results obtained through the application of three different finite element models of the ANSYS software. Mallek et al. [28-30] used an efficient nonlinear shell element to analyze multi-layered shells in modified FSDT without SCF [28]. An embedded piezoelectric layer using 3D-shell model based on a discrete double directors shell element was used to lay up shell structures [29]. An original attempt was made to conduct a geometric nonlinear analysis of functionally graded carbon nanotube reinforced composite (FG-CNTRC) structures, with surface-bonded active layers, based on the Kirchhoff shell theory[30]. Mellouli et al. [31-32] proposed a geometrically nonlinear meshfree analysis of 3D shell structures using the double director shell theory with finite rotations. The radial point interpolation method (RPIM) was employed for the construction of shape functions. The discrete system of equations was obtained by incorporating the interpolations into the weak form. Nguyen et al. [33-36]  developed a novel generalized quasi-3D and refined plate theory (RPT) in different studies based on the new transverse shear function to involve the normal deformation in the displacement field and mathematically represented all existing transverse shear functions of HSDT models by a unique polynomial formulation. The weak form was constructed by using the principle of virtual displacements and the variational approach was proposed to derive the strong forms. The numerical results were obtained by using the modern isogeometric analysis. Nguyen et al. [37] proposed a new kinematic function in developing refined hypothesis. This kinematic function satisfies zero transverse stress conditions at the top and bottom face of the plate. In the work reported in the present paper, the new HSDT model using the kinematic function of Nguyen et al. [37] and the C⁰ continuous element developed in the reference [26] have been employed for analysis of LCPs.

The key novelty of the present study is to investigate the new non-polynomial type shear deformation theory in comparison with existing shear deformation theories. The theory involves inverse trigonometric function in the different in-plane directions. The theory considers the effect of trans-verse shear deformation with five unknowns. The  linear and Von Karman’s non-linearity based models inculcate the finite element method approach for analysis . The FE formulations allow the symbolic form to be coded in MATLAB [37]. The numerical results are presented for uniformly distributed load. Non-linear static deflection under uniformly distributed load is limited in the literature(UDL).

The numerical results are computed for different stacking sequences of laminates and compared with similar studies. Few cases are also solved by the FE software ANSYS for central displacement and stresses in LCPs.

2.     The Mathematical Formulation

The LCP illustrated in Fig.1 consists of N number of layers with dimensions (a × b × h) in the Cartesian coordinate system x-y-z. The boundaries of the complete plate coincide with the lines x=0-a and y=0-b. All the layers of elastic material are perfectly united on surface area. The plate is subjected to transverse load q (x, y) on the top surface of the plate (z = h/2) in the form of UDL.

2.1. The displacement fields

According to current HSDT, the axial displacements are u and v along the x and y directions, respectively, and the transverse displacement, w along the z-direction at any point in the plate written in the functional form are shown in Eq. 1. [7]:

Fig. 1. Schematic diagram of  LCP

where the is a shape function which determines the variation of transverse shear stress through the thickness of the plate and u₀, v₀ and w₀ are the mid-surface in-plane displacements along x, y, and z axes, correspondingly while  being the rotations of the transverse normal about y and x axes correspondingly. The form of the displacement field in Eq.1 allows reduction of the 3-D problem to a mid-plane 2-D problem. Once u₀, v₀, w_0, are identified, the displacements of random points (x, y & z) in the continuum can be determined using Eq.1. The von Karman strains associated with the displacement field in static loading are computed using the Green-Lagrange strain-displacement relations, in the view of the small strain but moderate rotations assumption.

The state of strains at any point in the plate is defined by the in-plane strains and transverse shear strains as shown in Eq.2. [7]. Here the suffix ‘L’ and ‘NL’ stand for the linear and nonlinear part, respectively. The strains and are the normal strains vector and transverse shear strain vector correspondingly, while is an in-plane strain.

2.2. The Kinematic Function

The kinematic functions mentioned in CLPT, FSDT, and Reddy’s HSDT and a new  one [27] for present inverse trigonometric shear deformation theory (ITSDT) is shown in Eq.3. The present kinematic function satisfies the requirement of shear stresses to be zero on the top and bottom surfaces.

2.3. The Constitutive Relations

The relationship needed between the stresses and strains [40] in the directions of the orthotropic axes, in the typical lamina of the laminate is shown in Eq.4.

where Q_ij contains the material stiffness coefficient [17] shown in Eq.5

In this  are the Young’s moduli in the fiber orientation (direction-1) and in-plane normal (direction-2) directions respectively. And  are the shear moduli in the planes respectively, with Poisson’s ratios . There are four independent material properties:  and their reciprocal relations are specified by . The stress-strain relations are accordingly decisive and are the foundation for the stiffness and stress analysis of an individual lamina subjected to UDL in its own surface.

3.     Equilibrium Equations

The summation of the total internal work ( ) and total external work ( ) in accordance with the principle of virtual work (PVW) is written as shown in Eq.6. In further interpretation  is shown in Eq.6.

where, v is the undefined volume while and  are the stress vector and virtual strain vector due to virtual displacement ( ) respectively. And q is the transverse load on the top of the plate.

3.1. Stress resultants

The stress resultants computed by integrating stresses through the thickness direction of the plate [8] are shown in Eq. 7.

Solving Eq.7, force-strain relation and moment-curvature relationships [8] are obtained, shown in Eq.8.

where:

N =In-plane force resultants: (Nx, N_y and Nxy), M=Bending moments (Mx, My and Mxy), A=Extensional stiffness matrix, relates in-plane forces of the in-plane strains, B=Coupling stiffness matrix, which couples the forces and moments to the mid-plane strain-curvature, D=Bending moment stiffness matrix., (Relates bending moments to the plate curvature),  =Membrane strains. Eq. 9., k = Curvature strains. Eq. 9.

where i, j=1, 2, 6

The total virtual work equation [3] is then written in Eq.10.

where, P_x and P_y are the stated global edge tractions in the x and y directions and q is the transverse UDL.

3.2. The finite element models

FEM is a powerful numerical technique for the solution of differential and integral equations in the geometrically complex behaviour of composite structures. Firstly the LCP domain is discretised in sub-domain using eight nodes serendipity elements shown in Fig.2. Further the FE model of plate is developed using PV displacements for calculation of deformations.

The C⁰ continuity is used in the present theory for its FE approximation. Adopting shape functions (Ni), the displacement vector d within the finite element given as a function of eight discrete points is given as:

The shape functions (Ni) and their derivatives, which are needed for the computations of strains. The shape functions N₁ to N₈ are computed based on local coordinates at specific node and are shown in Eq.12 [42-43].

Fig.2. Eight nodes serendipity element

The nodal displacements within the element in-summation notations are given by Eq.13.

Incorporating Eq. 13 in Eq. 10, PVW equation is acquired in terms of shape functions and global displacement. [43]

B₀=In-plane components, B_b=Bending components

BS= Transverse shear components, W_s = Transverse load components, W_e=In-plane edge load components.

The Newton-Raphson method (NRM) is adopted to solve the assembled nonlinear equilibrium equations (Eq.14). The solution is initialized using residual equation, Eq. 15, by initial guess value , while  being tangent stiffness matrix:

The convergence of the solution is achieved for assuming three decimal accuracy [41-42].

Gaussian Quadrature numerical integration respective three-point (3P) method and two-point (2P) method are used for bending and shear term solution of general form are shown in Eq.17 [41].

while the 2P and 3P wightages are, 2P: N=2, (xN, k): ±0.05773 & (wN, k):1, 1. 3P: N=3, (xN, k): ±0.7754, 0 & (wN, k):0.555, 0.8888.

4.     Numerical Results and Discussions

In the present investigation, central displacements and stresses are presented for symmetric and anti-symmetric cross-ply laminates shown in Fig. 3. They consist of graphite-epoxy lamina subjected to the UDL using the present hypothesis-based FE model derived in the preceding sections. The formulation and accuracy of the present FEM-MATLAB code is verified with the closed-form solution of Zenkour [5] and other diverse theories in the literature.

Following material properties [3] of Graphite-epoxy layers are considered for the analysis.

Similarly, geometrical information on the overall plate is taken as:  a=b=1, S=4 to 100. The symbols have their usual meaning.

Firstly, linear bending analysis of plate was performed, followed by geometrically nonlinear bending analysis.

Loading: The given transversely load (UDL) as shown in Fig.4 is expanded in a double Fourier series shown in Eq.18 [25].

where qo denotes the load intensity at the center of the plate and for UDL, m, n = 1, 3, 5….

4.1. Linear bending behavior

In a linear bending analysis relation between applied forces and displacements is linear while the stiffness matrix of LCP structure is constant. The non-dimensional transverse displacement and stresses quantities [26] shown in Eq.19 are used for presenting the results in graphical and tabular form:

Fig.3. Configuration of plate laminates

Table 1. Boundary Conditions (SS-BCs) [18]

Fig. 4. Square LCP subjected to UDL

An assessment of non-dimensional transverse displacement and stresses are displayed in Table 2 for (0°/90°) antisymmetric LCP of equal-thickness lamina subjected to UDL.  The present results are compared with those stated by Sayyad and Ghugal [15], Reddy [3], Mindlin [2], Kirchhoff [1] and 3D elasticity solution of Zenkour [5]. The investigation of Table 2 reveals that the present results are in exceptional agreement with a 3D elasticity solution [19]. It is also to be  noted that the present results are improved even more than the renowned theory by Reddy [3], FSDT by Mindlin [2] and CLPT by Kirchhoff [1].       The CLPT of Kirchhoff[1]  predicts the results as transverse shear deformation is ignored. Similarly, Table 3 shows a comparison of non-dimensional central displacements and various stresses for (0°/90°/0°) symmetric LCP of equal-thickness lamina subjected to UDL. It is observed from Table 3 that the present kinematic function computes outstanding results compared to those presented by Zenkour [5] for symmetric lamination.

The results of non-dimensional central displacements and stresses are presented in Table 4 for (0°/90°/0°/90°) antisymmetric laminated composite plates of the equal - thickness lamina with S=4, 10, 20, 50 & 100. The comparison of the results with previous studies [1-3, 15] and particular  one [5] as shown in Table 2-4, shows the quality of the present results . From the convergence analysis, it is clear that the outcome differs from other outcomes[1-3, 15] but up to 4 percent proximity to the Eaxct [5]. The results reported here are also in excellent agreement with the 3D elasticity solution presented by Zenkour [5]. It is also observed that the present kinematic trigonometric function shows improvement over other functions shown in Eq.3 and those expressed by Sayyad and Ghugal [15].

From the results shown in Table 2-4, it has been noted that results of non-dimensional transverse displacements  for (0°/90°) at S=4 are quite close to results obtained by Zenkour [5] which is not the case when is compared. This is a fact that non-dimensional transverse displacement is independent on thickness parameters and  varies through the thickness. This is also due to the FEM discrete technique and does not give a continuous solution. The average disparity in the present results is calculated in the form of % Error as shown in the Tables 2-4  compared to the Zenkour [5] results. The results show more accuracy of displacement outcomes than stress outcomes. There are two kinds of plate theories: correspondingly based on displacement and stress. In  the current study, displacement is assumed in the form of polynomial, while in the latter, stress is assumed in  the form of a polynomial. The central displacements are precisely estimated as the  present theory is based on displacements. Although the results of stresses are even more noticeable in some cases. The CPT [1] it gives is deprived of results due to the assumption in the displacement field sides. Figsures 5-14 shows the variation of non-dimensional displacements and stresses across the thickness of the LCP.

The non-dimensional stresses are plotted along the x-axis and plate thickness (z/h) along the y-axis for symmetric and antisymmetric LCPs with schemes (0°/90°), (0°/90°/0°), (0°/90°/90°/0°), (0°/90°/0°/90°) and aspect ratio, S=4 & 10. It can be seen that the results predicted by the present theory are almost the same as that of Zenkour [5].

Fig. 5Through thickness variation of non-dimensional in-plane normal stress ( ) for (0°/90°) antisymmetric LCP

Fig. 6 Through thickness variation of non-dimensional transverse shear stress ( ) for (0°/90°) antisymmetric LCP

Fig. 7Through thickness variation of non-dimensional in-plane shear stress ( ) for (0°/90°) antisymmetric LCP

Fig. 8Through thickness variation of non-dimensional in-plane normal stress ( ) for (0°/90°/0⁰) antisymmetric LCP

It may be easily inferred from Figs. 5-17, that all the variations of stresse have similar behaviours as specified by Sayyad and Ghugal [15] and the exact solution of Zenkour [5].

It can be seen in Fig. 5 that there are distinct discontinuities at the midplane of the plate. The stresses along the fibre axis is more than that of seen in the perpendicular direction of the fibre axis as an effect of E₁/E₂=25 (High elastic modulus along the fibre).

The Figs. 6-14 show out of plane stresses and those accurately satisfy the required boundary conditions at the free (lowermost) and loaded (top) surfaces. The  top surface of plate encounters compressive stress (- ) while the bottom face experiences positive stress (+ )  shown in Fig 8.

As seen from Figs. 5-16, the out-of-plane transverse normal and shear stresses accurately satisfy the continuity conditions at the interfaces. As an exception, as is noticed in Fig. 6,  the maximum value of transverse shear stress occurs at the middle surface of the plate. In any particular case, the stresses in the specific lamina are more for S=4 than for S=10. From  observing  Fig. 8, it is seen that in-plane stress ( ) variation is linear and parabolic, respectively, for S=4 and S=10. That also applies to the rest of the cases .

By analyzing the bending characteristic of LCP the analysis of these results determines the performance and efficacy of this approach.

These analyses show that non-dimensional transverse displacement decreases in the sequence of stacking in  order of antisymmetric (0⁰/90⁰), symmetric (0⁰/90⁰/0⁰), symmetric (0⁰/90⁰/90⁰/0⁰) and antisymmetric (0⁰/90⁰/0⁰/90⁰) square LCPs.

Fig. 9Through thickness variation of non-dimensional transverse shear stress ( ) for (0°/90°/0⁰) symmetric LCP

Fig. 10Through thickness variation of non-dimensional transverse shear stress ( ) for (0°/90°/0⁰) antisymmetric LCP

Fig. 11Through thickness variation of non-dimensional in plane normal stress ( ) for (0°/90°/90°/0°) symmetric LCP

Tables 3-5 and Figs. 5–14, show that the results from the present elements agree well with the analytical solution and the exact solution results. The Figs. 15-16, validate non-dimensional central displacement  versus x/a subjected to UDL for S=4 & 10.

Incidentally, the displacements vanish along the edge, as does the need for BCs shown in Figs.17-18. There is not much variation in central deflection results found for S>>25 and therefore the increase in aspect ratio, S is redundant.

In addition to analyses based on FE MATLAB , the  structural analysis is carried out using  ANSYS workbench  to  solve linear or nonlinear problems by  modeling LCPs by solid element. ANSYS provides ANSYS composite Prep-Post to facilitate building finite element models (FEMs) and access results. The SOLID185 element is used for the linear case.  The linear bending analysis of LCPs carried out in the ANSYS environment is shown in Fig. 19-25. The plate is discretized using Solid 185 finite elements.

The Hashin failure criteria are used in the 2-dimensional approach to calculate stress at nodes. The contour plots for linear bending analysis of LCP analyzed in ANSYS environment is shown in Figs. 19-23.

Fig. 12Through thickness variation of non-dimensional in plane normal stress ( ) for (0°/90°/90°/0°) symmetric LCP

Fig. 13Through thickness variation of non-dimensional transverse shear stress ( ) for (0°/90°/90°/0°) symmetric LCP

Fig. 14 Through thickness variation of non-dimensional transverse shear stress ( ) for (0°/90°/90°/0°) symmetric LCP

Fig. 15 Variation of non-dimensional central deflection ( )  of LCP for S=4.

Figure 19 shows transverse deflection  for S=10 whereas by FEM-MATLAB  (present) and by Zenkour [5]. The percentages of error while compared with Zenkour [5] are 1.2422 % (Present ANSYS) and 1.2422 % (Present EEM-ATLAB). Fig.21 shows the non-dimensional shear stresses  computed by ANSYS as 0.1090 and by FEM-ANSYS codes as 0.1029.

Thus the results of the formulations presented are validated with those of the ANSYS analysis and there is a reasonably good agreement between them. The different cases shown in Figs 20-23 are simulated in the ANSYS. From these results it can be observed that ANSYS also provides a reasonably accurate solution with less computational effort.

Fig.16. Variation of non-dimensional central deflection ( )  of LCP for S=10

Fig.17 The descretized (0°/90°)  antisymmetric LCP using Solid 185 Element

Fig. 18. The SS-BCs (0°/90°) antisymmetric  LCP subjected to transverse load (UDL)

Fig. 19. The non-dimensional central displacement  of the (0°/90°)  antisymmetric LCP  using ANSYS

Fig. 20. The non-dimensional inplane shear stress  in the (0°/90°) antisymmetric LCP using ANSYS

Fig. 21. The non-dimensional central  displacement of the plate of  (0°/90°/0⁰) symmetric LCP using ANSYS

Fig. 22. The non-dimensional in-plane stress  in the (0°/90°/0⁰)  symmetric LCP  using ANSYS

4.2. Geometrically nonlinear behavior

Having recognized the reliability of the FE model built for the linear bending, the analysis of LCPs drawn-out for geometric nonlinear analysis is performed by considering nonlinearity terms. In nonlinear analysis, the relation between applied forces and displacements is a nonlinear one. The geometric nonlinear effects which originate as deformations are large. The stiffness matrix of the plate material varies with changes in configuration during the load process.

The non-dimensional quantities for deflection and stresses are shown in Eq.20, and used to present the results in graphic and tabular form [26]. For the study, a plate with the aspect ratio S= 20 and 40 is considered.

The equilibrium equations are solved at each load step by using a Newton Raphson (NR) method with a convergence tolerance for nodal displacements by solving nonlinear FE equations.  In all the examples, loads 0 to 250 are considered with load increments of 50.

A convergence study is conducted and presented in Table 5. Based on this convergence study, it is concluded that an 8×8 mesh is sufficient.  Table 5 clearly shows that the performance of the present FE formulations are very desireable in terms of solution accuracy. The central displacements are noteworthy for this mesh size (8×8). The present results are compared with FEM results of Kant & Kommineni (26) and with Ron & Hinton (8).

The results are also computed for few nonlinear bending cases in ANSYS environment. The SOLID186, quadratic serendipity element is used in the nonlinear case. The present modification is to improve the stress predictions within these finite elements. For (0°/90°/ 0°/90⁰)antisymmetric LCP with aspect ratio, S=20 results are  shown in Fig. 26. The non-dimensional transverse displacement,  shown in Table 5, by FEM-MATLAB (present), Kant [26], Ren [8] and ANSYS (present) are 0.7582(-1.057% Error), 0.7520(-1.866 % Error), 0.7663 and 0.780 (1.787% Error) respectively. The comparison of the results in Table 5  show  that in the geometrically nonlinear analysis  ANSYS also gives results with reasonable agreement.

Table 6 and 7 show comparisons of geometrically linear and nonlinear displacements obtained for (0°/90°/0°/90°) antisymmetric LCP. In nonlinear bending analysis, the coupling between lamina increases to plate stiffening due to the von Karman nonlinearity and thus decreases in central displacements compared to linear bending analysis. This is because in the case of the linear plate central displacement relies on bending and the same in the case of geometrically nonlinear on both extension and bending.

Fig. 23. The non-dimensional central displacement  of the (0°/90°/90⁰/0⁰) symmetric LCP using ANSYS

Fig. 24. The descretized (0°/90°/0⁰/90⁰)  antisymmetric LCP using Solid 186

Fig. 25. The non-dimensional transverse displacement  of the (0°/90°/0⁰/90⁰)  LCP  using ANSYS (S=20)

Fig. 26 Load-deflection diagram at specific load points in geometrically nonlinear analysis.

The analysis of the results in Table 6-7, shows that the findings reported are in acceptable agreement with the results of Kant and Kommineni [26] and with Y.X. Zhang and Kim [7] and other researchers. The geometrically nonlinear analysis creates central displacement and stress solutions that are more realistic. Fig. 26, presents the nonlinear central displacement  in the plate as a function of the load parameter. The displacement turns out to be closely independent for larger loads. As the number of lamina in the plate increases, the plate becomes increasingly stiffer. The non-dimensional central displacements,  initially increases, then increase nonlinearly with each increase in load.  The displacement   significantly depends on the plate aspect ratio (S) at small loads.

5.     Conclusions

Based on the current theory and the nonlinearity of von Karman, an FE model is developed. In this research, inverse TSDT is used on the basis of the modern kinematic method for the geometrically linear and nonlinear analysis of LCPs. The kinematic function satisfies traction free BCs and also it does not for SCF unlike FSDT. The programming codes are written in the MATLAB environment based on FEM formulations and are significant in plate analysis with SS-BCs. The Newton Raphson method is used to solve the nonlinear equations. Numerical results are obtained for different symmetric and antisymmetric lamination schemes. Results are also computed by using FE software ANSYS for limited cases. From the numerical results it is seen that the present new hypothesis calculates  central displacements and stresses more precisely, compared to other HSDTs available in the literature. The authors conclude that this article would pose numerical methodology as a guide for two-dimensional plate theories and methods.

Acknowledgements

The results reported herein were obtained during an research work carried out at Sanjivani College of Engineering, Kopargaon (SCOEK), Savitribai Phule Pune University (SPPU), Pune (India). We are thankful to both SCOEK and SPPU Pune for this opportunity.

L: Linear (L) and NL: Geometrically nonlinear analysis.

References

[1]   Kirchhoff  G. About the balance and movement of an infinitely thin elastic rod. Journal of Pure and Applied Mathematics 1850; 40:51-88.

[2]   Mindlin R. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. ASME Journal of Applied Mechanics 1951; 18:31–38.

[3]   Reddy J. A Simple Higher-Order Theory for Laminated Composite Plates. Journal of Applied Mechanics 1984; 51:745-752.

[4]   Pagano N. Exact Solutions for Rectangular Bidirecional Composites and Sandwich Plates. Journal of Composite materials 1970;4, 20-34.

[5]   Zenkour A. Three-dimensional Elasticity Solution for Uniformly Loaded Cross-ply Laminates and Sandwich Plates. Journal of Sandwich Structures & Material 2007;9:213–238.

[6]   Panda S R, Natarajan R. Finite element analysis of laminated composite plates. International journal for numerical methods in engineering 1979;14: 69-79.

[7]   Zhang Y,  Kim K. Geometrically nonlinear analysis of laminated composite plates by two new displacements-based quadrilateral plate elements. Composite Structures 2006; 72: 301–310.

[8]   Ren J, Hinton H. The finite element analysis of homogeneous and laminated composite plates using a simple higher order theory. Communications in applied numerical methods 1986;2: 217-228.

[9]   B Pandya, T Kant. Flexural analysis of laminated composites using refined higher-order C⁰ plate bending elements. Computer Methods in Appl. Mech. Eng. 1988; 66(2): 173-198 .

[10] Savithri S, Varadan T. Large deflection analysis of laminated composite plates. HT. Journal of  Non-linear Mechanics 1993;28: 1-12.

[11] Qinm Q. Nonlinear analysis of thick plates by HT FE approach. Computers & Structures 1995; 61:227-281.

[12] Reddy J, Roman A, Filipa M. Finite Element Analysis of Composite Plates and Shells. Encyclopedia of Aerospace Engineering 2010.

[13] Reddy B, Reddy R, Kumar S, Reddy K. Bending analysis of laminated composite plates using finite element method. International Journal of engineering, Science and Technology 2012;4:177-190.

[14] Goswami S, Becker W. A New Rectangular Finite Element Formulation Based on Higher Order Displacement Theory for Thick and Thin Composite and Sandwich Plates. World Journal of Mechanics 2013;3:194-201.

[15] Sayyad A S, Ghugal Y. A new shear and normal deformation theory for isotropic, transversely isotropic, laminated composite and sandwich plates. International Journal of Mechanics and Materials in Design 2014;10:247-267.

[16] Ghugal Y, Sayyad A S. Stress analysis of thick laminated plates using trigonometric shear deformation theory. International Journal of Applied Mechanics 2013;5(1): 1350003 1-23

[17] Sayyad A S, Ghugal Y. Flexure of cross-ply laminated plates using equivalent single layer trigonometric shear deformation theory. Structural Engineering & Mechanics 2014; 51(5):867-891.

[18] Sayyad A S, Ghugal Y. On the Buckling of Isotropic, Transversely Isotropic and Laminated Composite Rectangular Plates. International Journal of Structural Stability and Dynamics 2014; 14(6), 1-32.

[19] A Fereidoon, A Mohyeddin, M Sheikhi, H Rahmani. Bending analysis of functionally graded annular sector plates by extended Kantorovich method. Composites: Part B 2012; 43: 2172–2179.

[20] A Fereidoon, M Asghardokht seyedmahalle, A Mohyeddin. Bending analysis of thin functionally graded plates using generalized differential quadrature method. Arch Appl Mech; 2011; 81:1523–1539.

[21] Morteza Rezvani, Ahmad Ghasemi Ghalebahman. Interlaminar stresses in symmetric cross-ply composite laminates using Layerwise theory. Modares Mechanical Engineering 2015; 14(1):56-66.

[22] Mantari J L, Oktem A S, Guedes Soares C. A new higher order shear deformation theory for sandwich and composite laminated plates. Composites Part B: Engineering 2012; 43(3), 1489–1499.

[23] Mantari J L, Oktem A S, Guedes Soares C . A new trigonometric shear deformation theory for isotropic, laminated composite and sandwich plates. International Journal of Solids and Structures 2012; 49(1), 43–53.

[24] Arya H, Shimpi R P, Naik N K. A zigzag model for laminated composite beams. Composite Structures 2002; 56(1):21–24.

[25] Sayyad A S, Ghugal Y. On the free vibration analysis of laminated composite and sandwich plates: A review of recent literature with some numerical results. Composite Structures 2015; 129:177–201.

[26] Kant T, Kommineni J R. C⁰ finite element geometrically nonlinear analysis of fibre reinforced composite and sandwich laminates based on a higher-order theory. Comput Struct 1992;45:511–20.

[27] Rakočević M, & Popović S.Bending analysis of simply supported rectangular laminated composite plates using a new computation method based on analytical solution of layerwise theory. Archive of Applied Mechanics 2017; 88(5): 671–689.

[28] Mallek H, Jrad H, Wali M, Dammak F. Geometrically nonlinear finite element simulation of smart laminated shells using a modified first-order shear deformation theory. Journal of Intelligent Material Systems and Structures 2019;30(4): 517-535.

[29] Mallek H, Jrad H, Wali M, Dammak F. Piezoelastic response of smart functionally graded structure with integrated piezoelectric layers using discrete double directors shell element. Composite Structures 2019;210:354-366.

[30] Mallek H, Jrad H, Wali M, Dammak F. Geometrically non-linear analysis of FG-CNTRC shell structures with surface-bonded piezoelectric layers. Computer Methods in Applied Mechanics and Engineering Comput. Methods Appl. Mech. Engrg. 2019; 347:679–699. 

[31] Hana Mellouli, Hanen Jrad, Mondher Wali, Fakhreddine Dammak. Geometrically nonlinear meshfree analysis of 3D-shell structures based on the double directors shell theory with finite rotations.Steel and Composite Structures 2019;31(4): 397-408.

[32] Mellouli H, Jrad H, Wali M, Dammak F. Meshfree implementation of the double director shell model for FGM shell structures analysis. Engineering Analysis with Boundary Elements 2019; 99: 111–121.

[33] Nguyen T. N., Thai C H, Nguyen-Xuan, H. On the general framework of high order shear deformation theories for laminated composite plate structures: A novel unified approach. International Journal of Mechanical Sciences 2016;110:242–255.

[34] Nguyen T N, Ngo T D, Nguyen-Xuan H. (2017). A novel three-variable shear deformation plate formulation: Theory and Isogeometric implementation. Computer Methods in Applied Mechanics and Engineering 2017;326:376–401.

[35] Nguyen N-T, Hui D, Lee J,Nguyen-Xuan H. An efficient computational approach for size-dependent analysis of functionally graded nanoplates. Computer Methods in Applied Mechanics and Engineering 2015;297:191–218.

[36] Nguyen H X, Nguyen T N, Abdel-Wahab, M, Bordas S P A, Nguyen-Xuan, H, Vo T P. A refined quasi-3D isogeometric analysis for functionally graded microplates based on the modified couple stress theory. Computer Methods in Applied Mechanics and Engineering 2017;313:904–940.

[37] Roque C. Symbolic and numerical analysis of plates in bending using Matlab. Journal of Symbolic Computation 2014;61-62: 3–11.

[38] Nguyen T K, Truong Phong, Nguyen T, Vo T P, Thai H T. Vibration and buckling analysis of functionally graded sandwich beams by a new higher-order shear deformation theory. Composites Part B: Engineering 2015;76: 273–285. 

[39] Hashin Z. Failure criteria for unidirectional fibre composites, ASME Journal of Applied Mechanics, Vol. 47 (2), 1980, pp 329-334.

[40] Stephen Timoshenko, J M Gere.Theory of Elastic Stability. 1961; McGraw-Hill Book Company.

[41] Ochoa O O, Reddy J N.Finite Element Analysis of Composite Laminates. 1992;Kluwer academic publisher.

[42] Ever J Barbero. Finite Element Analysis of Composite Materials Using ANSYS. 2014;CRC Press.

[43] Nachiketa Tiwari. Courses on Finite element Method. to Indian Institute of Technology, Kanpur(India). NPTEL chapters.