Open access peer-reviewed chapter

# Bending of Laminated Composite Plates in Layerwise Theory

Written By

Marina Rakočević

Submitted: January 27th, 2017 Reviewed: May 31st, 2017 Published: December 20th, 2017

DOI: 10.5772/intechopen.69975

From the Edited Volume

## Lamination

Edited by Charles A. Osheku

Chapter metrics overview

View Full Metrics

## Abstract

Determination of stress‐strain state in contemporary laminated composite plates containing layers with continuous unidirectional fibers requires the application of refined plate theories, which include layerwise theory. In contrast to homogeneous isotropic plates, heterogeneity of the anisotropic structure of laminated composite plates often leads to the appearance of imperfections in the connection between the layers. Mathematical models, which are formed on the assumption that the plate is homogeneous and isotropic, cannot properly include irregularities that can occur at the level of the layer in the process of manufacture, transportation, installation, or exploitation. Mathematical models of layerwise theory allow defining a more realistic stress‐strain state through the thickness of the plate, where consideration is carried out at the level of the layer. Additionally, this model makes possible to include delaminations that might occur on the connection between the individual layers. In this chapter, Reddy's layerwise theory is applied in order to determine equations for the problem of bending of laminated composite plates. The bending equations are solved by applying analytical method by means of double trigonometric series, as well as by using numerical methods based on the finite elements. This chapter presents examples for both applied approaches.

### Keywords

• laminated plate
• bending
• layerwise theory
• analytical solution
• finite element method

## 1. Introduction

Composite materials are obtained by a combination of two or more materials which together have characteristics that separately usually cannot be achieved with reasonable costs. During the last decades, a sudden increase of the use of modern composite materials is based on the achieved advanced properties which include their lower weight and higher strength and resistance compared to the classical and conventional composite materials. It is clear that the modern composite materials are materials of the future, because they allow the designers to choose the appropriate characteristics of the elements structure depending on the applicable problem. Due to its complexity, the researches of modern materials involve several scientific fields, and each of them, in their own way, contributes to their development.

Composite plates, made of layers of the various material and geometric characteristics, are called laminated plates. Ply or layer is the basic element of the laminated plate and it is made of fibers installed into the matrix. The fibers can be discontinuous, continuous unidirectional, bidirectional, woven unidirectional, or woven bidirectional. In the unidirectional laminated composites, the fibers are load‐bearing elements of the layer, while the matrix has a role to protect the fibers from external influences, to hold the fibers together and to perform uniform distribution of the influences to each of the fibers. The materials used for fibers have better material property and greater capacity compared to the matrix, and the geometrical characteristics of the fiber cross section are significantly smaller than its length. Materials used for fibers can be aluminum, copper, iron, nickel, steel, titanium, or organic materials such as glass, carbon, and graphite. A layer with unidirectional fibers has significantly better characteristics in fiber direction than in a direction perpendicular to the fiber.

Heterogeneity of anisotropic laminated composite plates often causes the appearance of a large number of imperfections that can occur at the level of the laminated plate or at the level of the layer, as well as at the local level of the fiber/matrix. General deformation of laminated plates is often defined by complex coupling between the axial deformation, bending, and shear deformation. In laminated composite plates for smaller aspect ratio, the importance of shear deformation is higher than in the corresponding homogeneous isotropic plates.

At the level of the layer, composites often contain concentrations of transverse stresses near the geometrical and material imperfections which lead to damage. Discontinuity between adjacent layers can occur in the stages of production, transportation, installation, or exploitation. This phenomenon in the literature is known as delamination. Delamination and its increase during exploitation is one of the most important parameters that affect the bearing capacity of the plate. Plates with delamination have reduced bearing capacity to bending. In the analysis of these plates, it is necessary to choose the mathematical models which can successfully incorporate delamination in the calculations.

Mathematical models for these particular problems need to determine the real stress‐strain state in the laminated plate, which requires the application of more accurate theories. In addition, it is important to find a balance between the desired accuracy and calculation costs.

The literature that studies the problems of laminated composite plates is significant and extensive. Available mathematical models for the laminated plates are based on assumptions of Equivalent Single Layer (ESL) theories and assumptions of LayerWise (LW) theories, see Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. ESL theories consider plate as one layer with equivalent stiffnesses. The most commonly applied ESL theories of laminated composite plates are Classical Laminated Plate Theory (CLPT), First-Order Shear Deformation Theory (FSDT), Third-Order Shear Deformation Theory (TSDT), and other high‐order shear theories. ESL theories are recommended to apply for thin and very thin laminated plates. More accurate ESL plate theories introduce shear deformation in the calculation which increases accuracy compared to the classical laminated theory. Layerwise theories consider laminated composite plate at the level of the layer. LW theories are recommended for moderately thick and thick laminated plates and for plates with a pronounced anisotropic behavior for which it is necessary to define a more realistic internal stresses distribution, as well as where it is necessary to include delamination and the other geometrical and material imperfections in calculation. On the basis of LW theories, many calculation procedures, analytical, and numerical methods were created in order to accurately analyze mentioned stresses of the laminated composite plates.

This chapter describes the application of layerwise theories in solving problems of bending of laminated plates, whereby the three‐dimensional (3D) problem of elasticity is reduced by solving differential equations with two‐dimensional (2D) variables in the plane plate. Additionally, the approximation of the displacement through the thickness of the plate is conducted with one‐dimensional interpolation functions of coordinate perpendicular to the plane of the layer. Mathematical model is based on the Generalized Laminated Plate Theory (GLPT) developed by Reddy and Reddy et al. [10, 11, 12, 13, 14].

Based on the adopted assumptions of the layerwise theory, bending equations are obtained by applying the principle of virtual displacements. In this chapter, the governing differential equations of LW theory are solved by applying analytical and numerical methods. Closed analytical solution of bending equations in LW theory is done by using the double trigonometric series of the simple supported rectangular plate with the angle of the layers orientation 0 and 90°, see Refs. [10, 18, 22, 23, 24, 25, 26].

Analytical solution for this mechanical problem is difficult to determine in case of complex plate geometry, arbitrary angles of layers orientation, arbitrary boundary conditions, as well as for various forms of nonlinearity. These are the reasons why the applied Finite Element Method (FEM) as one of the most used numerical methods was used, see Refs. [1, 10, 11, 12, 13, 14, 21]. The idea of the FEM is a representation of the domain as the sets of appropriate simple geometrical shapes called the finite elements. Discretization of the domain with finite elements makes this method a valuable tool in solving a large number of complex engineering problems.

Since the FEM is a method of discretization, errors can be classified into two categories: geometric discretization errors and interpolation errors. Accuracy of the results obtained by applying FEM depends on the choice of interpolation function in the plane (x, y) and interpolation function through the thickness of the plate. It also depends on the number of elements and the number of nodes of the adopted network in the plane plate and perpendicular to the plane plate. The increase of the number of finite element nodes gives rise to the increase of the number of degrees of freedom, making this calculation uneconomic. Regardless of the pointed disadvantages, FEM is much more general and more acceptable for the analysis of laminated composite plates than analytical methods. In this chapter, the variable kinematic finite element (with C0 continuity through the plate thickness) based on the assumptions of the partial LW theory is presented. The element is hierarchical and can include delamination. This finite element can be expanded in order to obtain a more complex finite element, by simply adding neglected elongation perpendicular to the plane of the plate.

## 2. Theoretical backgrounds

Laminated composite plate, loaded by transversal load, is considered in Cartesian coordinate system (x, y, z). The plate includes N orthotropic layers wherein each layer carries a predetermined direction in the plane of the plate (x, y). Layer orientation angle is defined by angle θ, angle between continuous fibers, and axis x of the global coordinate system of the plate.

It is assumed that the dilation εz perpendicular to the middle plane of the plate is neglected which results in a constant deflection through the plate thickness w(x, y, z) = w(x, y). A displacement field is determined by components:

u(x,y,z)=u(x,y)+U(x,y,z)ν(x,y,z)=ν(x,y)+V(x,y,z)w(x,y,z)=w(x,y)E1

Displacement components u(x, y, z) and v(x, y, z) are obtained as a sum of displacements of the middle plane of the plate u(x, y), v(x, y) and the additional displacements U(x, y, z) and V(x, y, z) through the plate thickness. Reduction of a 3D problem to a 2D problem is carried out by introducing the following relations:

U(x,y,z)=J=1nuJ(x,y)ψJ(z)V(x,y,z)=J=1nvJ(x,y)ψJ(z)E2

where uJ and vJ are the displacement components of Jth node through the plate thickness in the directions of x‐ and y‐axes, respectively; ψJ(z) is the one‐dimensional (1D) interpolation function.

The governing equations of the laminated composite plates in case of bending are solved by displacement method of analysis, so that the primary variables are displacement components defined by Eqs. (1) and (2). The primary variables can be conveniently displayed in a vector form:

{Δ}={uvw}{ΔJ}={uJvJ}E3

where Δ are the components of the displacement vector of the middle plane plate and ΔJ are the components of the displacement vector of the Jth node through the plate thickness.

By relations (Eq. (2)), the discretization through the plate thickness is performed, wherein the uJ and vJ are unknown nodal displacements and ψJ(z) are adopted general 1D interpolation functions that can be polynomials. The degree of a polynomial is related to the number of nodes in which determined unknown nodal displacements occur, where

n=pN+1E4

where p is the polynomial degree; n is the number of nodes through the plate thickness; N is the number of layers.

For linear interpolation p = 1, the number of nodes through the plate thickness is one larger than the number of layers n = N+1, and the linear interpolation functions are determined by

ϕJ(z)={ψ1J1(z)=zzJ1zJzJ1,zJ1<z<zJψ2J(z)=zJ+1zzJzJ1,zJ<z<zJ+1}E5

where zJ is the coordinate of the Jth node which is located between the layers of J and J‐1.

Graphical display of the laminated composite plates with linear interpolation through the plate thickness is shown in Figure 1. By analyzing separate layers, this LW theory can be understood as a first‐order shear deformation theory for each layer of the observed laminated plate. A layer of the layered plates has an effect on the previous and the next layer of the laminate.

Based on Eqs. (1) and (2), strain components can be expressed as follows:

{ε}={u,xv,yu,y+v,xw,xw,y}{εJ}={U,xJV,yJU,yJ+V,xJUJVJ}E6

where the members of the vector ε determine the deformation of the middle plane of the plate, and the members of the vector εJ determine the deformation through the thickness of the plate. The short form of Eq. (3) is defined by

{ε}=[H]{Δ},{εJ}=[H¯]{ΔJ}E7

where members of matrix H and H¯ are the first derivative of the interpolation functions with respect to the global coordinates x and y.

For an orthotropic material, the constitutive equations of the kth layer are determined by

{σ}(k)=[ Q¯11Q¯12Q¯1600Q¯12Q¯22Q¯2600Q¯16Q¯26Q¯6600000Q¯44Q¯45000Q¯45Q¯55 ](k){ε}E8

In the previous equation, Q¯ij are members of the transformed plane stress‐reduced stiffnesses of kth layer which are obtained as follows; see Ref. [1]:

Q¯11=Q11cos4θ+2(Q12+2Q66)sin2θcos2θ+Q22sin4θQ¯12=(Q11+Q224Q66)sin2θcos2θ+Q12(sin4θ+cos4θ)Q¯22=Q11sin4θ+2(Q12+2Q66)sin2θcos2θ+Q22cos4θQ¯16=(Q11Q122Q66)cos3θsinθ+(Q12Q22+2Q66)sin3θcosθQ¯26=(Q11Q122Q66)sin3θcosθ+(Q12Q22+2Q66)cos3θsinθQ¯66=(Q11+Q222Q122Q66)sin2θcos2θ+Q66(sin4θ+cos4θ)Q¯44=Q44cos2θ+Q55sin2θQ¯45=(Q55Q44)cosθsinθQ¯55=Q55cos2θ+Q44sin2θE9

where θ is the angle between the fibers and the x‐axis of the global coordinate system of the plate and Qij is the plane stress‐reduced stiffnesses of kth layer (including materials constants).

Transverse stresses σxzk and σyzk obtained using Eq. (8) have a constant value along the kth layer. In order to determine a more realistic distribution of transverse shear stresses through the thickness approximate method that is given in Refs. [23, 26] can be used.

The internal forces in a cross section of the laminated plate are defined as integrals of stresses

(Nx,Ny,Nxy,Qx,Qy)=0h(σxx,σyy,σxy,σxz,σyz)dz(NxJ,NyJ,NxyJ)=oh(σxx,σyy,σxy)ϕJ(z)dz(QxJ,QyJ)=0h(σxz,σyz)ϕ,zJdzE10

where σx, σy, σxy, σxz, σyz are stress components.

The internal forces can be expressed in terms of the displacements by inserting Eq. (6) into Eq. (7) and then in the integrals in Eq. (10):

{N}=[A]{ε}+J=1n[BJ]{εJ}{NJ}=[BJ]{ε}+J=1n[DIJ]{εJ}E11

where [A], [BJ], and[DIJ] are stiffnesses of a laminated plate

Aij=J=1Nzkzk+1Q¯ij(k)dz,(i,j=1,2,6,4,5)E12
BijJ=k=1Nzkzk+1Q¯ij(k)ψJdz(i,j=1,2,6)BijJ=k=1Nzkzk+1Q¯ij(k)ψ,zJdz(i,j=4,5)E13
DijJI=k=Nzkzk+1Q¯ij(k)ψJψIdz(i,j=1,2,6)DijJI=k=1Nzkzk+1Q¯ij(k)ψ,zJψ,zIdz(i,j=4,5)E14

For the adopted linear interpolation functions (Eq. (5)) from Eqs. (12)(14) can be obtained

Aij=k=1NQ¯ijK(zk1zk),(i,j=1,2,6,4,5)BijJ=Q¯ijJ1(zJzJ1)2+Q¯ijJ(zJ+1zJ)2(i,j=1,2,6)BijJ=Q¯ijJ1Q¯ijJ(i,j=4,5)forJ=IDijIJ=Q¯ijJ1(zJzJ1)3+Q¯ijJ(zJ+1zJ)3(i,j=1,2,6)DijIJ=Q¯ijJ1(zJzJ1)+Q¯ijJ(zJ+1zJ)(i,j=4,5)forJ=I±1DijIJ=DijJI=Q¯ijJ(zJ+1zJ)6(i,j=1,2,6)DijIJ=DijJI=Q¯ijJ(zJ+1zJ)(i,j=4,5)E15

It is noted that the elements of matrix A are determined for the laminated plate as a whole and do not depend on the adopted discretization through the thickness of the plate.

The governing equations of motion can be derived by using the principle of displacements. For zero values of stiffnesses A16 = A26 = A45 = B16J = B26J = B45J = D16JI = D26JI = D45JI, the equations take the form:

A11u,xx+A12v,yx+A66(u,yy+v,xy)+J=1n[B11JU,xxJ+B12JV,yxJ+B66J(U,yyJ+V,xyJ)]=0A12u,xy+A22v,yy+A66(u,yx+v,xx)+J=1n[B12JU,xyJ+B22JV,yyJ+B66J(U,yxJ+V,xxJ)]=0A55w,xx+A44w,yy+J=1n[B55JU,xJ+B44JV,yJ]+q=0B11Ju,xx+B12Jv,xy+B66J(u,yy+v,xy)B55Jw,x+J=1n[D11JIU,xxJ+D12JIV,yxJ+D66JI(U,yyJ+V,xyJ)D55JIUI]=0B11Ju,xy+B21Jv,yy+B66J(u,yx+v,xx)B44Jw,y+J=1n[D12JIU,xyJ+D22JIV,yyJ+D66JI(U,yxJ+V,xxJ)D44JIVI]=0E16

where q is load perpendicular to the middle plane of the plate.

The system of Eq. (16) contains (3+2n) differential equations with the unknown displacements of the middle plane of the plate (u, v, w), and the displacements through the plate thickness UJ and VJ, where J=1, …, n. These equations can be solved by using analytical or numerical and approximate methods.

This chapter shows the analytical solution that is obtained by using the double trigonometric series for simply supported plate, as well as the numerical solution obtained by using the finite element method as one of the currently most applied numerical methods.

## 3. Analytical solution

### 3.1. Navier’s solution in layerwise theory

Simply supported rectangular laminated plate a x b with N orthotropic layers with angles of orientation θ=0° or 90° is considered. Displacements (Eqs. (1) and (2)) are given in the form of double trigonometric series

u=m,nXmncosmπxasinnπybv=m,nYmnsinmπxacosnπybw=m,nWmnsinmπxasinnπybUJ=m,nRmnJcosmπxasinnπybVJ=m,nSmnJsinmπxacosnπybE17

where α=mπa;β=nπbJ=1,,n

Boundary conditions for simply supported plate are met:

u=w=VJ=Nx=NxJ=0x=0,a J=1,nv=w=UJ=Ny=NyJ=0y=0,b J=1,nE18

Transverse load is displayed in the form of double Fourier series:

q(x,y)=m,nQmnsinmπxasinnπybE19

where Qmn is the Furie’s coefficients which depend on the type of transverse load.

When derivatives of the displacement (Eq. (17)) are substituted into Eq. (16), for each pair (m,n), a system of equations is obtained:

[[K][KJ][KJ]T[CJI]]{XmnYmnWmnRmnJSmnJ}={00Qmn00}E20

Elements of the submatrices K, KJ, and KJI are given in the function of the stiffnesses of laminated plate (Eq. (15)) and variables α and β [26].

By solving Eq. (20), the unknown coefficients Xmn, Ymn, Wmn, Rmn, and Smn for each pair (m, n) are determined. Additionally, the application of Eq. (17) and the summing of the unknown displacements are determined.

The strains are determined using Eq. (6) and then the stresses using Eq. (8). Stresses in the plane of the plate are determined by the expressions:

σxx(x,y,z)=m,n{ [Q¯11α(Xmn+J=1nRmnJϕJ(z))+Q¯12β(Ymn+J=1nSmnJϕJ(z)) ]sinαx sinβy}σyy(x,y,z)=m,n{ [ Q¯12α(Xmn+J=1nRmnJϕJ(z))+Q¯22β(Ymn+J=1nSmnJϕJ(z)) ]sinαx sinβy }σ(x,y,z)=Q¯66m,n{ [ β(Xmn+J=1nRmnJϕJ(z))+α(Ymn+J=1SmnJϕJ(z)) ]cosαx sinβy}E21

Shear stresses σxz and σy are constant through the layer J and can be obtained with

σxzj(x,y,z)=m,n{Q¯45j{1hk[(Smnj+1Smnj)+βWmn]sinαx cosβy}+Q¯55j{1hk[(Rmnj+1Rmnj)+αWmn]cosαx sinβy}}σxzj(x,y,z)=m,n{Q¯44j{1hk[(Smnj+1Smnj)+βWmn]sinαx cosβy}+Q¯45j{1hk[(Rmnj+1Rmnj)+αWmn]cosαx sinβy}}E22

A more realistic change of shear stresses σxz and σy through the plate thickness can be obtained using approximate procedures. In this chapter, approximate procedure as described in Ref. [23] was used.

### 3.2. Numerical examples

For a simply supported rectangular plate, stresses are determined in dimensionless form:

σ¯xx=1qs2σxx,σ¯yy=1qs2σyy,σ¯xz=1qsσxz,σ¯yz=1qsσyzE23

where s=a/h, a=b are plate dimensions in the plane of the plate; h is the plate thickness

All layers of the laminated plate have the same thickness and material properties. For predefined material characteristics of layers E1, E2, G12, G21, G23, ν12, and ν13, changes of dimensionless stresses σxx, σyy, σyy, and σxz are displayed in the points with coordinates A=1,105662(a/2) and B=1,894338(a/2)

σ¯xx=σ¯xx(A,A,z),σ¯yy=σ¯yy(A,A,z),σ¯xz=σ¯xz(B,A,z),σ¯yz=σ¯yz(A,B,z)E24

Example 1: The square laminated composite plate with five layers 0°/90°/0°/90°/0° and the aspect ratio a/h=4 loaded by a bi‐sinusoidal is considered. Characteristics of the layers are E1/E2=25, E2=1, G12=G13=0.5, G23=0.2, ν12=ν13=0.25. The change of dimensionless normal stresses σxx, σyy and shear stresses σyz, σxz are displayed in Figures 2 and 3.

Example 2: The square laminated composite plate with six layers 0o/90o/0o/90o/0o/90o and the aspect ratio a/h=10 loaded by uniformly distributed load of intensity q is considered. The material properties of the layers are as in Example 1. Figures 4 and 5 graphically show the change of dimensionless stresses σxx, σyy, σyz, and σxz through the thickness of the plate.

## 4. Finite element method

### 4.1. Stiffness matrix of the kinematic variable finite element

The components of displacement of finite element can be written as a combination of two‐dimensional interpolation function ϕi and nodal displacements ui, vi, wi, ui j, vi j

(u,v,w,uj,vj)=i=1m(ui,vi,wi,uij,vij)ϕiE25

Interpolation functions ϕi are defined in the natural coordinates (ξ,η). The isoparametric formulation of the finite element is adopted, so that the geometry is interpolated the same way as the primary variables (displacements, see Eq. (25)):

x=i=1mxiϕi(ξ,η)y=i=1myiϕi(ξ,η)E26

where x, y are point coordinates and xi, yi are coordinates of the finite element.

According to Eqs. (6) and (25), matrices H and H¯ are defined as follows:

[H](5x3m)=[ϕ1x00ϕ2x00ϕmx000ϕ1y00ϕ2y00ϕmy0ϕ1yϕ1x0ϕ2yϕ2x0ϕmyϕmx000ϕ1x00ϕ2x00ϕmx00ϕ1y00ϕ2y00ϕmy]E27
[H¯](5x2m)=[ϕ1x0ϕ2x0ϕmx00ϕ1y0ϕ2y0ϕmyϕ1yϕ1xϕ2yϕ2xϕmyϕmxϕ10ϕ20ϕm00ϕ10ϕ20ϕm]E28

where m is the number of nodes in the plane of the plate.

The first derivatives of the ϕi with respect to the global coordinates x, y can be obtained by using the local coordinates ξ, η as follows:

{ϕixϕiy}=[x,ξy,ξx,ηy,η]1{ϕ,ξiϕ,ηi}=[J]1{ϕ,ξiϕ,ηi}E29

where J is the Jacobian matrix of the transformation in the natural coordinates:

[J]=[xξyξxηyη]E30

The second derivatives of the interpolation functions, which are used during determinations of the shear stresses perpendicular to the plane plate, are

{ϕ,xxϕ,yyϕ,xy}=[J1]1({ϕi,ξξϕi,ηηϕi,ξη}[J2]{ϕi,xϕi,y})E31

where

[J1]=[x,ξ2y,ξ22y,ξ2x,ξ2x,η2y,η22y,η2x,η2x,ηξy,ηξx,ηξ+y,ηξ]E32
[J2]=[x,ξξy,ξξx,ηηy,ηηx,ηξy,ηξ]E33

According to the principle of stationary potential energy, the first variation of the potential energy must be zero for equilibrium conditions. The first variation of the potential energy of the finite element is

δΠe=VeεTQ¯δεdVE34

By using Eq. (7), the strains can be written in the following form:

{ε}=[H]{Δ}+[H¯]j=1nϕJ{ΔJ}E35

where ϕ J are the interpolation functions perpendicular to the plane of the plate or its derivatives:

ϕJ(z)=[ψj(z)ψj(z),z]E36

Inserting Eq. (35) into Eq. (34), the first variation of the potential energy of the finite element becomes

δΠe=Ve(ΔTHTH¯TJ=1nϕJΔJT)Q¯(HδΔ+H¯I=1nϕIδΔI)dVE37

After multiplication, the stiffness matrix of the finite element takes the form:

ke=Ve(HTQ¯H+J=1nHTQ¯ϕJH¯+I=1nH¯TQ¯ϕIH+I=1nJ=1nH¯TQ¯ϕIϕJ H¯)dV=Ae(HTAH+J=1nHTBJH¯+I=1nH¯TBIH+I=1nJ=1nH¯TDIJH¯)dAE38

The stiffness matrix of the finite element can be written in a matrix from as

ke=[k11k112k212kn12k121k1122k1222k1n22k221k1222k2222k2n22kn21kn122kn222knn22]E39

From the structure of the matrix (Eq. (39)), four submatrices can be noticed, where elements (see Eq. (38)) are determined as follows:

k11=AeHTAHdAi=1,,nki12=AeHTBiH¯dAi=1,,nki21=AeH¯TBiHdAi=1,,nkij22=AeH¯TDijH¯dAi,j=1,,nE40

The area integral can be written in the form of the natural coordinates:

dA=dxdy=detJdξdηE41

as follows:

Ae()dA==1+11+1()detJdξdηE42

where det J is the determinant of the Jacobian matrix of the transformation according to Eqs. (26) and (30):

|J|=|i=1mxiϕξi=1myiϕξi=1mxiϕηi=1myiϕη|E43

Thereafter, the submatrices of the stiffness matrix (Eq. (39)) are given by

k11=1+11+1HTAHdetJdξdηi=1,,nki12=1+11+1HTBiH¯detJdξdηi=1,,nki21=1+11+1H¯TBiTHdetJdξdηi=1,,nkij22=1+11+1H¯TDijH¯detJdξdηi,j=1,,nE44

### 4.2. Interpolation functions and numerical integration

Interpolation or shape functions in finite element analysis depend on the dimensionality of the problem and type of elements used for discretization of the problem. Layerwise theory introduces two types of interpolation functions. The first type is the interpolation functions in the plane (x,y) that can be selected from the family of well‐known two‐dimensional shape functions. The degree and type of interpolation functions depend on the type of finite element and the number of nodes in which variables are defined. The second type is the interpolation functions through the plate thickness. These interpolation functions are the functions of one variable (z) and can be linear, quadratic, or higher order functions.

For the models of rectangular finite element, it is possible to implement a wide range of standard 2D interpolation functions in the plane of the plate, as well as 1D interpolation functions through the thickness of the layer. At our disposal are 2D finite elements such as E4—four‐node Lagrange rectangular element, E8—eight‐node Serendipity rectangular element, E9—nine‐node Lagrange rectangular element, E12—12‐node Serendipity rectangular element, and E16—16‐node Lagrange rectangular element. Each of these elements can be combined with one or more 1D Lagrange elements: L—linear element, Q—square element, C—cubic element, in order to create a wide spectrum of a variety of finite elements based on the layerwise theory. We use notation E16‐L4 to denote 2D element with 16 nodes combined with four linear 1D elements through the thickness of the plate.

For laminated plate containing N layers with adopted linear interpolation through the thickness (number of nodes is n=N + 1) and the interpolation with m nodes in the plane of the plate, the total number of degrees of freedom of one finite element is (3 + 2n)m. For the finite element displayed in Figure 6, m=8 and N=3, which leads to the number of degrees of freedom (3+2 × 4) × 8=88. It is clear that the number of degrees of freedom of the finite element is rapidly growing with the increasing degree of interpolation functions.

To determine the stiffness matrix of elements, it is necessary to perform numerical integration. If Gauss‐Legendre integration is used, then the integral is calculated by

I(ξ,η)=1+11+1f(ξ,η)dξdη=i=1mj=1mωiωjf(ξ¯,η¯)E45

where ωi and ωj are the weights of the integration, and (ξ¯,η¯) coordinates of the points at which the integration is performed. The weights and the position of integration points can be found in the wide literature for the finite element method.

### 4.3. Governing equations

Governing equations in the finite element method are calculated from the principle of virtual displacements. The first variation of the potential of the system of finite elements is

Π=δUTKUδUTFE46

By applying the principle of virtual displacements, the system of algebraic equations of the finite element method is obtained:

KU=FE47

where U is the vector of the nodal displacements, F is the vector of the nodal forces, and K is the stiffness matrix, all for system of assembly of finite elements formed in the global coordinate system.

Methods of forming the stiffness matrix K and the vectors U and F are generally known and can be found in the expanded literature that deal with the finite element method. Solutions of the system of algebraic equations (Eq. (46)) are the nodal displacements.

### 4.4. Numerical examples

A square four‐layer simply supported laminated plate with the arrangement of layers 60°/‐45°/‐45°/60° with boundary conditions u = v = w = 0 for x = 0, a and y = 0, b, loaded by the uniformly distributed load intensity q, was considered. The adopted finite element is E4–L1 with layers of the same thickness and the same material properties. Due to the demanding and very extensive process of calculation with a large number of unknowns, a relatively rough mesh of finite elements with 4×4 finite elements is adopted.

Material characteristics of the layers are E1=25, E2=1, G12=G13=0.5, G23=0.2, ν1213=0.25.

Stresses are determined in a dimensionless form according to Eq. (23) described in Section 2 of this chapter, wherein the shear stress in the plane of the plate is also given in the same dimensionless form as the normal stresses. Stresses are calculated in points with coordinates:

σ¯xx(0.4471688a;0.4471688a)σ¯yy(0.4471688a;0.4471688a)σ¯xy(0.05283122a;0.05283122a)σ¯yz(0.4471688a;0.9471688a)σ¯xz(0.9471688a;0.4471688a)

Changing of dimensionless stresses through the thickness of the laminated plate is presented in Figures 79.

## 5. Conclusion

This chapter refers to the application of Reddy's layerwise theory, based on the assumption that dilatation perpendicular to the middle plane of the laminated plate is neglected. In the bending equations for laminated plates, the primary variables of node J are displacements through the thickness of the plate uJ(x,y,z) and vJ(x,y,z) and deflections wJ(x,y,z)=w(x,y).

Analytical solution of the bending equations in the layerwise theory is presented for simply supported rectangular laminated plate which contains orthogonal layers. A graphical representation of the distribution of dimensionless stresses σxx, σyy, σyz, and σxz through the thickness for the square plate with five layers 0°/90°/0°/90°/0° and the aspect ratio a/h=4, which is loaded by bi‐sinusoidal load, as well as for six‐layer plate 0°/90°/0°/90°/0°/90° with the aspect ratio a/h=10, which is loaded with uniformly distributed load, is presented. Obtained results indicate that the analytical solution provides determination of more realistic change of interlayer stresses, and it also allows the incorporation of delamination in the calculation procedures. However, the analytical solution is difficult to define for complex geometry of the plates, for arbitrary angle of orientation layers, for arbitrary boundary conditions, as well as for various forms of nonlinearity. These are the reasons why in layerwise theory the use of numerical method such as finite element method is usually recommended.

In Chapter 4, a graphical representation of the distribution of dimensionless stresses σxx, σyy, σxy, σyz, and σxz through the thickness of the square, simply supported rectangular laminated plate with four layers 60°/‐45°/‐45°/60° and the aspect ratio a/h=10, which is loaded by uniformly distributed load, is presented. The accuracy of the solution, obtained by using the finite element method, is highly influenced by the selection of interpolation functions in the plane (2D), by the selection of interpolation functions through the thickness (1D) of the plate. Additionally, as in general for this calculation method, a number of finite elements of the analyzed domain highly influence the accuracy of obtained solution. However, arbitrary increase of the number of finite element nodes in the plane of the plate and through the thickness of the plate suddenly drastically increases the number of degrees of freedom, which gives rise to the uneconomical nature of the finite element method for more complex mechanical and engineering problems. Regardless of this deficiency, the finite element method is a more general and acceptable for the analysis of laminated composite plates than the proposed analytical method.

Aside the drawbacks and disadvantages of the analytical solution, the same should not be excluded from the engineering and scientific practice, since it represents the ultimate convergence and accuracy criteria for any other numerical or approximate methods, including the finite element method presented in this chapter.

## References

1. 1. Reddy JN. Mechanics of Laminated Composite Plate: Theory and Analysis. New York, NY: CRC Press; 1997. ISBN 0‐8493‐3101‐3
2. 2. Whitney JM. The effect of transverse shear deformation in the bending of laminates plates. Journal of Composite Materials. 2004;3:534-547
3. 3. Reissner E. A consistent treatment of transverse shear deformations in laminated anisotropic plates. The American Institute of Aeronautics and Astronautics Journal. 1972;10(5):716‐718
4. 4. Reddy JN. Energy and Variational Methods in Applied Mechanics. New York, NY: John Wiley; 1984
5. 5. Christensen RM, Lo KH, Wu EM. A high‐order theory of plate deformation, part1: Homogeneous plates. Journal of Applied Mechanics. 1977;44(7):636‐668
6. 6. Christensen RM, Lo KH, Wu EML. A high‐order theory of plate deformation, part 2: Laminated plates. Journal of Applied Mechanics. 1977;44(4):669‐676
7. 7. Vuksanović DJ. Linear analysis of laminated composite plates using single layer higher‐order discrete models. Composite Structures. 2000;48:205‐211
8. 8. Reddy JN. A simple higher‐order theory for laminated composite plates. Journal of Applied Mechanics. 1984;51:745‐752
9. 9. Reddy JN. A refined nonlinear theory of plates with transverse shear deformation. International Journal of Solids and Structure. 1984;20(9):881‐906
10. 10. Reddy JN, Barbero EJ, Teply JL. A plate bending element based on a generalized laminated plate theory. International Journal for Numerical Methods in Engineering. 1989;28:2275‐2292
11. 11. Reddy JN, Barbero EJ, Teply JL. An accurate determination of stresses in thick laminates using a generalized plate theory. International Journal for Numerical Methods in Engineering. 1990;29:1‐14
12. 12. Redyy JN, Barbero EJ. Modelling of thick composites using a layerwise laminate theory. International Journal for Numerical Methods in Engineering. 1993;36:655‐677
13. 13. Reddy JN, Robbins DH. Theories and computational models for composite laminates. Applied Mechanics Review. 1994;47(6):147‐165
14. 14. Reddy JN. Theory and analysis of laminated composite plates. Mechanics of Composite Materials and Structures. 1999;361:1‐79
15. 15. Ferreira AJM, Fasshauer GE, Batra RC, Rodrigues JD. Static deformations and vibration analysis of composite and sandwich plates using a layerwise theory and RBF‐PS discretizations with optimal shape parameter. Composite Structures. 2008;86:328‐343
16. 16. Vuksanović Ð, Rakočević M. Kompozitne ploče sa delaminacijama‐analitičko rješenje za slobodno oslonjene ploče. In: Simpozijum o istraživanjima i primjeni savremenih dostignuća u našem gradevinarstvu u oblasti materijala i konstrukcija (JUDIMK); Oktobar 2002. pp. 131‐136
17. 17. Vuksanović Ð, Rakočević M. Opšta teorija laminarnih ploča–analitičko rješenje za slobodno oslonjene ploće, In: Simpozijum o istraživanjima i primjeni savremenih dostignuća u našem graevinarstvu u oblasti materijala i konstrukcija, JUDIMK; Oktobar 2002. pp. 125‐130
18. 18. Vuksanović Ð, Ćetković M. Analytical solution for multilayer plates using general layerwise plate theory. Facta Universitates Series. Architecture and Civil Engineering. 2005;3(2):121‐136
19. 19. Cetković M, Vuksanović D. Bending, free vibrations and buckling of laminated composite and sandwich plates using a layerwise displacement model. Composite Structures. 2009;88(2):219‐227
20. 20. Cetković M, Vuksanović D. Large deflection analysis of laminated composite plates using layerwise displacement model. Structural Engineering & Mechanics. 2011;40(2):257‐277
21. 21. Rakočević M. Stress in multi‐layered composite plates. Gradevinar. 2005;57(7):503‐509
22. 22. Rakočević M. Analysis of laminated composite plate. Gradevinar. 2011;63(9/10):819‐825
23. 23. Rakocevic M. Approximate procedure for calculation of shear stresses σxz and σyz. Journal of Applied Engineering Science. 2012;10(1):37‐42
24. 24. Rakočević M. Analysis of simply supported laminated composite plates using double trigonometric series. In: Proceedings of the 16th International Symposium of Macedonian Association of Structural Engineers (MASE 2016); October 2015; Ohrid, Macedonia
25. 25. Rakocevic M. Analytical solution for simply supported laminated composite plate based on partial layerwise theory. Journal of Applied Engineering Science. 2016;14(1):102‐108
26. 26. Rakočević M, Popović S, Ivanišević N. A computational method for laminated composite plates based on layerwise theory. Composites Part B. 2017;122:202‐218

Written By

Marina Rakočević

Submitted: January 27th, 2017 Reviewed: May 31st, 2017 Published: December 20th, 2017