Bending of Laminated Composite Plates in Layerwise Theory

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.

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

$\begin{array}{l}u{\displaystyle \left(x,y,z\right)}=u{\displaystyle \left(x,y\right)}+U{\displaystyle \left(x,y,z\right)}\hfill \\ \nu {\displaystyle \left(x,y,z\right)}=\nu {\displaystyle \left(x,y\right)}+V{\displaystyle \left(x,y,z\right)}\hfill \\ w{\displaystyle \left(x,y,z\right)}=w{\displaystyle \left(x,y\right)}\hfill \end{array}$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:

$\begin{array}{c}U{\displaystyle \left(x,y,z\right)}={\displaystyle \sum _{J=1}^n{u}^J{\displaystyle \left(x,y\right)}{\psi }^J{\displaystyle \left(z\right)}}\\ V{\displaystyle \left(x,y,z\right)}={\displaystyle \sum _{J=1}^n{v}^J{\displaystyle \left(x,y\right)}{\psi }^J{\displaystyle \left(z\right)}}\end{array}$E2

where u^J and v^J 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:

${\displaystyle \left\{\Delta \right\}}={\displaystyle \left\{\begin{array}{c}u\\ v\\ w\end{array}\right\}}\hspace{1em}{\displaystyle \left\{{\Delta }^J\right\}}={\displaystyle \left\{\begin{array}{c}{u}^J\\ v^J\end{array}\right\}}$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 u^J and v^J 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+1$E4

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

${\varphi }^J{\displaystyle \left(z\right)}={\displaystyle \left\{\begin{array}{c}{\psi }_1^{J-1}{\displaystyle \left(z\right)}=\frac{z-z_{J-1}}{z_J-z_{J-1}},z_{J-1}<z<z_J\\ {\psi }_2^J{\displaystyle \left(z\right)}=\frac{z_{J+1}-z}{z_J-z_{J-1}},z_J<z<z_{J+1}\end{array}\right\}}$E5

where z_J 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.

Figure 1.

3D geometry of the laminated composite plates and discretization through the plate thickness with linear interpolation.

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

${\displaystyle \left\{\epsilon \right\}}={\displaystyle \left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{u}_{,x}\\ v_{,y}\\ u_{,y}+v_{,x}\end{array}\\ w_{,x}\end{array}\\ w_{,y}\end{array}\right\}}\hspace{1em}{\displaystyle \left\{{\epsilon }^J\right\}}={\displaystyle \left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{U}_{,x}^J\\ V_{,y}^J\\ U_{,y}^J+V_{,x}^J\end{array}\\ U^J\end{array}\\ V^J\end{array}\right\}}$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

$\left\{\epsilon \right\}={\displaystyle \left[H\right]}\left\{\Delta \right\},\hspace{0.5em}\left\{{\epsilon }^J\right\}={\displaystyle \left[\overline{H}\right]}\left\{{\Delta }^J\right\}$E7

where members of matrix H and $\overline{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

${\left\{\sigma \right\}}^{(k)}={\left[\begin{array}{ccccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& {\overline{Q}}_{16}& 0& 0\\ {\overline{Q}}_{12}& {\overline{Q}}_{22}& {\overline{Q}}_{26}& 0& 0\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& {\overline{Q}}_{66}& 0& 0\\ 0& 0& 0& {\overline{Q}}_{44}& {\overline{Q}}_{45}\\ 0& 0& 0& {\overline{Q}}_{45}& {\overline{Q}}_{55}\end{array}\right]}^{(k)}\left\{\epsilon \right\}$E8

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

$\begin{array}{l}\begin{array}{l}{\overline{Q}}_{11}=Q_{11}{\text{cos}}^4\theta +2(Q_{12}+2Q_{66}){\text{sin}}^2\theta {\text{cos}}^2\theta +Q_{22}{\text{sin}}^4\theta \hfill \\ {\overline{Q}}_{12}=(Q_{11}+Q_{22}-4Q_{66}){\text{sin}}^2\theta {\text{cos}}^2\theta +Q_{12}{\displaystyle \left({\text{sin}}^4\theta +{\text{cos}}^4\theta \right)}\hfill \end{array}\hfill \\ \begin{array}{l}{\overline{Q}}_{22}=Q_{11}{\text{sin}}^4\theta +2(Q_{12}+2Q_{66}){\text{sin}}^2\theta {\text{cos}}^2\theta +Q_{22}{\text{cos}}^4\theta \hfill \\ {\overline{Q}}_{16}=(Q_{11}-Q_{12}-2Q_{66}){\text{cos}}^3\theta \text{sin}\theta +(Q_{12}-Q_{22}+2Q_{66}){\text{sin}}^3\theta \text{cos}\theta \hfill \\ {\overline{Q}}_{26}=(Q_{11}-Q_{12}-2Q_{66}){\text{sin}}^3\theta \text{cos}\theta +(Q_{12}-Q_{22}+2Q_{66}){\text{cos}}^3\theta \text{sin}\theta \hfill \end{array}\hfill \\ \begin{array}{l}{\overline{Q}}_{66}=(Q_{11}+Q_{22}-2Q_{12}-2Q_{66}){\text{sin}}^2\theta {\text{cos}}^2\theta +Q_{66}{\displaystyle \left({\text{sin}}^4\theta +{\text{cos}}^4\theta \right)}\hfill \\ {\overline{Q}}_{44}=Q_{44}{\text{cos}}^2\theta +Q_{55}{\text{sin}}^2\theta \hfill \\ \begin{array}{l}{\overline{Q}}_{45}=(Q_{55}-Q_{44})\text{cos}\theta \text{sin}\theta \hfill \\ {\overline{Q}}_{55}=Q_{55}{\text{cos}}^2\theta +Q_{44}{\text{sin}}^2\theta \hfill \end{array}\hfill \end{array}\hfill \end{array}$E9

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

Transverse stresses ${\sigma }_{xz}^k$ and ${\sigma }_{yz}^k$ 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

$\begin{array}{l}{\displaystyle \left(N_x,N_y,N_{xy},Q_x,Q_y\right)}={\displaystyle {\int }_0^h{\displaystyle \left({\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy},{\sigma }_{xz},{\sigma }_{yz}\right)}}dz\\ {\displaystyle \left(N_x^J,N_y^J,N_{xy}^J\right)}={\displaystyle {\int }_o^h{\displaystyle \left({\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy}\right)}}{\varphi }^J{\displaystyle \left(z\right)}dz\\ {\displaystyle \left(Q_x^J,Q_y^J\right)}={\displaystyle {\int }_0^h{\displaystyle \left({\sigma }_{xz},{\sigma }_{yz}\right)}}{\varphi }_{,z}^{J}dz\end{array}$E10

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

$\begin{array}{l}{\displaystyle \left\{N\right\}}={\displaystyle \left[A\right]}{\displaystyle \left\{\epsilon \right\}}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[B^J\right]}{\displaystyle \left\{{\epsilon }^J\right\}}}\hfill \\ {\displaystyle \left\{N^J\right\}}={\displaystyle \left[B^J\right]}{\displaystyle \left\{\epsilon \right\}}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[D^{IJ}\right]}{\displaystyle \left\{{\epsilon }^J\right\}}}\hfill \end{array}$E11

where ${\displaystyle \left[A\right]}$, ${\displaystyle \left[B^J\right]}$, and${\displaystyle \left[D^{IJ}\right]}$ are stiffnesses of a laminated plate

$A_{ij}={\displaystyle \sum _{J=1}^N\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}{\overline{Q}}_{ij}^{(k)}dz},(i,j=1,2,6,4,5)$E12

$\begin{array}{c}{B}_{ij}^J={\displaystyle \sum _{k=1}^N\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}{\overline{Q}}_{ij}^{(k)}{\psi }^{J}dz}(i,j=1,2,6)\\ B_{ij}^J={\displaystyle \sum _{k=1}^N\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}{\overline{Q}}_{ij}^{(k)}{\psi }_{,z}^{J}d}z(i,j=4,5)\end{array}$E13

$\begin{array}{l}{D}_{ij}^{JI}={\displaystyle \sum _{k=}^N\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}{\overline{Q}}_{ij}^{(k)}{\psi }^J{\psi }^{I}dz}(i,j=1,2,6)\\ D_{ij}^{JI}={\displaystyle \sum _{k=1}^N\underset{z_k}{\overset{z_{k+1}}{{\displaystyle \int }}}{\overline{Q}}_{ij}^{(k)}{\psi }_{,z}^J{\psi }_{,z}^{I}dz}(i,j=4,5)\end{array}$E14

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

$\begin{array}{l}{A}_{ij}={\displaystyle \sum _{k=1}^N{{\displaystyle \overline{Q}}}_{ij}^K(z_{k-1}-z_k),}(i,j=1,2,6,4,5)\hfill \\ B_{ij}^J={{\displaystyle \overline{Q}}}_{ij}^{J-1}\frac{(z_J-z_{J-1})}{2}+{{\displaystyle \overline{Q}}}_{ij}^J\frac{(z_{J+1}-z_J)}{2}(i,j=1,2,6)\hfill \\ B_{ij}^J={\overline{Q}}_{ij}^{J-1}-{\overline{Q}}_{ij}^J(i,j=4,5)\hfill \\ forJ=I\hfill \\ D_{ij}^{IJ}={{\displaystyle \overline{Q}}}_{ij}^{J-1}\frac{(z_J-z_{J-1})}{3}+{{\displaystyle \overline{Q}}}_{ij}^J\frac{(z_{J+1}-z_J)}{3}(i,j=1,2,6)\\ D_{ij}^{IJ}=\frac{{\overline{Q}}_{ij}^{J-1}}{(z_J-z_{J-1})}+\frac{{\overline{Q}}_{ij}^J}{(z_{J+1}-z_J)}(i,j=4,5)\hfill \\ forJ=I\pm 1\\ D_{ij}^{IJ}=D_{ij}^{JI}={\overline{Q}}_{ij}^J\frac{(z_{J+1}-z_J)}{6}(i,j=1,2,6)\hfill \\ D_{ij}^{IJ}=D_{ij}^{JI}=-\frac{{\overline{Q}}_{ij}^J}{(z_{J+1}-z_J)}(i,j=4,5)\hfill \end{array}$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 A₁₆ = A₂₆ = A₄₅ = B₁₆^J = B₂₆^J = B₄₅^J = D₁₆^JI = D₂₆^JI = D₄₅^JI, the equations take the form:

$\begin{array}{l}{A}_{11}{u}_{,xx}+A_{12}{v}_{,yx}+A_{66}{\displaystyle \left(u_{,yy}+v_{,xy}\right)}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[B_{11}^J{U}_{,xx}^J+B_{12}^J{V}_{,yx}^J+B_{66}^J{\displaystyle \left(U_{,yy}^J+V_{,xy}^J\right)}\right]}}=0\\ A_{12}{u}_{,xy}+A_{22}{v}_{,yy}+A_{66}{\displaystyle \left(u_{,yx}+v_{,xx}\right)}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[B_{12}^J{U}_{,xy}^J+B_{22}^J{V}_{,yy}^J+B_{66}^J{\displaystyle \left(U_{,yx}^J+V_{,xx}^J\right)}\right]}}=0\\ A_{55}{w}_{,xx}+A_{44}{w}_{,yy}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[B_{55}^J{U}_{,x}^J+B_{44}^J{V}_{,y}^J\right]}}+q=0\\ B_{11}^J{u}_{,xx}+B_{12}^J{v}_{,xy}+B_{66}^J{\displaystyle \left(u_{,yy}+v_{,xy}\right)}-B_{55}^J{w}_{,x}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[D_{11}^{JI}{U}_{,xx}^J+D_{12}^{JI}{V}_{,yx}^J+D_{66}^{JI}{\displaystyle \left(U_{,yy}^J+V_{,xy}^J\right)}-D_{55}^{JI}{U}^I\right]}}=0\\ B_{11}^J{u}_{,xy}+B_{21}^J{v}_{,yy}+B_{66}^J{\displaystyle \left(u_{,yx}+v_{,xx}\right)}-B_{44}^J{w}_{,y}+{\displaystyle \sum _{J=1}^n{\displaystyle \left[D_{12}^{JI}{U}_{,xy}^J+D_{22}^{JI}{V}_{,yy}^J+D_{66}^{JI}{\displaystyle \left(U_{,yx}^J+V_{,xx}^J\right)}-D_{44}^{JI}{V}^I\right]}}=0\end{array}$E16

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 U^J and V^J, 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

$\begin{array}{l}\begin{array}{l}u={\displaystyle \sum _{m,n}^{\infty }{X}_{mn}\text{cos}\frac{m\pi x}{a}\text{sin}\frac{n\pi y}{b}}\\ v={\displaystyle \sum _{m,n}^{\infty }{Y}_{mn}\text{sin}\frac{m\pi x}{a}\text{cos}\frac{n\pi y}{b}}\\ w={\displaystyle \sum _{m,n}^{\infty }{W}_{mn}\text{sin}\frac{m\pi x}{a}\text{sin}\frac{n\pi y}{b}}\end{array}\hfill \\ \begin{array}{l}{U}^J={\displaystyle \sum _{m,n}^{\infty }{R}_{mn}^J\text{cos}\frac{m\pi x}{a}\text{sin}\frac{n\pi y}{b}}\\ V^J={\displaystyle \sum _{m,n}^{\infty }{S}_{mn}^J\text{sin}\frac{m\pi x}{a}\text{cos}\frac{n\pi y}{b}}\end{array}\hfill \end{array}$E17

where $\alpha =\frac{m\pi }{a};\hspace{1em}\beta =\frac{n\pi }{b}\hspace{1em}J=1,\mathrm{\dots },n$

Boundary conditions for simply supported plate are met:

$\begin{array}{l}u=w=V^J=N_x=N_x^J=0x=0,aJ=1,n\\ v=w=U^J=N_y=N_y^J=0y=0,bJ=1,n\end{array}$E18

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

$q(x,y)={\displaystyle \sum _{m,n}^{\infty }{Q}_{mn}sin\frac{m\pi x}{a}sin\frac{n\pi y}{b}}$E19

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

${\displaystyle \left[\begin{array}{cc}{\displaystyle \left[K\right]}& {\displaystyle \left[K^J\right]}\\ {{\displaystyle \left[K^J\right]}}^T& {\displaystyle \left[C^{JI}\right]}\end{array}\right]}{\displaystyle \left\{\begin{array}{c}{X}_{mn}\\ Y_{mn}\\ W_{mn}\\ R_{mn}^J\\ S_{mn}^J\end{array}\right\}}={\displaystyle \left\{\begin{array}{c}0\\ 0\\ Q_{mn}\\ 0\\ 0\end{array}\right\}}$E20

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

By solving Eq. (20), the unknown coefficients X_mn, Y_mn, W_mn, R_mn, and S_mn 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:

$\begin{array}{l}{\sigma }_{xx}\left(x,y,z\right)=-{\displaystyle \sum _{m,n}^{\infty }\left\{\left[{\overline{Q}}_{11}\alpha (X_{mn}+{\displaystyle \sum _{J=1}^n{R}_{mn}^J{\varphi }^J\left(z\right)})+{\overline{Q}}_{12}\beta (Y_{mn}+{\displaystyle \sum _{J=1}^n{S}_{mn}^J{\varphi }^J\left(z\right)})\right]\text{sin}\alpha x\text{sin}\beta y\right\}}\\ {\sigma }_{yy}\left(x,y,z\right)=-{\displaystyle \sum _{m,n}^{\infty }\left\{\left[{\overline{Q}}_{12}\alpha (X_{mn}+{\displaystyle \sum _{J=1}^n{R}_{mn}^J{\varphi }^J\left(z\right)})+{\overline{Q}}_{22}\beta (Y_{mn}+{\displaystyle \sum _{J=1}^n{S}_{mn}^J{\varphi }^J\left(z\right)})\right]\text{sin}\alpha x\text{sin}\beta y\right\}}\\ \sigma \left(x,y,z\right)=-{\overline{Q}}_{66}{\displaystyle \sum _{m,n}^{\infty }\left\{\left[\beta (X_{mn}+{\displaystyle \sum _{J=1}^n{R}_{mn}^J{\varphi }^J\left(z\right)})+\alpha (Y_{mn}+{\displaystyle \sum _{J=1}^{\infty }{S}_{mn}^J{\varphi }^J\left(z\right)})\right]\text{cos}\alpha x\text{sin}\beta y\right\}}\end{array}$E21

Shear stresses σ_xz and σ_y are constant through the layer J and can be obtained with

$\begin{array}{l}{\sigma }_{xz}^j{\displaystyle \left(x,y,z\right)}={\displaystyle \sum _{m,n}^{\infty }{\displaystyle \left\{{\overline{Q}}_{45}^j{\displaystyle \left\{\frac{1}{h_k}{\displaystyle \left[{\displaystyle \left(S_{mn}^{j+1}-S_{mn}^j\right)}+\beta W_{mn}\right]}\text{sin}\alpha x\text{ cos}\beta y\right\}}+{\overline{Q}}_{55}^j{\displaystyle \left\{\frac{1}{h_k}{\displaystyle \left[{\displaystyle \left(R_{mn}^{j+1}-R_{mn}^j\right)}+\alpha W_{mn}\right]}\text{cos}\alpha x\text{ sin}\beta y\right\}}\right\}}}\\ {\sigma }_{xz}^j{\displaystyle \left(x,y,z\right)}={\displaystyle \sum _{m,n}^{\infty }{\displaystyle \left\{{\overline{Q}}_{44}^j{\displaystyle \left\{\frac{1}{h_k}{\displaystyle \left[{\displaystyle \left(S_{mn}^{j+1}-S_{mn}^j\right)}+\beta W_{mn}\right]}\text{sin}\alpha x\text{ cos}\beta y\right\}}+{\overline{Q}}_{45}^j{\displaystyle \left\{\frac{1}{h_k}{\displaystyle \left[{\displaystyle \left(R_{mn}^{j+1}-R_{mn}^j\right)}+\alpha W_{mn}\right]}\text{cos}\alpha x\text{ sin}\beta y\right\}}\right\}}}\end{array}$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:

${\overline{\sigma }}_{xx}=\frac{1}{q{s}^2}{\sigma }_{xx},{\overline{\sigma }}_{yy}=\frac{1}{q{s}^2}{\sigma }_{yy},{\overline{\sigma }}_{xz}=\frac{1}{qs}{\sigma }_{xz},{\overline{\sigma }}_{yz}=\frac{1}{qs}{\sigma }_{yz}$E23

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 E₁, E₂, G₁₂, G₂₁, G₂₃, ν_12, and ν₁₃, 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)

${\overline{\sigma }}_{xx}={\overline{\sigma }}_{xx}(A,A,z),\hspace{1em}{\overline{\sigma }}_{yy}={\overline{\sigma }}_{yy}(A,A,z),\hspace{1em}{\overline{\sigma }}_{xz}={\overline{\sigma }}_{xz}(B,A,z),\hspace{1em}{\overline{\sigma }}_{yz}={\overline{\sigma }}_{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 E₁/E₂=25, E₂=1, G₁₂=G₁₃=0.5, G₂₃=0.2, ν₁₂=ν₁₃=0.25. The change of dimensionless normal stresses σ_xx, σ_yy and shear stresses σ_yz, σ_xz are displayed in Figures 2 and 3.

Figure 2.

${\overline{\sigma }}_{xx}$ and ${\overline{\sigma }}_{yy}$ through the thickness of the plate 0^o/90^o/0^o/90^o/0^o, bi‐sinusoidal load.

Figure 3.

${\overline{\sigma }}_{yz}$ and ${\overline{\sigma }}_{xz}$ through the thickness of the plate 0^o/90^o/0^o/90^o/0^o, bi‐sinusoidal load.

Example 2: The square laminated composite plate with six layers 0^o/90^o/0^o/90^o/0^o/90^o 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.

Figure 4.

${\overline{\sigma }}_{xx}$ and ${\overline{\sigma }}_{yy}$ through the thickness of the plate 0^o/90^o/0^o/90^o/0^o/90°.

Figure 5.

${\overline{\sigma }}_{yz}$ and ${\overline{\sigma }}_{xz}$ through the thickness of the plate 0^o/90^o/0^o/90^o/0^o/90^o.

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 u_i, v_i, w_i, u_i ^j, v_i ^j

${\displaystyle \left(u,v,w,u_j,v_j\right)}={\displaystyle \sum _{i=1}^m{\displaystyle \left(u_i,v_i,w_i,u_i^j,v_i^j\right)}{\varphi }_i}$E25

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

$\begin{array}{l}x={\displaystyle \sum _{i=1}^m{x}_i}{\varphi }_i{\displaystyle \left(\xi ,\eta \right)}\\ y={\displaystyle \sum _{i=1}^m{y}_i{\varphi }_i{\displaystyle \left(\xi ,\eta \right)}}\end{array}$E26

where x, y are point coordinates and x_i, y_i are coordinates of the finite element.

According to Eqs. (6) and (25), matrices H and $\overline{H}$ are defined as follows:

${{\displaystyle \left[H\right]}}_{{\displaystyle \left(5x3m\right)}}={\displaystyle \left[\begin{array}{cccccccccc}\frac{\partial {\varphi }_1}{\partial x}& 0& 0& \frac{\partial {\varphi }_2}{\partial x}& 0& 0& \dots & \frac{\partial {\varphi }_m}{\partial x}& 0& 0\\ 0& \frac{\partial {\varphi }_1}{\partial y}& 0& 0& \frac{\partial {\varphi }_2}{\partial y}& 0& \dots & 0& \frac{\partial {\varphi }_m}{\partial y}& 0\\ \frac{\partial {\varphi }_1}{\partial y}& \frac{\partial {\varphi }_1}{\partial x}& 0& \frac{\partial {\varphi }_2}{\partial y}& \frac{\partial {\varphi }_2}{\partial x}& 0& \dots & \frac{\partial {\varphi }_m}{\partial y}& \frac{\partial {\varphi }_m}{\partial x}& 0\\ 0& 0& \frac{\partial {\varphi }_1}{\partial x}& 0& 0& \frac{\partial {\varphi }_2}{\partial x}& \dots & 0& 0& \frac{\partial {\varphi }_m}{\partial x}\\ 0& 0& \frac{\partial {\varphi }_1}{\partial y}& 0& 0& \frac{\partial {\varphi }_2}{\partial y}& \dots & 0& 0& \frac{\partial {\varphi }_m}{\partial y}\end{array}\right]}$E27

${{\displaystyle \left[\overline{H}\right]}}_{{\displaystyle \left(5x2m\right)}}={\displaystyle \left[\begin{array}{ccccccc}\frac{\partial {\varphi }_1}{\partial x}& 0& \frac{\partial {\varphi }_2}{\partial x}& 0& \dots & \frac{\partial {\varphi }_m}{\partial x}& 0\\ 0& \frac{\partial {\varphi }_1}{\partial y}& 0& \frac{\partial {\varphi }_2}{\partial y}& \dots & 0& \frac{\partial {\varphi }_m}{\partial y}\\ \frac{\partial {\varphi }_1}{\partial y}& \frac{\partial {\varphi }_1}{\partial x}& \frac{\partial {\varphi }_2}{\partial y}& \frac{\partial {\varphi }_2}{\partial x}& \dots & \frac{\partial {\varphi }_m}{\partial y}& \frac{\partial {\varphi }_m}{\partial x}\\ {\varphi }_1& 0& {\varphi }_2& 0& \dots & {\varphi }_m& 0\\ 0& {\varphi }_1& 0& {\varphi }_2& \dots & 0& {\varphi }_m\end{array}\right]}$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:

${\displaystyle \left\{\begin{array}{c}\frac{\partial {\varphi }_i}{\partial x}\\ \frac{\partial {\varphi }_i}{\partial y}\end{array}\right\}}={{\displaystyle \left[\begin{array}{cc}x{,}_{\xi }& y{,}_{\xi }\\ x{,}_{\eta }& y{,}_{\eta }\end{array}\right]}}^{-1}{\displaystyle \left\{\begin{array}{c}\varphi {,}_{\xi }^i\\ \varphi {,}_{\eta }^i\end{array}\right\}}={{\displaystyle \left[J\right]}}^{-1}{\displaystyle \left\{\begin{array}{c}\varphi {,}_{\xi }^i\\ \varphi {,}_{\eta }^i\end{array}\right\}}$E29

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

${\displaystyle \left[J\right]}={\displaystyle \left[\begin{array}{cc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }\end{array}\right]}$E30

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

${\displaystyle \left\{\begin{array}{c}\varphi {,}_{xx}\\ \varphi {,}_{yy}\\ \varphi {,}_{xy}\end{array}\right\}}={{\displaystyle \left[J_1\right]}}^{-1}{\displaystyle \left({\displaystyle \left\{\begin{array}{c}{\varphi }_i{,}_{\xi \xi }\\ {\varphi }_i{,}_{\eta \eta }\\ {\varphi }_i{,}_{\xi \eta }\end{array}\right\}}-{\displaystyle \left[J_2\right]}{\displaystyle \left\{\begin{array}{c}{\varphi }_i{,}_x\\ {\varphi }_i{,}_y\end{array}\right\}}\right)}$E31

where

${\displaystyle \left[J_1\right]}={\displaystyle \left[\begin{array}{ccc}x{,}_{\xi }^2& y{,}_{\xi }^2& 2y{,}_{\xi }^{2}x{,}_{\xi }^2\\ x{,}_{\eta }^2& y{,}_{\eta }^2& 2y{,}_{\eta }^{2}x{,}_{\eta }^2\\ x{,}_{\eta \xi }& y{,}_{\eta \xi }& x{,}_{\eta \xi }+y{,}_{\eta \xi }\end{array}\right]}$E32

${\displaystyle \left[J_2\right]}={\displaystyle \left[\begin{array}{cc}x{,}_{\xi \xi }^{}& y{,}_{\xi \xi }^{}\\ x{,}_{\eta \eta }^{}& y{,}_{\eta \eta }^{}\\ x{,}_{\eta \xi }& y{,}_{\eta \xi }\end{array}\right]}$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

$\delta {\Pi }_e={\displaystyle \underset{V_e}{\int }{\epsilon }^T\overline{Q}\delta \epsilon dV}$E34

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

${\displaystyle \left\{\epsilon \right\}}={\displaystyle \left[H\right]}{\displaystyle \left\{\Delta \right\}}+{\displaystyle \left[\overline{H}\right]}{\displaystyle \sum _{j=1}^n{\varphi }^J{\displaystyle \left\{{\Delta }^J\right\}}}$E35

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

${\varphi }^J{\displaystyle \left(z\right)}={\displaystyle \left[\begin{array}{c}{\psi }^j(z)\\ {\psi }^j(z){,}_z\end{array}\right]}$E36

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

$\delta {\Pi }_e={\displaystyle \underset{V_e}{\int }{\displaystyle \left({\Delta }^T{H}^T{\overline{H}}^T{\displaystyle \sum _{J=1}^n{\varphi }^J{\Delta }^{JT}}\right)}\cdot \overline{Q}\cdot {\displaystyle \left(H\delta \Delta +\overline{H}{\displaystyle \sum _{I=1}^n{\varphi }^I\delta {\Delta }^I}\right)}dV}$E37

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

$\begin{array}{l}{k}_e={\displaystyle \underset{V_e}{\iiint }(H^T\overline{Q}H+{\displaystyle \sum _{J=1}^n{H}^T\overline{Q}{\varphi }^J\overline{H}+{\displaystyle \sum _{I=1}^n{\overline{H}}^T\overline{Q}{\varphi }^{I}H}+{\displaystyle \sum _{I=1}^n{\displaystyle \sum _{J=1}^n{\overline{H}}^T\overline{Q}{\varphi }^I{\varphi }^J\overline{H})}}}dV=}\\ {\displaystyle \underset{A_e}{\iint }(H^{T}AH+{\displaystyle \sum _{J=1}^n{H}^T{B}^J\overline{H}+{\displaystyle \sum _{I=1}^n{{\displaystyle \overline{H}}}_{}^T{B}^{I}H}+{\displaystyle \sum _{I=1}^n{\displaystyle \sum _{J=1}^n{{\displaystyle \overline{H}}}^T{D}^{IJ}\overline{H})}}}dA}\end{array}$E38

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

$k_e={\displaystyle \left[\begin{array}{ccccc}{k}^{11}& k_1^{12}& k_2^{12}& \dots & k_n^{12}\\ k_1^{21}& k_{11}^{22}& k_{12}^{22}& \dots & k_{1n}^{22}\\ k_2^{21}& k_{12}^{22}& k_{22}^{22}& \dots & k_{2n}^{22}\\ \dots & \dots & \dots & \dots & \dots \\ k_n^{21}& k_{n1}^{22}& k_{n2}^{22}& \dots & k_{nn}^{22}\end{array}\right]}$E39

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

$\begin{array}{l}{k}^{11}={\displaystyle \underset{A_e}{\iint }{H}^{T}AHdA}i=1,\mathrm{\dots },n\\ k_i^{12}={\displaystyle \underset{A_e}{\iint }{H}^T{B}^i\overline{H}dA}i=1,\mathrm{\dots },n\\ k_i^{21}={\displaystyle \underset{A_e}{\iint }{\overline{H}}^T{B}^{i}HdA}i=1,\mathrm{\dots },n\\ k_{ij}^{22}={\displaystyle \underset{A_e}{\iint }{\overline{H}}^T{D}^{ij}\overline{H}dA}i,j=1,\mathrm{\dots },n\end{array}$E40

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

$dA=dxdy=\det Jd\xi d\eta$E41

as follows:

${\displaystyle \underset{A_e}{\iint }(\dots )dA}={\displaystyle \underset{=1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}(\dots )\det Jd\xi d\eta }}$E42

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

$|J|=\left|\begin{array}{cc}{\displaystyle \sum _{i=1}^m{x}_i\frac{\partial \varphi }{\partial \xi }}& {\displaystyle \sum _{i=1}^m{y}_i\frac{\partial \varphi }{\partial \xi }}\\ {\displaystyle \sum _{i=1}^m{x}_i\frac{\partial \varphi }{\partial \eta }}& {\displaystyle \sum _{i=1}^m{y}_i\frac{\partial \varphi }{\partial \eta }}\end{array}\right|$E43

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

$\begin{array}{l}{k}^{11}={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}{H}^{T}AH\det Jd\xi d\eta }}i=1,\mathrm{\dots },n\\ k_i^{12}={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}}}{H}^T{B}^i{\overline{H}}_{}^{}\det Jd\xi d\eta i=1,\mathrm{\dots },n\\ k_i^{21}={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}{\overline{H}}^T{B}^{iT}H\det Jd\xi d\eta }}i=1,\mathrm{\dots },n\\ k_{ij}^{22}={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}{\overline{H}}^T{D}^{ij}\overline{H}\det Jd\xi d\eta }}i,j=1,\mathrm{\dots },n\end{array}$E44

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.

Figure 6.

Rectangular finite element E8–L3.

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(\xi ,\eta )={\displaystyle \underset{-1}{\overset{+1}{\int }}{\displaystyle \underset{-1}{\overset{+1}{\int }}f(\xi ,\eta )d\xi d\eta ={\displaystyle \sum _{i=1}^m{\displaystyle \sum _{j=1}^m{\omega }_i{\omega }_{j}f(\overline{\xi },\overline{\eta })}}}}$E45

where ω_i and ω_j are the weights of the integration, and $(\overline{\xi },\overline{\eta })$ 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

$\partial \Pi =\delta {\mathit{U}}^{\mathit{T}}\mathit{K}\mathit{U}-\delta {\mathit{U}}^{\mathit{T}}\mathit{F}$E46

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

$\mathit{K}\mathit{U}=\mathit{F}$E47

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 E₁=25, E₂=1, G₁₂=G₁₃=0.5, G₂₃=0.2, ν₁₂=ν₁₃=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:

$\begin{array}{l}{\overline{\sigma }}_{xx}{\displaystyle \left(0.4471688a;0.4471688a\right)}{\overline{\sigma }}_{yy}{\displaystyle \left(0.4471688a;0.4471688a\right)}{\overline{\sigma }}_{xy}{\displaystyle \left(0.05283122a;0.05283122a\right)}\\ {\overline{\sigma }}_{yz}{\displaystyle \left(0.4471688a;0.9471688a\right)}{\overline{\sigma }}_{xz}{\displaystyle \left(0.9471688a;0.4471688a\right)}\end{array}$

Changing of dimensionless stresses through the thickness of the laminated plate is presented in Figures 7–9.

Figure 7.

${\overline{\sigma }}_{xx}$ and ${\overline{\sigma }}_{yy}$ through the thickness of the plate.

Figure 8.

${\overline{\sigma }}_{xy}$ and ${\overline{\sigma }}_{xz}$ through the thickness of the plate.

Figure 9.

${\overline{\sigma }}_{yz}$ through the thickness of the plate.

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 u^J(x,y,z) and v^J(x,y,z) and deflections w^J(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.