Chapter 3 - Plates and shells

III.1 Overview

Plates and shells are common structures in the field of applied mechanics, their main property being that they support transverse loading through bending action (and is some cases as we will see, through membrane and shear effects). Finite element formulation for plate and shells has a higher degree of complexity than other types of elements. The complexity arises from the nature of high order differential equations. To overcome the complex formulations, different solutions have been proposed. First, the classical plate bending theory is presented along with its finite element formulation. Then other types of plate theories follow, such as shear deformable plate elements, with displacement degrees of freedom only, mixed elements and hybrid elements. Then the type of shell derived, is the one which is the combination of plate bending and plane stress elements, with a proper coordinate transformation from local to global coordinate axis, while using the drilling degree of freedom to bypass singularity problem in the transformation matrix. The second shell element is developed by degenerating from 3-D solid elements, but for the sake of simplicity we focus on the first type of shell only.

III.2 Classical plate theory (Kirchoff-Love)

The basic assumptions for the classical Kirchhoff plate bending theory are very similar to those for the Euler-Bernoulli beam theory. One of the most important assumptions for both theories is that a straight line normal to the midplane in the un-deformed case remains normal after deformation. This leads to the fact that transverse deformation is not going to be taken into account, or as we might say: neglected.

In-plane displacements $ u $ and $ v $ can be expressed as

$$u=-z \frac{\partial w}{\partial x}$$

$$v=-z \frac{\partial w}{\partial y}$$

Where $ x $ and $ y $ are the in-plane axes located at the midplane of the plate, and $ z $ is perpendicular to the plane defined by these (in direction of the thickness). As follows $ u $ and $ v $ are displacements along, respectively, the $ x $- and $ y $- axes, while $ w $ is the transverse displacement (deflection) along the $ z $-axis.

Because we neglect the transverse shear deformation, in-plane strains can be written in terms of displacements

$$\{\begin{array}{lll}{\epsilon }_x& {\epsilon }_y& {\gamma }_{xy}\}\end{array}=-z\left\{\begin{array}{ccc}{\kappa }_x& {\kappa }_y& {\kappa }_{xy}\end{array}\right\}$$

where

$$\{\kappa {\}}^T=\left\{\begin{array}{lll}{\kappa }_x& {\kappa }_y& {\kappa }_{xy}\end{array}\right\}=z\left\{\begin{array}{lll}\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{x}^2}& \frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }{y}^2}& 2\frac{{\mathrm{\partial }}^{2}w}{\mathrm{\partial }x\mathrm{\partial }y}\end{array}\right\}$$

is called curvature.

Assuming plane stress condition for plate bending and substituting (3) and (4) into the constitutive equation. This states the relationship between the stresses and strains. For an isotropic material, the constitutive equation becomes

$$\{\sigma \}=[D]\{\epsilon \}$$

By substituting (3) and (4)

$$\{\sigma \}=-z[D]\{\kappa \}$$

in which

$$\{\sigma \}=\{\begin{array}{lll}{\sigma }_x& {\sigma }_y& {{\tau }_{xy}\}}^T\end{array}$$

and the material property matrix $ [D] $ becomes

$$[D]=\frac{E}{1-v^2}\left[\begin{array}{ccc}1& v& 0\\ v& 1& 0\\ 0& 10& \frac{1-v}{2}\end{array}\right]$$

Moments are defined as

$$\{M\}={\int }_{h/2}^{h/2}\{\sigma \}zdz$$

where

and $ h $ is the plate thickness. Substituting Eq. (6) into Eq. (9) gives the relationship between moment and curvature.

$$\{M\}=-[\tilde{D}]\{\kappa \}$$

where

$$[\widetilde{D}]=\frac{h^{3}}{12}[D]$$

Equilibrium equations are obtained from the free-body diagram. Moment equilibriums about the $ y $- and $ x $-axes and force equilibrium about the $ z $-axis yield, after neglecting the higher order terms,

$$\frac{\partial M_{x}}{\partial x}+\frac{\partial M_{x y}}{\partial y}-Q_{x}=0$$

$$\frac{\partial M_{x y}}{\partial x}+\frac{\partial M_{y}}{\partial y}-Q_{y}=0$$

$$\frac{\partial Q_{x}}{\partial x}+\frac{\partial Q_{y}}{\partial y}+p=0$$

where $ Q_{x} $ and $ Q_{y} $ are the shear forces and $ p $ is the distributed pressure loading. Elimination of the shear forces from Eqs. (12) to (14) gives

$$\frac{\partial^{2} M_{x}}{\partial x^{2}}+2 \frac{\partial^{2} M_{x y}}{\partial x \partial y}+\frac{\partial^{2} M_{y}}{\partial y^{2}}+p=0$$

Combining Eqs. (4), (10) and (15) finally produces the biharmonic governing equation for plate bending in terms of transverse displacement $ w $.

$$\frac{\partial^{4} w}{\partial x^{4}}+2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}}+\frac{\partial^{4} w}{\partial y^{4}}=\frac{p}{D_{r}}$$

where $ D_{e}=\frac{E h^{3}}{12\left(1-v^{2}\right)} $ is the plate rigidity.

This is often written as

$$\nabla^{4} w=\frac{q}{D_{r}}$$

Summing up, the deflection theory of plates attributed to Kirchhoff is based on the following assumptions:

• The x–y plane coincides with the middle plane of the plate in the un-deformed geometry.
• The lateral dimension of the plate is at least 10 times its thickness.
• The vertical displacement of any point of the plate can be taken equal to that of the point (below or above it) in the middle plane.
• A vertical element of the plate before bending remains perpendicular to the middle surface of the plate after bending.
• Strains are small: deflections are less than the order of (1/100) of the span length.
• The strain of the middle surface is zero or negligible.

III.3 Classical plate bending elements

III.3.1 Three-node element

Deriving a three-node element based on the classical plate bending theory. Shown in the following figure. Each node has three degrees of freedom: displacement in $ z $ direction – $ w $; a rotation around the $ x $-axis, $ w_{y} $, (derivative of $ w $ with respect to $ y $); and a rotation about the $ y $-axis, $ w_{x} $ (derivative of $ w $ with respect to $ x $).

The displacement function $ w $ is assumed to be

$$\begin{array}{c}{w(x, y)=a_{1}+a_{2} x+a_{3} x+a_{4} x^{2}+a_{5} x y+a_{6} y^{2}+a_{7} x^{3}+a_{8}\left(x^{2} y+x y^{2}\right)+a_{9} y^{3}=} \\ {=[X]{a}}\end{array}$$

where

$$[X]=\left[\begin{array}{llllllll}{1} & {x} & {y} & {x^{2}} & {x y} & {y^{2}} & {x^{3}} & {\left(x^{2} y+x y^{2}\right)} & {y^{3}}\end{array}\right]$$

and

$$\{a\}=\{\begin{array}{llllllll}{a}_1& a_2& a_3& a_4& a_5& a_6& a_7& {a_8\}}^T\end{array}$$

With constants $ a_{i} $ are to be replaced by nodal variables $ w $, $ w_{x} $ and $ w_{y} $. Taking derivatives of the displacement function with respect to $ x $ and $ y $ yields

$$\frac{\partial w}{\partial x}=a_{2}+2 a_{4} x+a_{5} y+3 a_{7} x^{2}+a_{8}\left(2 x y+y^{2}\right)$$

$$\frac{\partial w}{\partial y}=a_{3}+a_{5} x+2 a_{6} y+2 a_{8}\left(x^{2}+2 x y\right)+3 a_{9} y^{2}$$

Evaluation of equations from (18) through (22) at the three nodal points gives the following equation:

$$\{d\}=[\tilde{X}]\{a\}$$

where

$$\{d\}={\left\{w_1{\left(w_x\right)}_1{\left(w_y\right)}_1{w}_2{\left(w_x\right)}_2{\left(w_y\right)}_2{w}_3{\left(w_x\right)}_3{\left(w_y\right)}_3\right\}}^T$$

and

$$[\tilde{X}]=\left[\begin{array}{ccccccccc}1& x_1& y_1& x_1^2& x_1{y}_1& y_1^2& x_1^3& x_1^2{y}_1+x_1{y}_1^2& y_1^3\\ 0& 1& 0& 2x_1& y_1& 0& 3x_1^2& 2x_1{y}_1+y_1^2& 0\\ 0& 0& 1& 0& x_1& 2y_1& 0& x_1^2+2x_1{y}_1& 3y_1^2\\ 1& x_2& y_2& x_2^2& x_2{y}_2& y_2^2& x_2^3& x_2^2{y}_2+x_2{y}_2^2& y_2^3\\ 0& 1& 0& 2x_2& y_2& 0& 3x_2^2& 2x_2{y}_2+y_2^2& 0\\ 0& 0& 1& 0& x_2& 2y_2& 0& x_2^2+2x_2{y}_2& 3y_2^2\\ 1& x_3& y_3& x_3^2& x_3{y}_3& y_3^2& x_3^3& x_3^2{y}_3+x_3{y}_3^2& y_3^3\\ 0& 1& 0& 2x_3& y_3& 0& 3x_3^2& 2x_3{y}_3+y_3^2& 0\\ 0& 0& 1& 0& x_3& 2y_3& 0& x_3^2+2x_3{y}_2& 3y_3^2\end{array}\right]$$

Inverting equations (23) and substituting the results into Eq. (18) yields the following

$$\{w\}=[H]\{d\}$$

Where the row vector of shape functions of size $ 9 * 1 $ is computed from

$$[H]=[X][X]^{-1}$$

In-plane strains are computed from Eq. (26) as

$$\{\epsilon \}=[B]\{d\}$$

where

$$\{\epsilon \}=\{\begin{array}{lll}{\epsilon }_x& {\epsilon }_y& {{\gamma }_{xy}\}}^T\end{array}$$

$$[B]=-z[L][\tilde{X}]^{-1}$$

$$[L]=\left[\begin{array}{ccccccccc}0& 0& 0& 2& 0& 10& 6x& 2y& 0\\ 0& 0& 0& 0& 0& 2& 0& 2x& 6y\\ 0& 0& 0& 0& 2& 0& 0& 4(x+y)& 0\end{array}\right]$$

Substitution of the strain-nodal displacement relation, Eq. (28) into the expression for the element stiffness matrix,

$$\left[K_{e}\right]=\int_{\Omega^{e}}[B]^{T}[D][B] d \Omega=[B]^{T}[D][B] A$$

yields

$$\begin{array}{l}{\quad\left[K^{e}\right]=\int_{\Omega^{e}} \int_{z}[B]^{T}[D][B] d z d \Omega} \\ {=[\tilde{X}]^{-T} \int_{\Omega^{e}} \int_{z} Z^{2}[L]^{T}[D][L] d z d \Omega[\tilde{X}]^{-1}} \\ {\quad=[\tilde{X}]^{-T} \int_{\Omega^{e}}[L]^{T}[\widetilde{D}][L] d \Omega[\tilde{X}]^{-1}}\end{array}$$

where

$$[\widetilde{D}]=\frac{h^{3}}{12}[D]$$

Here $ [D] $ is the constitutive matrix of the plane stress condition, $ \Omega^{e} $ is the two-dimensional element domain in the $ x y $-plane, and $ h $ is the thickness of the plate. The element domain is the triangular shape as seen in the previous figure. The element stiffness matrix is of size $ 9 × 9 $ and the corresponding element nodal vector is given in Eq. (24).

III.3.2 Four-node element

Consider the rectangular plate element shown in the following figure. The element has four nodes and 12 dof in total. A trial function for the unknown $ w(x,y) $ will contain 12 parameters:

$$w(x, y)=\alpha_{1}+\alpha_{2} x+\alpha_{3} y+\alpha_{4} x^{2}+\alpha_{5} x y+\alpha_{6} y^{2}+\alpha_{7} x^{3}+\alpha_{8} x^{2} y+\alpha_{9} x y^{2}+ \ \alpha_{10} y^{3}+\alpha_{11} y x^{3}+\alpha_{12} x y^{3}$$

$$\frac{\partial w}{\partial x}=\alpha_{2}+2 \alpha_{4} x+\alpha_{5} y+3 \alpha_{7} x^{2}+2 \alpha_{8} x y+\alpha_{9} y^{2}+3 \alpha_{11} x^{2} y+\alpha_{12} y^{3}$$

$$\frac{\partial w}{\partial y}=\alpha_{3}+\alpha_{5} x+2 \alpha_{6} y+\alpha_{8} x^{2}+2 \alpha_{9} y x+3 \alpha_{10} y^{2}+\alpha_{11} x^{3}+3 \alpha_{12} y^{2} x$$

It can be seen that both $ w (x,y) $ and its derivatives are defined by cubic polynomials. As a cubic is uniquely defined by four constants, the two end values of the displacements and slopes will therefore define the displacements uniquely along any boundary. However, this is not the case for the derivatives, since only two end values of the slopes exist, the cubic is not specified uniquely. And in general, a discontinuity of normal slope will occur. The function is therefore called “non-conforming.”

III.4 Shear deformable plate element (Reissner-Mindlin plate)

The Reissner-Mindlin plate theory includes the effect of transverse shear deformation like the Timoshenko beam theory. Hence, a plane normal to the midplane of the plate before deformation does not remain normal to the midplane any longer after deformation. The internal energy expression for the shear deformable plate should include transverse shear energy as well as bending energy. The internal energy is expressed as

$$U=\frac{1}{2}{\int }_V{\left\{{\sigma }_b\right\}}^T\left\{{ϵ}_b\right\}dV+\frac{\kappa }{2}{\int }_V{\left\{{\sigma }_s\right\}}^T\left\{{\epsilon }_s\right\}dV$$

where

$$\left\{{\sigma }_b\right\}=\{\begin{array}{lll}{\sigma }_x& {\sigma }_y& {{\tau }_{xy}\}}^T\end{array}$$

$$\left\{{ϵ}_b\right\}={\left\{{ϵ}_x{ϵ}_y{\gamma }_{xy}\right\}}^T$$

are the bending stress and strain components while

$$\left\{{\sigma }_s\right\}={\left\{{\tau }_{xz}{\tau }_{yz}\right\}}^T$$

$$\left\{{ϵ}_s\right\}={\left\{{\gamma }_{xz}{\gamma }_{yz}\right\}}^T$$

are the transverse shear components. In addition, $ \kappa $ is the shear energy correction factor and equal to $ \frac{5}{6} $ .

Substitution of the constitutive equations for both bending and shear components yields

$$U=\frac{1}{2}{\int }_V{\left\{{ϵ}_b\right\}}^T\left[D_b\right]\left\{{ϵ}_b\right\}dV+\frac{\kappa }{2}{\int }_V{\left\{{ϵ}_s\right\}}^T\left[D_s\right]\left\{{ϵ}_s\right\}dV$$

in which

$$\left[D_b\right]=\frac{E}{1-v^2}\left[\begin{array}{ccc}1& v& 0\\ v& 1& 0\\ 0& 0& \frac{1-v}{2}\end{array}\right]$$

is the constitutive equation for the plane stress condition and

$$\left[D_s\right]=\left[\begin{array}{ll}G& 0\\ 0& G\end{array}\right]$$

Further, $ V $ is the three-dimensional domain which is equal to $ 𝑑Ω×𝑑𝑧 $. The $ xyz $ coordinate axis is the same as in the figure.

In order to get the derivate element stiffness matrix for the shear deformable plate, we need to express the strains in terms of nodal variables. The in-plane displacements are given as

$$u=-z \theta_{x}(x, y)$$

$$v=-z \theta_{y}(x, y)$$

and the transverse displacement is

$$w=w(x, y)$$

where $ \theta_{x} $ and $ \theta_{y} $ are rotations of the midplane about the $ y $ and $ x $ axes, respectively. The midplane is assumed to have no in-plane deformation. In the case of the shear deformable plate,

$$\theta_{x}=\frac{\partial w}{\partial x}-\gamma_{x z}$$

$$\theta_{y}=\frac{\partial w}{\partial y}-\gamma_{y z}$$

where $ \gamma $ is the angle caused by the shear deformation.

Because the transverse displacement $ w $ and slope $ \theta $ are independent, we need shape functions to interpolate them independently. As a result, the shear deformable element requires $ C^{0} $ compatibility. Isoparametric shape functions are used for the plate element formulation as we will see in the following. The transverse displacement and slopes are interpolated as

$$w=\sum_{i=1}^{n} H_{i}(\xi, \eta) w_{i}$$

$$\theta_{x}=\sum_{i=1}^{n} H_{i}(\xi, \eta)\left(\theta_{x}\right)_{i}$$

$$\theta_{y}=\sum_{i=1}^{n} H_{i}(\xi, \eta)\left(\theta_{y}\right)_{i}$$

Here $ n $ is the number of nodes per element and the same shape functions are used for the displacement and slope interpolations. As for the sake of simplicity, bilinear isoparametric shape functions are used. However, higher-order shape functions can and are used in the same manner. Bending and shear strains are computed from displacements.

$$\left\{{\epsilon }_b\right\}=-z\left[B_b\right]\left\{d^e\right\}$$

$$\left\{{\epsilon }_s\right\}=\left[B_s\right]\left\{d^e\right\}$$

where

$$\left[B_b\right]=\left[\begin{array}{cccccccccccc}\frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }x}& 0& 0& \frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }x}& 0& 0& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }x}& 0& 0& \frac{\mathrm{\partial }{H}_4}{\mathrm{\partial }x}& 0& 0\\ 0& \frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }y}& 0& 0& \frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }y}& 0& 0& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }y}& 0& 0& \frac{\mathrm{\partial }{H}_4}{\mathrm{\partial }y}& 0\\ \frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }y}& \frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }x}& 0& \frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }y}& \frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }x}& 0& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }y}& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }x}& 0& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }y}& \frac{\mathrm{\partial }{H}_4}{\mathrm{\partial }x}& 0\end{array}\right]$$

$$\left[B_s\right]=\left[\begin{array}{cccccccccccc}-H_1& 0& \frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }x}& -H_2& 0& \frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }x}& -H_3& 0& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }x}& -H_4& 0& \frac{\mathrm{\partial }{H}_4}{\mathrm{\partial }x}\\ 0& -H_1& \frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }y}& 0& -H_2& \frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }y}& 0& -H_3& \frac{\mathrm{\partial }{H}_3}{\mathrm{\partial }y}& 0& -H_4& \frac{\mathrm{\partial }{H}_4}{\mathrm{\partial }y}\end{array}\right]$$

$$\begin{array}{cccccccccccc}\{{\left({\theta }_x\right)}_1& {\left({\theta }_y\right)}_1& w_1& {\left({\theta }_x\right)}_2& {\left({\theta }_y\right)}_2& w_2& {\left({\theta }_x\right)}_3& {\left({\theta }_y\right)}_3& w_3& {\left({\theta }_x\right)}_4& {\left({\theta }_y\right)}_4& w_4\}\end{array}$$

Substitution of Eqs. (54) and (55) into the energy expression Eq. (43) yields for each plate element

$$U=\frac{1}{2}{\left\{d^e\right\}}^T{\int }_{{\mathrm{\Omega }}^e}{\int }_z{\left[B_b\right]}^T\left[D_b\right]\left[B_b\right]dzd\mathrm{\Omega }\left\{d^e\right\}+\frac{\kappa }{2}{\left\{d^e\right\}}^T{\int }_{{\mathrm{\Omega }}^e}{\int }_z{\left[B_s\right]}^T\left[D_s\right]\left[B_s\right]dzd\mathrm{\Omega }\left\{d^e\right\}$$

As a result, the element stiffness matrix for plate bending in case of thick plates can be expressed as

$$\left[K^e\right]=\frac{h^3}{12}{\int }_{{\mathrm{\Omega }}^e}{\left[B_b\right]}^T\left[D_b\right]\left[B_b\right]d\mathrm{\Omega }+\kappa h{\int }_{{\mathrm{\Omega }}^e}{\left[B_s\right]}^T\left[D_s\right]\left[B_s\right]d\mathrm{\Omega }$$

where $ h $ is the plate thickness.

III.5 Plate element with displacement degrees of freedom

We have the plate bending developed in this section represented on the following figure, where $ x $, $ y $ and $ z $ describe global coordinates of the plate and $ u $, $ v $ and $ w $ are the displacements; $ h $ is the plate thickness. The $ xy $ plane is parallel to the midplane prior to deflection.

The displacement of any point can be described with

$$u=u(x, y, z)$$

$$v=v(x, y, z)$$

$$w=w(x, y, z)$$

Meaning that the in-plane displacements $ u $ and $ v $ vary through the plate thickness as well as within the $ xy $-plane while the transverse displacements $ w $ remains constant through the plate thickness. In order to interpolate displacements with the help of shape functions and nodal displacements, two kinds of interpolations are needed: one interpolation within the $ xy $-plane and the other in the $ z $-axis.

In the $ xy $-plane during interpolation, shape functions $ N_{i}(x,y) $ are used where subscript $ i $ varies depending on the number of nodes on the $ xy $-plane. On the other hand, shape functions $ H_{j}(z) $ are used for interpolation along the $ z $-axis, where subscript $ j $ varies depending on the number of nodes along the plate thickness. Because two in-plane displacements are functions of $ x $, $ y $ and $ z $, both shape functions are used, while transversal displacement uses only shape functions $ N_{i}(x,y) $. With the help of isoparametric elements with mapping of the $ \xi\eta $-plane onto the $ xy $-plane and the $ \zeta $-axis to the $ z $-axis, the three displacements can be expressed as

$$u=\sum_{i=1}^{N_{1}} \sum_{j=1}^{N_{2}} N_{i}(\xi, \eta) H_{j}(\zeta) u_{i j}$$

$$v=\sum_{i=1}^{N_{1}} \sum_{j=1}^{N_{2}} N_{i}(\xi, \eta) H_{j}(\zeta) v_{i j}$$

$$w=\sum_{i=1}^{N_{1}} N_{i}(\xi, \eta) w_{i}$$

in which $ N_{1} $ and $ N_{2} $ are the numbers of nodes in the $ xy $-plane (or $ \xi\eta $-plane) and $ z $-axis ($ \zeta $-axis), respectively. In addition, the first subscript for $ u $ and $ v $ denotes the node numbering in terms of the $ xy $-plane ($ \xi\eta $-plane) and the second subscript indicates the node numbering in terms of the $ z $-axis ($ \zeta $-axis).

In the following for further studying this plate type, we set up a case $ N_{1}=4 $ and $ N_{2}=2 $; that is, four-node quadrilateral shape functions are used for the $ xy $-plane interpolation and linear shape functions are employed for the $ z $-axis interpolation. Nodal displacements $ u_{i1} $ and $ v_{i1} $ are displacements on the bottom surface of the plate element, while displacements $ u_{i2} $ and $ v_{i2} $ are located on the top surface. As we can see in Eqs.(64),(65) through Eq.(66) there is no rotational degree of freedom for this kind of bending element. This is represented in the previous figure, for ease of understanding.

In this present formulation, we have the both bending strain energy and transverse shear strain energy present. We can express bending and transverse shear energy in the following manner (with the help of displacements):

$$\left\{{ϵ}_b\right\}=\left\{\begin{array}{l}{ϵ}_x\\ {ϵ}_y\\ {\gamma }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\mathrm{\partial }}{\mathrm{\partial }x}& 0& 0\\ 0& \frac{\mathrm{\partial }}{\mathrm{\partial }y}& 0\\ \frac{\mathrm{\partial }}{\mathrm{\partial }y}& \frac{\mathrm{\partial }}{\mathrm{\partial }x}& 0\end{array}\right]\left\{\begin{array}{l}u\\ v\\ w\end{array}\right\}$$

$$\left\{{ϵ}_s\right\}=\left\{\begin{array}{l}{\gamma }_{yz}\\ {\gamma }_{xz}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\mathrm{\partial }}{\mathrm{\partial }z}& 0& \frac{\mathrm{\partial }}{\mathrm{\partial }x}\\ 0& \frac{\mathrm{\partial }}{\mathrm{\partial }z}& \frac{\mathrm{\partial }}{\mathrm{\partial }y}\end{array}\right]\left\{\begin{array}{l}u\\ v\\ w\end{array}\right\}$$

in the case ${\epsilon_{b}} $ is the bending strains and ${\epsilon_{s}} $ is the transverse shear strains. The normal strain along the plate thickness ${\epsilon_{z}} $ is omitted here. Substitution of displacements, Eqs.(64), (65) through (66), into the kinematic equations, Eq.(67) and (68), with $ N_{1}=4 $ and $ N_{2}=2 $ expresses both bending and shear strains in the following way:

$$\left\{{ϵ}_b\right\}=\left\{\begin{array}{l}{ϵ}_x\\ {ϵ}_y\\ {\gamma }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\mathrm{\partial }}{\mathrm{\partial }x}& 0& 0\\ 0& \frac{\mathrm{\partial }}{\mathrm{\partial }y}& 0\\ \frac{\mathrm{\partial }}{\mathrm{\partial }y}& \frac{\mathrm{\partial }}{\mathrm{\partial }x}& 0\end{array}\right]\left\{\begin{array}{l}u\\ v\\ w\end{array}\right\}$$

where bending part is

$$\left[B_b\right]=\left[\left[B_{b1}\right]\left[B_{b2}\right]\left[B_{b3}\right]\left[B_{b4}\right]\right]$$

$$\left[B_{bi}\right]=\left[\begin{array}{ccccc}{H}_1\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }x}& 0& H_2\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }x}& 0& 0\\ 0& H_1\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }y}& 0& H_2\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }y}& 0\\ H_1\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }y}& H_1\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }x}& H_2\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }y}& H_2\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }x}& 0\end{array}\right]$$

$$\left\{d_i^e\right\}=\left\{\begin{array}{lllll}{u}_{i1}& v_{i1}& u_{i2}& v_{i2}& w_i\end{array}\right\}$$

$$\left\{{ϵ}_s\right\}=\left[B_s\right]\left\{d^e\right\}$$

where shear part is

$$[\begin{array}{lllll}{B}_s]& =[\left[B_{s1}\right]& \left[B_{s2}\right]& \left[B_{s3}\right]& \left[B_{s4}\right]]\end{array}$$

$$\left[B_{si}\right]=\left[\begin{array}{ccccc}{N}_i\frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }z}& 0& N_i\frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }z}& 0& \frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }z}\\ 0& N_i\frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }z}& 0& N_i\frac{\mathrm{\partial }{H}_2}{\mathrm{\partial }z}& \frac{\mathrm{\partial }{H}_1}{\mathrm{\partial }y}\end{array}\right]$$

The constitutive equation for the isotropic material is

$$\left\{{\sigma }_b\right\}=\left[D_b\right]\left\{{ϵ}_b\right\}$$

where

$$\left\{{\sigma }_b\right\}=\{\begin{array}{lll}{\sigma }_x& {\sigma }_y& {{\tau }_{xy}\}}^T\end{array}$$

$$\left[D_b\right]=\frac{E}{1-v^2}\left[\begin{array}{ccc}1& v& 0\\ v& 1& 0\\ 0& 0& \frac{1-v}{2}\end{array}\right]$$

for the bending components and

$$\left\{{\sigma }_s\right\}=\left[D_s\right]\left\{{ϵ}_s\right\}$$

where

$$\left\{{\sigma }_s\right\}=\{\begin{array}{ll}{\tau }_{yz}& {{\tau }_{xz}\}}^T\end{array}$$

$$D_s=\frac{E}{2(1+v)}\left[\begin{array}{ll}1& 0\\ 0& 1\end{array}\right]$$

where Eq.(83) is the material property matrix for the plane stress condition as usually assumed for the plate bending theory.

For a unidirectional fibrous composite, the material property matrices become

$$D_b=\left[\begin{array}{ccc}{D}_{11}& D_{12}& 0\\ D_{12}& D_{22}& 0\\ 0& 0& D_{33}\end{array}\right]$$

in which

$$D_{11}=\frac{E_{1}}{1-v_{12} v_{21}}$$

$$D_{12}=\frac{E_{1} v_{21}}{1-v_{12} v_{21}}$$

$$D_{22}=\frac{E_{2}}{1-v_{12} v_{21}}$$

$$D_{33}=G_{12}$$

and

$$\left[D_s\right]=\left[\begin{array}{cc}{G}_{11}& 0\\ 0& G_{12}\end{array}\right]$$

Here, $ 1 $ and $ 2 $ denote the longitudinal and transverse directions of the unidirectional composite. $ E $ is the elastic modulus, $ G_{ij} $ is the shear modulus of the $ i-j $ plane and $ v_{ij} $ is the Poisson’s ratio for strain in the $ j $-direction when stressed in the $ i $-direction. There are five independent material properties for Eq.(83) through (88) because of the reciprocal relation $ \frac{v_{12}}{E_{1}}=\frac{v_{21}}{E_{2}} $.

The total potential energy can be expressed as

$$\Pi=U-W$$

where the internal strain energy $ U $ consists of two parts

$$U=U_{b}+U_{s}$$

The bending strain energy $ U_{b} $ is

$$U_b=\frac{1}{2}{\int }_{\mathrm{\Omega }}{\left\{{\sigma }_b\right\}}^T\left\{{ϵ}_b\right\}d\mathrm{\Omega }$$

and the transverse shear strain energy is

$$U_s=\frac{1}{2}{\int }_{\mathrm{\Omega }}{\left\{{\sigma }_s\right\}}^T\left\{{ϵ}_s\right\}d\mathrm{\Omega }$$

where $ \Omega $ is the plate domain. After finite element discretization, substitution of the previous equations into Eqs.(91) and (92) gives

$$U_b=\sum _e\frac{1}{2}{\left\{d^e\right\}}^T{\int }_{{\mathrm{\Omega }}^e}{\left[B_b\right]}^T\left[D_b\right]\left[B_b\right]d\mathrm{\Omega }\left\{d^e\right\}$$

and

$$U_s=\sum _e\frac{1}{2}{\left\{d^e\right\}}^T{\int }_{{\mathrm{\Omega }}^e}{\left[B_s\right]}^T\left[D_s\right]\left[B_s\right]d\mathrm{\Omega }\left\{d^e\right\}$$

where summation is performed over the total number of finite elements and superscript $ e $ indicates each element. Kinematic matrices $ \left[B_{b}\right] $ and $ \left[B_{s}\right] $ are provided in Eqs.(70) and (75) while constitutive matrices $ \left[D_{b}\right] $ and $ \left[D_{s}\right] $ are given in Eq.(79) and (82) for the isotropic material, and in Eqs.(83) and (88) for the unidirectional composite.

The external work is written as

$$W=\{d{\}}^T\{F\}$$

in which $ {d} $ is the system nodal displacement vector and $ {F} $ is the system force vector. There is no rotational degree of freedom for the present element, the external moment is applied to the force vector as a couple applied on the top and bottom nodes of the plate element shown in the previous figure.

Invoking the stationary value of the total potential energy yields the finite element matrix equation. The element stiffness matrix can be expressed as

$$\left[K_{e}\right]=\int_{\Omega^{e}}\left[B_{b}\right]^{T}\left[D_{b}\right]\left[B_{b}\right] d \Omega+\int_{\Omega^{e}}\left[B_{s}\right]^{T}\left[D_{s}\right]\left[B_{s}\right] d \Omega$$

III.6 Shell made of in-plane and bending elements

Shells usually are made up of curved surfaces, but if we consider the radius of curvature infinitely large, then the geometry of the shell becomes a flat plate. When a shell is divided into a number of small finite elements, each element may be considered as a flat plate with a slightly different orientation is space. Thus, each element may be modeled as a plate bending element. However, there is the problem of different orientation in connecting plate elements, which cause in-plane deformation in the neighbouring elements.

As a result, a shell element can be obtained as a combination of a plate bending element and a plane stress element in the same way that a 2-D frame element is obtained from a beam bending element and a bar element. In the case of plate bending an element has a transverse defection and two bending rotations as degrees of freedom per node, while the plane stress element has two in-plane displacements per node. This gives us five degrees of freedom per node, three displacements and two rotations, the stiffness matrix of this kind of shell element is given as

$$\left[\begin{array}{cc}\left[K_b\right]& [0]\\ [0]& \left[K_m\right]\end{array}\right]\left\{\begin{array}{l}{d}_b\\ d_m\end{array}\right\}=\left\{\begin{array}{l}\left\{F_b\right\}\\ \left\{F_m\right\}\end{array}\right\}$$

where we have $ K $, $ d $ and $ F $ representing stiffness matrix, nodal displacement/rotation vector and nodal force/moment vector respectively. The matrices are made up of two parts, one for the plate bending and one for plate stretching. Subscripts $ b $ and $ m $ denote bending and membrane (stretching) deformations of the shell element.

When shell elements are oriented differently, for example the corners of a cube or box-type shell structure, a bending rotation in one element becomes a so-called drilling rotation in another element. Therefore, to assemble the element matrices and vectors into the system matrix and vector, the drilling rotational degree of freedom $ \theta_{z} $ should be included. Resulting from this, the degrees of freedom of the element matrix and vector shall increase by one (to include the drilling degree of freedom per node). This leads to equation (97) to be modified into

$$\left[\begin{array}{ccc}\left[K_b\right]& [0]& 0\\ \left[\begin{array}{c}0\end{array}\right]& \left[\begin{array}{c}{K}_m\end{array}\right]& 0\\ 0& 0& 0\end{array}\right]\left\{\begin{array}{c}\left\{d_b\right\}\\ \left\{d_m\right\}\\ {\theta }_z\end{array}\right\}=\left\{\begin{array}{c}\left\{F_b\right\}\\ \left\{F_m\right\}\\ 0\end{array}\right\}$$

The matrix and vectors in equation (98) are expressed in terms of a local coordinate system that has $ x $- and $ y $-axes along the midplane of each flat shell element and the $ z $-axis normal to the plane. To assemble those matrices and vectors into the system matrix and vectors, the nodal degrees of freedom in terms of each local axis must be transformed into the corresponding nodal degrees of freedom in terms of the global coordinate axes common to all elements. Having the transformation matrix called $ \left[T\right] $,

$$\left\{d^{\text{local}}\right\}=[T]\left\{d^{\text{global}}\right\}$$

Matrix $ \left[T\right] $ transforms the global degrees of freedom into the local degrees of freedom. It consists of direction cosines between the global and local coordinate systems.

At every node we have the relation between local and global degrees of freedom expressed as

$$\left\{\begin{array}{l}{u}^{\text{local}}\\ v^{\text{local}}\\ w^{\text{local}}\\ {\theta }_x^{\text{local}}\\ {\theta }_y^{\text{local}}\\ {\theta }_z^{\text{local}}\end{array}\right\}=\left[\begin{array}{cccccc}{l}_{11}& l_{11}& l_{11}& 0& 0& 0\\ l_{11}& l_{11}& l_{11}& 0& 0& 0\\ l_{11}& l_{11}& l_{11}& 0& 0& 0\\ 0& 0& 0& l_{11}& l_{11}& l_{11}\\ 0& 0& 0& l_{11}& l_{11}& l_{11}\\ 0& 0& 0& l_{11}& l_{11}& l_{11}\end{array}\right]\left\{\begin{array}{l}{u}^{\text{global}}\\ v^{\text{global}}\\ w^{\text{global}}\\ {\theta }_x^{\text{global}}\\ {\theta }_y^{\text{global}}\\ {\theta }_z^{\text{global}}\end{array}\right\}$$

where $ l_{ij} $ is the direction cosine between local axis $ x_{i} $ and the global axis $ X_{j} $. This relationship can be used for each node. Thus, the transformation matrix $ \left[T\right] $ for a four-node element becomes

$$[T]=\left[\begin{array}{cccc}\left[T_d\right]& 0& 0& 0\\ 0& \left[T_d\right]& 0& 0\\ 0& 0& \left[T_d\right]& 0\\ 0& 0& 0& \left[T_d\right]\end{array}\right]$$

where matrix $ \left[T_{d}\right] $ is that used in Eq. (100) having a size of $ 6×6 $.

Using the transformation matrix, the transformed stiffness matrix and load vector are given below:

$$\left[K^{g l o b a l}\right]=[T]^{T}\left[K^{\text {local}}\right][T]$$

$$\left\{F^{\text{global}}\right\}=[T{]}^T\left\{F^{\text{local}}\right\}$$