Document Type : Research Paper
Authors
Department of Mechanical Engineering, Sanjivani College of Engineering, Koprgaon,423 603, Savitribai Phule Pune University, Pune, India
Abstract
Keywords

Mechanics of Advanced Composite Structures 7(2020) 271 – 285


Semnan University 
Mechanics of Advanced Composite Structures journal homepage: http://MACS.journals.semnan.ac.ir 
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
keywords 

ABSTRACT 
Laminated composite plate Geometric nonlinearity New kinematic function Uniformly distributed load Finite Element Method

This paper presents a finite element method (FEM) for linear and geometrically nonlinear behaviours of cross ply square laminated composite plates (LCPs) subjected to a uniform distributed load (UDL) with simply supported boundary conditions (SSBCs). The original MATLAB codes were written to achieve a finite element (FE) solution for bending of the plate. In geometrically nonlinear analysis, changes in geometry take place when large deflection exists to consequently provide nonlinear changes in the material stiffness and affect the constitutive and equilibrium equations. The Von Karman form nonlinear strain displacement relations and a new inverse trigonometric shear deformation hypothesis were used for deriving the FE model. Here, inplane displacements made use of an inverse trigonometric shape function to account for the effect of transverse shear deformation. This hypothesis fulfilled the traction free BCs and disrupted the necessity of the shear correction factor (SCF). Overall the plate was discretized using the eightnode isoparametric serendipity element. The equilibriums governing equations associated boundary conditions were obtained by using the principle of virtual work (PVW). The numerical results were obtained for central deflections, inplane stresses and transverse shear stresses for different stacking sequences of cross ply laminates. The results were also computed by the FE software ANSYS for limited cases. The results obtained showed an acceptable agreement with the results previously published. The findings suggested the future use of a new FE model for linear and nonlinear laminated composite plate deformation. 




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 higherorder shear deformation theory (HSDT).
In the past few decades, several hypotheses have been proposed by wellknown 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 threedimensional (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 SSBCs. FEM by Panda and Natarajan[6] was used for analysis of laminated composite plates in which stressresultants were computed instead of shear stresses. The two new displacementbased 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 interlaminae stress through the plate thickness using C^{0} 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 HTFE model is suggested by Qin [11] and studied nonlinear analysis for simplification of coupling detachments between inplane and outofplane 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_{1}/E_{2} and a/b (aspect ratio) on the central deflection, in plane and transverse stresses through the thickness of the plate. These remarks [1213] were used in the present study in the selection of material properties. Goswami and Becker [14] proposed a new nonpolynomial HSDT by introducing a new inverse parable and secant kinematic function and achieved remarkable results. Sayyad and Ghugal [15] proposed bidirectionalbending analysis of plates analytically by Navier’s technique using the trigonometric shear and normal deformation theories. Sayyad and Ghugal [1618] 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. [1920] 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. [2223] 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 higherorder terms in Green’s strain vector in the sense of von Karman. A simple C^{0 } FE formulation was presented with a total Lagrangian approach and a 9nodes 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. [2830] used an efficient nonlinear shell element to analyze multilayered shells in modified FSDT without SCF [28]. An embedded piezoelectric layer using 3Dshell 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 (FGCNTRC) structures, with surfacebonded active layers, based on the Kirchhoff shell theory[30]. Mellouli et al. [3132] 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. [3336] developed a novel generalized quasi3D 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^{0} 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 nonpolynomial type shear deformation theory in comparison with existing shear deformation theories. The theory involves inverse trigonometric function in the different inplane directions. The theory considers the effect of transverse shear deformation with five unknowns. The linear and Von Karman’s nonlinearity 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. Nonlinear 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.
The LCP illustrated in Fig.1 consists of N number of layers with dimensions (a × b × h) in the Cartesian coordinate system xyz. The boundaries of the complete plate coincide with the lines x=0a and y=0b. 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.
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 zdirection at any point in the plate written in the functional form are shown in Eq. 1. [7]:
Fig. 1. Schematic diagram of LCP
(1) 
where the is a shape function which determines the variation of transverse shear stress through the thickness of the plate and u_{0}, v_{0} and w_{0} are the midsurface inplane 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 3D problem to a midplane 2D problem. Once u_{0}, v_{0}, 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 GreenLagrange straindisplacement relations, in the view of the small strain but moderate rotations assumption.
(2) 
The state of strains at any point in the plate is defined by the inplane 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 inplane strain.
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.
(3) 
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.
(4) 
where Q_{ij }contains the material stiffness coefficient [17] shown in Eq.5
(5) 
In this are the Young’s moduli in the fiber orientation (direction1) and inplane normal (direction2) 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 stressstrain 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.
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.
[3]

(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.
The stress resultants computed by integrating stresses through the thickness direction of the plate [8] are shown in Eq. 7.
(7) 
Solving Eq.7, forcestrain relation and momentcurvature relationships [8] are obtained, shown in Eq.8.
(8) 
where:
N =Inplane force resultants: (Nx, N_{y} and Nxy), M=Bending moments (Mx, My and Mxy), A=Extensional stiffness matrix, relates inplane forces of the inplane strains, B=Coupling stiffness matrix, which couples the forces and moments to the midplane straincurvature, D=Bending moment stiffness matrix., (Relates bending moments to the plate curvature), =Membrane strains. Eq. 9., k = Curvature strains. Eq. 9.
(9) 
where i, j=1, 2, 6
The total virtual work equation [3] is then written in Eq.10.
(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.
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 subdomain 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^{0} 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:
(11) 
The shape functions (Ni) and their derivatives, which are needed for the computations of strains. The shape functions N_{1} to N_{8 }are computed based on local coordinates at specific node and are shown in Eq.12 [4243].
Fig.2. Eight nodes serendipity element
(12) 
The nodal displacements within the element insummation notations are given by Eq.13.
(13) 
Incorporating Eq. 13 in Eq. 10, PVW equation is acquired in terms of shape functions and global displacement. [43]

(14) 
B_{0}=Inplane components, B_{b}=Bending components
BS= Transverse shear components, W_{s} = Transverse load components, W_{e}=Inplane edge load components.
The NewtonRaphson 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:
(15) 
The convergence of the solution is achieved for assuming three decimal accuracy [4142].
(16) 
Gaussian Quadrature numerical integration respective threepoint (3P) method and twopoint (2P) method are used for bending and shear term solution of general form are shown in Eq.17 [41].

(17) 
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.
In the present investigation, central displacements and stresses are presented for symmetric and antisymmetric crossply laminates shown in Fig. 3. They consist of graphiteepoxy lamina subjected to the UDL using the present hypothesisbased FE model derived in the preceding sections. The formulation and accuracy of the present FEMMATLAB code is verified with the closedform solution of Zenkour [5] and other diverse theories in the literature.
Following material properties [3] of Graphiteepoxy 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].
(18) 
where qo denotes the load intensity at the center of the plate and for UDL, m, n = 1, 3, 5….
In a linear bending analysis relation between applied forces and displacements is linear while the stiffness matrix of LCP structure is constant. The nondimensional transverse displacement and stresses quantities [26] shown in Eq.19 are used for presenting the results in graphical and tabular form:
(19) 
Fig.3. Configuration of plate laminates
At y, x=0, a 
At x, y=0, b 


Fig. 4. Square LCP subjected to UDL
An assessment of nondimensional transverse displacement and stresses are displayed in Table 2 for (0°/90°) antisymmetric LCP of equalthickness 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 nondimensional central displacements and various stresses for (0°/90°/0°) symmetric LCP of equalthickness 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 nondimensional 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 [13, 15] and particular one [5] as shown in Table 24, shows the quality of the present results . From the convergence analysis, it is clear that the outcome differs from other outcomes[13, 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].
Table 2. Nondimensional displacements and stresses for (0°/90°) antisymmetric LCP 


S 
Hypothesis

% Error 
% Error 
% Error 
% Error 
% Error 


4 
Present 
3.108 
1.583 
1.298 
8.662 
0.177 
11.77 
0.117 
0.244 
0.609 
0.238 
14.48 


Sayyad and Ghugal [15] 
2.998 
5.057 
1.260 
6.444 
0.139 
12.32 
0.110 
0.239 
2.845 
0.239 
14.33 


Reddy[3] 
3.070 
2.767 
1.269 
7.187 
0.131 
17.35 
0.1070 
0.241 
2.032 
0.241 
13.62 


Mindlin[2] 
3.008 
4.743 
1.063 
10.16 
0.125 
20.88 
0.099 
0.191 
22.35 
0.191 
31.54 


Kirchhoff[1] 
1.695 
46.31 
1.076 
9.096 
0.126 
20.18 
0.0934 
– 
– 
– 
– 


Zenkour[5] 
3.158 
– 
1.184 
– 
0.159 
– 
– 
0.246 
– 
0.279 
– 


10 
Present 
1.967 
1.811 
1.142 
5.193 
0.131 
0.816 
0.1029 
0.269 
9.634 
0.259 
4.596 


Sayyad and Ghugal [15] 
1.907 
1.294 
1.105 
1.813 
0.130 
0.538 
0.0978 
0.266 
8.130 
0.266 
7.258 


Reddy[3] 
1.917 
0.760 
1.104 
1.740 
0.127 
2.00 
0.0977 
0.264 
7.317 
0.264 
6.451 


Mindlin[2] 
1.905 
1.397 
1.053 
3.011 
0.126 
2.692 
0.0961 
0.194 
21.13 
0.194 
21.77 


Kirchhoff[1] 
1.695 
12.24 
1.076 
0.893 
0.126 
2.384 
0.0934 
– 
– 
– 
– 


Zenkour[5] 
1.932 
– 
1.086 
– 
0.130 
– 
– 
0.246 
– 
0.248 
– 


Table 3. Nondimensional displacements and stresses for (0°/90°/0^{0}) symmetric LCP 

S 
Hypothesis 
% Error 
% Error 
% Error 
% Error 
% Error 

4 
Present 
2.914 
4.257 
1.050 
6.4475 
0.118 
1.615 
0.109 
0.355 
19.69 
0.443 
8.814 

Sayyad and Ghugal[15] 
2.893 
4.941 
1.034 
7.8902 
0.113 
8.077 
0.109 
0.357 
19.26 
0.435 
10.43 

Reddy[3] 
2.909 
4.425 
1.017 
9.3685 
0.103 
16.80 
0.109 
0.353 
20.28 
0.442 
9.081 

Mindlin[2] 
2.353 
22.66 
0.654 
41.704 
0.085 
31.17 
0.073 
0.228 
48.37 
0.342 
29.587 

Kirchhoff[1] 
0.666 
78.11 
0.807 
28.079 
0.030 
75.20 
0.042 
– 
– 
– 
– 

Zenkour[5] 
3.043 
– 
1.122 
– 
0.123 
– 
– 
0.442 
– 
0.486 
– 

10 
Present 
1.091 
5.390 
0.854 
0.0189 
0.052 
1.512 
0.0594 
0.445 
0.137 
0.345 
0.137 

Sayyad and Ghugal [15] 
1.095 
5.069 
0.843 
0.0312 
0.051 
3.591 
0.059 
0.460 
0.135 
0.346 
0.135 

Reddy[3] 
1.090 
5.537 
0.839 
0.0359 
0.048 
9.073 
0.059 
0.440 
0.141 
0.344 
0.141 

Mindlin[2] 
0.964 
16.43 
0.772 
0.1134 
0.044 
16.44 
0.051 
0.253 
0.343 
0.263 
0.343 

Kirchhoff[1] 
0.666 
42.28 
0.807 
0.0725 
0.030 
41.96 
0.042 
– 
– 
– 
– 

Zenkour[5] 
1.153 
– 
0.870 
– 
0.052 
– 
– 
0.627 
– 
0.400 
– 

Table 4. Nondimensional displacements and stresses for (0°/90°/90^{0}/0^{0}) symmetric LCP 


S 
Hypothesis 
% Error 
% Error 
% Error 
% Error 
% Error 


4 
Present 
3.213 
0.948 
3.66 
0.942 
1.70 
0.107 
21.7 
0.384 
2.755 
0.574 
1.10 


Sayyad and Ghugal [15] 
2.857 
0.915 
7.01 
0.919 
4.12 
0.101 
25.6 
0.379 
4.091 
0.507 
12.5 


Zenkour [05] 
– 
0.984 
– 
0.958 
– 
0.136 
– 
0.395 
– 
0.580 
– 


10 
Present 
1.196 
0.816 
1.17 
0.555 
0 
0.058 
2.66 
0.568 
2.249 
0.442 
10.5 


Sayyad and Ghugal [15] 
1.109 
0.813 
1.62 
0.546 
0.02 
0.055 
6.83 
0.517 
6.871 
0.411 
16.7 


Zenkour [05] 
– 
0.826 
– 
0.559 
– 
0.060 
– 
0.555 
– 
0.494 
– 


20 
Present 
0.818 
0.823 
0.182 
0.413 
0.00 
0.454 
0.21 
0.601 
2.101 
0.443 
4.76 


Sayyad and Ghugal [15] 
0.793 
0.819 
0.26 
0.411 
0.01 
0.044 
1.97 
0.563 
8.252 
0.379 
18.6 


Zenkour [05] 
– 
0.822 
– 
0.416 
– 
0.045 
– 
0.614 
– 
0.466 
– 


50 
Present 
0.702 
0.827 
0.546 
0.364 
0.10 
0.040 
0 
0.618 
3.342 
0.454 
2.46 


Sayyad and Ghugal [15] 
0.698 
0.823 
0.060 
0.363 
0.32 
0.040 
0.24 
0.5805 
9.244 
0.369 
20.8 


Zenkour [05] 
– 
0.823 
– 
0.364 
– 
0.040 
– 
0.6396 
– 
0.466 
– 


100 
Present 
0.6853 
0.828 
0.595 
0356 
0.00 
0.0398 
0 
0.6389 
1.035 
0.457 
2.38 


Sayyad and Ghugal [15] 
0.6842 
0.824 
0.109 
0.356 
0.00 
0.0398 
0 
0.5833 
9.643 
0.367 
21.5 


Zenkour [05] 
– 
0.823 
– 
0.356 
– 
0.0398 
– 
0.6456 
– 
0.4690 
– 

















From the results shown in Table 24, it has been noted that results of nondimensional 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 nondimensional 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 24 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 514 shows the variation of nondimensional displacements and stresses across the thickness of the LCP.
The nondimensional stresses are plotted along the xaxis and plate thickness (z/h) along the yaxis 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 nondimensional inplane normal stress ( ) for (0°/90°) antisymmetric LCP
Fig. 6 Through thickness variation of nondimensional transverse shear stress ( ) for (0°/90°) antisymmetric LCP
Fig. 7Through thickness variation of nondimensional inplane shear stress ( ) for (0°/90°) antisymmetric LCP
Fig. 8Through thickness variation of nondimensional inplane normal stress ( ) for (0°/90°/0^{0}) antisymmetric LCP
It may be easily inferred from Figs. 517, 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_{1}/E_{2}=25 (High elastic modulus along the fibre).
The Figs. 614 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. 516, the outofplane 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 inplane 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 nondimensional transverse displacement decreases in the sequence of stacking in order of antisymmetric (0^{0}/90^{0}), symmetric (0^{0}/90^{0}/0^{0}), symmetric (0^{0}/90^{0}/90^{0}/0^{0}) and antisymmetric (0^{0}/90^{0}/0^{0}/90^{0}) square LCPs.
Fig. 9Through thickness variation of nondimensional transverse shear stress ( ) for (0°/90°/0^{0}) symmetric LCP
Fig. 10Through thickness variation of nondimensional transverse shear stress ( ) for (0°/90°/0^{0}) antisymmetric LCP
Fig. 11Through thickness variation of nondimensional in plane normal stress ( ) for (0°/90°/90°/0°) symmetric LCP
Tables 35 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. 1516, validate nondimensional 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.1718. 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 PrepPost 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. 1925. The plate is discretized using Solid 185 finite elements.
The Hashin failure criteria are used in the 2dimensional approach to calculate stress at nodes. The contour plots for linear bending analysis of LCP analyzed in ANSYS environment is shown in Figs. 1923.
Fig. 12Through thickness variation of nondimensional in plane normal stress ( ) for (0°/90°/90°/0°) symmetric LCP
Fig. 13Through thickness variation of nondimensional transverse shear stress ( ) for (0°/90°/90°/0°) symmetric LCP
Fig. 14 Through thickness variation of nondimensional transverse shear stress ( ) for (0°/90°/90°/0°) symmetric LCP
Fig. 15 Variation of nondimensional central deflection ( ) of LCP for S=4.
Figure 19 shows transverse deflection for S=10 whereas by FEMMATLAB (present) and by Zenkour [5]. The percentages of error while compared with Zenkour [5] are 1.2422 % (Present ANSYS) and 1.2422 % (Present EEMATLAB). Fig.21 shows the nondimensional shear stresses computed by ANSYS as 0.1090 and by FEMANSYS 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 2023 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 nondimensional central deflection ( ) of LCP for S=10
Fig.17 The descretized (0°/90°) antisymmetric LCP using Solid 185 Element
Fig. 18. The SSBCs (0°/90°) antisymmetric LCP subjected to transverse load (UDL)
Fig. 19. The nondimensional central displacement of the (0°/90°) antisymmetric LCP using ANSYS
Fig. 20. The nondimensional inplane shear stress in the (0°/90°) antisymmetric LCP using ANSYS
Fig. 21. The nondimensional central displacement of the plate of (0°/90°/0^{0}) symmetric LCP using ANSYS
Fig. 22. The nondimensional inplane stress in the (0°/90°/0^{0}) symmetric LCP using ANSYS
Having recognized the reliability of the FE model built for the linear bending, the analysis of LCPs drawnout 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 nondimensional 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.
(20) 
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^{0})antisymmetric LCP with aspect ratio, S=20 results are shown in Fig. 26. The nondimensional transverse displacement, shown in Table 5, by FEMMATLAB (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 nondimensional central displacement of the (0°/90°/90^{0}/0^{0}) symmetric LCP using ANSYS
Table 5. Nondimensional transverse displacements ( ) for (0°/90°/0°/90°) for S=20. 

S 
Mesh size 
Nonlinear 
% Error 
20 
Present (4×4) 
0.7482 

Present (6×6) 
07579 


Present (8×8) 
0.7582 


Present (10×10) 
0.7582 
1.057 

FEM by Kant & Kommineni [26] 
0.7520 
1.866 

FEM by Ren & Hinton[08] 
0.7663 
 

In ANSYS environment 
0.780 (Fig.25). 
1.787 

40 
Present (4×4) 
0.7281 

Present (6×6) 
0.7335 


Present (8×8) 
0.7357 


Present (10×10) 
0.7357 
0.982 

FEM by Kant & Kommineni [26] 
0.7308 
1.644 

FEM by Ren & Hinton [08] 
0.7430 
 

In ANSYS environment 
0.7290 
1.844 
Fig. 24. The descretized (0°/90°/0^{0}/90^{0}) antisymmetric LCP using Solid 186
Fig. 25. The nondimensional transverse displacement of the (0°/90°/0^{0}/90^{0}) LCP using ANSYS (S=20)
Fig. 26 Loaddeflection diagram at specific load points in geometrically nonlinear analysis.
Table 6. Nondimensional transverse displacements ( ) for (0°/90°/0°/90°) for S=20 

S=20 

L (Present) 
NL (Present) 
NL by Kant & Kommineni [26] (C^{0} FEM) 
NL Zhang [07] (FEM) 

0 
0 
0 
0 
0 
50 
0.50253 
0.3255 
0.320 
0.3270 
100 
0.72218 
0.4886 
0.486 
0.4949 
150 
1.00680 
0.6010 
0.592 
0.6048 
200 
1.26520 
0.6911 
0.680 
0.6948 
250 
1.53130 
0.7582 
0.752 
0.7663 
The analysis of the results in Table 67, 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 nondimensional 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.
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 SSBCs. 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 twodimensional 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.
Table 7. Nondimensional transverse displacements ( ) for (0°/90°/0°/90°) for S=40 

S=40 

L (Present) 
NL (Present) 
NL Kant & Kommineni [26] (C^{0} FEM) 
NL Zhang [07] (FEM) 
0 
0 
0 
0 
0.44405 
0.2981 
0.293 
0.2913 
0.68024 
0.4656 
0.464 
0.4606 
0.97431 
0.5840 
0.582 
0.5771 
1.23800 
0.6651 
0.664 
0.6668 
1.50877 
0.7358 
0.738 
0.7430 
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:5188.
[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 HigherOrder Theory for Laminated Composite Plates. Journal of Applied Mechanics 1984; 51:745752.
[4] Pagano N. Exact Solutions for Rectangular Bidirecional Composites and Sandwich Plates. Journal of Composite materials 1970;4, 2034.
[5] Zenkour A. Threedimensional Elasticity Solution for Uniformly Loaded Crossply 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: 6979.
[7] Zhang Y, Kim K. Geometrically nonlinear analysis of laminated composite plates by two new displacementsbased 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: 217228.
[9] B Pandya, T Kant. Flexural analysis of laminated composites using refined higherorder C^{0} plate bending elements. Computer Methods in Appl. Mech. Eng. 1988; 66(2): 173198 .
[10] Savithri S, Varadan T. Large deflection analysis of laminated composite plates. HT. Journal of Nonlinear Mechanics 1993;28: 112.
[11] Qinm Q. Nonlinear analysis of thick plates by HT FE approach. Computers & Structures 1995; 61:227281.
[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:177190.
[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:194201.
[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:247267.
[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 123
[17] Sayyad A S, Ghugal Y. Flexure of crossply laminated plates using equivalent single layer trigonometric shear deformation theory. Structural Engineering & Mechanics 2014; 51(5):867891.
[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), 132.
[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 crossply composite laminates using Layerwise theory. Modares Mechanical Engineering 2015; 14(1):5666.
[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^{0} finite element geometrically nonlinear analysis of fibre reinforced composite and sandwich laminates based on a higherorder 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 firstorder shear deformation theory. Journal of Intelligent Material Systems and Structures 2019;30(4): 517535.
[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:354366.
[30] Mallek H, Jrad H, Wali M, Dammak F. Geometrically nonlinear analysis of FGCNTRC shell structures with surfacebonded 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 3Dshell structures based on the double directors shell theory with finite rotations.Steel and Composite Structures 2019;31(4): 397408.
[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, NguyenXuan, 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, NguyenXuan H. (2017). A novel threevariable shear deformation plate formulation: Theory and Isogeometric implementation. Computer Methods in Applied Mechanics and Engineering 2017;326:376–401.
[35] Nguyen NT, Hui D, Lee J,NguyenXuan H. An efficient computational approach for sizedependent analysis of functionally graded nanoplates. Computer Methods in Applied Mechanics and Engineering 2015;297:191–218.
[36] Nguyen H X, Nguyen T N, AbdelWahab, M, Bordas S P A, NguyenXuan, H, Vo T P. A refined quasi3D 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;6162: 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 higherorder 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 329334.
[40] Stephen Timoshenko, J M Gere.Theory of Elastic Stability. 1961; McGrawHill 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.