Extended rotation-free plate and beam elements with shear deformation effects

Latest revision as of 11:34, 4 January 2019

Published in Int. Journal for Numerical Methods in Engineering Vol. 83 (2), pp.196-227, 2010
doi: 10.1002/nme.2836

Abstract

The paper describes a methodology for extending rotation-free plate and beam elements in order to accounting for transverse shear deformation effects. The ingredients for the element formulation are: a Hu-Washizu type mixed functional, a linear interpolation for the deflection and the shear angles over standard finite elements and a finite volume approach for computing the bending moments and the curvatures over a patch of elements. As a first application of the general procedure we present an extension of the 3-noded rotation-free basic plate triangle (BPT) originally developed for thin plate analysis to accounting for shear deformation effects of relevance for thick plates and composite laminated plates. The nodal deflection degrees of freedom (DOFs) of the original BPT element are enhanced with the two shear deformation angles. This allows to computing the bending and shear deformation energies leading to a simple triangular plate element with 3 DOFs per node (termed BPT+ element). For the thin plate case the shear angles vanish and the element reproduces the good behaviour of the original thin BPT element. As a consequence the element is applicable to thick and thin plate situations without exhibiting shear locking effects. The numerical solution for the thick case can be found iteratively starting from the deflection values for the Kirchhoff theory using the original thin BPT element. A 2-noded rotation-free beam element termed CCB+ applicable to slender and thick beams is derived as a particular case of the plate formulation. The examples presented show the robustness and accuracy of the BPT+ and the CCB+ elements for thick and thin plate and beam problems.

Keywords: Rotation-free triangle, rotation-free beam, thick and thin plates and beams, finite elements, shear deformation

1 INTRODUCTION

In past years some authors have derived a number of thin plate and shell triangular elements free of rotational degrees of freedom (DOFs) based on Kirchhoff theory [1]–[24]. This idea goes back to the original attempts to solve thin plate bending problems using finite difference schemes with the deflection as the only nodal variable [25]–[27]. In essence all methods attempt to express the curvature field over an element in terms of the displacements of a collection of nodes belonging to a patch of adjacent elements. Oñate and Cervera [6] proposed a general procedure of this kind combining finite element and finite volume concepts for deriving thin plate triangles and quadrilaterals with the deflection as the only nodal variable and presented a simple and competitive rotation-free three DOFs triangular element termed BPT (for Basic Plate Triangle). These ideas were extended in [12] to derive a number of rotation-free thin plate and shell triangles. The basic ingredients of the method are a mixed Hu-Washizu formulation, a standard discretization into three-noded triangles, a linear finite element interpolation of the displacement field within each triangle and a finite volume type approach for computing constant curvature and bending moment fields within appropriate non-overlapping control domains. The so called cell-centered and cell-vertex triangular domains yield different families of rotation-free plate and shell triangles. Both the BPT plate element and its extension to shell analysis (termed BST for Basic Shell Triangle) can be derived from the cell-centered formulation. Here the control domain is an individual triangle. The constant curvature field within a triangle is computed in terms of the displacements of the six nodes belonging to the four element patch formed by the chosen triangle and the three adjacent triangles. In the cell-vertex approach the control volume is the tributary domain of a node and this leads to a different family of rotation-free plate and shell triangles. The detailed derivation of the different rotation-free triangular plate and shell element families can be found in [12]. The three-noded BST element has been successfully extended to non-linear shell problems involving frictional-contact situations and dynamics [15,19,20]. Practical applications of the BST element to sheet stamping analysis are reported in [9,10,11,24].

The paper describes an extension of the original rotation-free thin BPT element to accounting for transverse shear deformation effects of relevance for thick plates and composite laminated plates. The nodal deflection DOFs of the original BPT element are enhanced with the two shear deformation angles. This allows to computing the bending and shear deformation energies leading to a simple triangular plate element with 3 DOFs per node (termed BPT+ element). For the thin plate case the shear angles vanish and the element reproduces the good behaviour of the original thin BPT element. As a consequence the element is applicable to thick and thin plate situations without exhibiting shear locking. It is interesting that the thick plate solution can be found iteratively starting from the deflection values obtained using the standard Kirchhoff theory and the original thin BPT element. The ingredients of the formulation are: a Hu-Washizu type mixed functional, a linear interpolation for the deflection and the shear angles over 3-noded triangles and a finite volume approach for computing the bending moments and the curvatures over a patch of elements. Details of the element formulation are given in the paper.

The rotation-free formulation described for the BPT+ element is taken as the starting point for deriving a two-noded rotation-free beam element with shear deformation effects (termed CCB+ element). The examples presented in the last part of the paper show the robustness and accuracy of the BPT+ and the CCB+ elements for thick and thin plate and beam problems.

2 BASIC THEORY

2.1 Reissner-Mindlin plate theory

Let us consider the plate of Figure 1. We will assume Reissner-Mindlin conditions to hold, i.e.

(1)

with

(2.a)

or

(2.b)

with

(2.c)

where ${\textstyle u,v,w}$ are the cartesian displacements ${\textstyle \theta _{x}}$ , ${\textstyle \theta _{y}}$ are the rotations and ${\textstyle \phi _{x}}$ , ${\textstyle \phi _{y}}$ are angles coinciding with the transverse shear deformations as shown below.

Figure 1: Sign convenion for the deflection and the rotations in a plate

The generalized bending and shear strain vectors are defined as [29,30]

{\boldsymbol {\varepsilon }}_{b}=[\kappa _{x},\kappa _{y},\kappa _{xy}]^{T}={\begin{bmatrix}\displaystyle {\partial \theta _{x} \over \partial x},&\displaystyle {\partial \theta _{y} \over \partial y},&\displaystyle \left({\partial \theta _{x} \over \partial y}+{\partial \theta _{y} \over \partial x}\right)\end{bmatrix}}^{T}={\boldsymbol {L}}_{b}{\boldsymbol {\theta }}

(3.a)

(3.b)

where ${\textstyle [\kappa _{x},\kappa _{y},\kappa _{xy}]}$ and ${\textstyle [\gamma _{xz},\gamma _{yz}]}$ denote the bending strains (typically called curvatures) and the transverse shear strains, respectively and

{\boldsymbol {L}}_{b}={\begin{bmatrix}\displaystyle {\partial  \over \partial x}&0\\[.3cm]0&\displaystyle {\partial  \over \partial y}\\[.3cm]\displaystyle {\partial  \over \partial y}&\displaystyle {\partial  \over \partial x}\\\end{bmatrix}}\quad ,\quad {\boldsymbol {L}}_{s}={\begin{bmatrix}\displaystyle {\partial  \over \partial x}&-1&0\\[.3cm]\displaystyle {\partial  \over \partial y}&0&-1\\\end{bmatrix}}\quad ,\quad {\boldsymbol {u}}=\left\{{\begin{array}{c}w\\\theta _{x}\\\theta _{y}\end{array}}\right\}

(4)

Substituting Eq.(2.a) into (3.b) gives

(5)

i.e. the angles ${\textstyle \phi _{x}}$ , ${\textstyle \phi _{y}}$ coincide (with opposite sign) with the transverse shear deformations. In the following ${\textstyle \phi _{x}}$ and ${\textstyle \phi _{y}}$ will be called“shear angles”.

The set of governing equations can be expressed in integral form starting from the following Hu-Washizu type functional [12,28]

{\begin{array}{ll}\Pi =&\displaystyle {1 \over 2}\left[\iint _{A}{\boldsymbol {\varepsilon }}_{b}^{T}{\boldsymbol {D}}_{b}{\boldsymbol {\varepsilon }}_{b}\,dA+\iint _{A}{\boldsymbol {\varepsilon }}_{s}^{T}{\boldsymbol {\sigma }}_{s}\,dA\right]+\iint _{A}[{\boldsymbol {L}}_{b}{\boldsymbol {\theta }}-{\boldsymbol {\varepsilon }}_{b}]^{T}{\boldsymbol {\sigma }}_{b}\,dA-\\&-\iint _{A}(qw+{\boldsymbol {\theta }}^{T}{\boldsymbol {m}})dA-\sum \limits _{i=1}^{n_{c}}w_{i}P_{i}\end{array}}

(6)

where ${\textstyle q}$ is the distributed vertical loading, ${\textstyle {\boldsymbol {m}}=[m_{x},m_{y}]^{T}}$ are distributed bending moments, ${\textstyle P_{i}}$ is the vertical force acting at point ${\textstyle i}$ , ${\textstyle n_{c}}$ is the number of points with external concentrated forces and ${\textstyle A}$ is the area of the plate. Concentrated bending moments have been excluded for simplicity.

In Eq.(6) ${\textstyle {\boldsymbol {\sigma }}_{b}=[M_{x},M_{y},M_{xy}]^{T}}$ is the bending moment vector, ${\textstyle {\boldsymbol {\sigma }}_{s}=[Q_{x},Q_{y}]^{T}}$ is the shear force vector and ${\textstyle {\boldsymbol {D}}_{b}}$ is the bending constitutive matrix given by (for the isotropy case)

(7.a)

The shear forces and the shear strains are assumed to be related point-wise by the standard constitutive equation

(7.b)

In Eqs.(7) ${\textstyle t}$ is the plate thickness and ${\textstyle E,\nu }$ and ${\textstyle G}$ are the Young modulus, the Poisson's ratio and the shear modulus, respectively and ${\textstyle \alpha }$ is the shear correction factor; ${\textstyle \alpha =5/6}$ is taken for an isotropic plate.

2.2 Using the deflection and the shear angles as main variables

Let us express Eqs.(3)–(6) in terms of the deflection ${\textstyle w}$ and the shear angles ${\textstyle \phi _{x}}$ and ${\textstyle \phi _{y}}$ . The resulting Reissner-Mindlin theory contains the standard expressions of Kirchhoff thin plate theory plus additional transverse shear deformation terms.

Substituting Eq.(2) into Eqs.(3) gives

{\boldsymbol {\varepsilon }}_{b}={\boldsymbol {L}}_{b}{\boldsymbol {\theta }}=\left[{\frac {\partial ^{2}w}{\partial x^{2}}}+{\partial \phi _{x} \over \partial x},{\frac {\partial ^{2}w}{\partial y^{2}}}+{\partial \phi _{y} \over \partial y},2{\frac {\partial ^{2}w}{x\partial y}}+\left({\partial \phi _{x} \over \partial y}+{\partial \phi _{y} \over \partial x}\right)\right]^{T}={\boldsymbol {L}}_{w}w+{\boldsymbol {L}}_{b}{\boldsymbol {\phi }}

(8)

with

(9)

Substituting Eqs.(2.b) and (8) into the functional of Eq.(6) gives

{\begin{array}{ll}\Pi =&\!\!\!\!\displaystyle {1 \over 2}\left[\iint _{A}{\boldsymbol {\varepsilon }}_{b}^{T}{\boldsymbol {D}}_{b}{\boldsymbol {\varepsilon }}_{b}\,dA-\iint _{A}{\boldsymbol {\phi }}^{T}{\boldsymbol {\sigma }}_{s}\right]+\iint _{A}[{\boldsymbol {L}}_{w}w+{\boldsymbol {L}}_{b}{\boldsymbol {\phi }}-{\boldsymbol {\varepsilon }}_{b}]^{T}{\boldsymbol {\sigma }}_{b}\,dA+\\[.4cm]&\displaystyle -\!\!\iint _{A}[qw+({\boldsymbol {\nabla }}w+{\boldsymbol {\phi }})^{T}{\boldsymbol {m}}]dA-\!\sum \limits _{i=1}^{n_{c}}w_{i}P_{i}\end{array}}

(10)

Note that by making the shear angles ${\textstyle {\boldsymbol {\phi }}}$ equal to zero we recover precisely the variational form of standard Kirchhoff thin plate theory [29].

Variation of ${\textstyle \Pi }$ with respect to ${\textstyle {\boldsymbol {\varepsilon }}_{b}}$ , ${\textstyle {\boldsymbol {\sigma }}_{b}}$ , ${\textstyle w}$ and ${\textstyle {\boldsymbol {\phi }}}$ leads to the following equations:

Bending constitutive equation

(11)

Relationship between bending strains, deflection and shear angles

(12)

Equilibrium equations

$\iint _{A}[{\boldsymbol {L}}_{w}\delta w]^{T}{\boldsymbol {\sigma }}_{b}\,dA-\iint _{A}[\delta wq+({\boldsymbol {\nabla }}\delta w)^{T}{\boldsymbol {m}}]\,dA-\sum \limits _{i=1}^{n_{c}}\delta w_{i}P_{i}=0$	(13.a)
$\iint _{A}\left[({\boldsymbol {L}}_{b}\delta {\boldsymbol {\phi }})^{T}{\boldsymbol {\sigma }}_{b}-\delta {\boldsymbol {\phi }}^{T}{\boldsymbol {\sigma }}_{s}\right]dA-\iint _{A}\delta {\boldsymbol {\phi }}^{T}{\boldsymbol {m}}\,dA=0$	(13.b)

These integral expressions, completed with the constitutive equation for the shear forces (Eq.(7.b)), are the basis for deriving the element equations.

3 THREE-NODED BASIC ROTATION-FREE PLATE TRIANGLE WITH SHEAR DEFORMATION EFFECTS

Let us consider an arbitrary discretization of the plate into standard three-noded triangles. We assume a linear interpolation of the deflection ${\textstyle w}$ and the shear angles ${\textstyle {\boldsymbol {\phi }}}$ within each element in terms of the nodal values in the standard manner [28]

(14)

where ${\textstyle N_{i}}$ are the linear shape functions of the three-noded triangle [28,30], ${\textstyle {\bar {(\cdot )}}}$ denotes nodal variables and

(15.a)

{\bar {\boldsymbol {\phi }}}^{e}=\left\{{\begin{array}{c}{\bar {\boldsymbol {\phi }}}_{1}\\{\bar {\boldsymbol {\phi }}}_{2}\\{\bar {\boldsymbol {\phi }}}_{3}\\\end{array}}\right\}\quad ,\quad {\bar {\boldsymbol {\phi }}}_{i}=\left\{{\begin{array}{c}{\bar {\phi }}_{x_{i}}\\{\bar {\phi }}_{y_{i}}\end{array}}\right\}

(15.b)

with

{\begin{array}{c}N_{i}=\displaystyle {1 \over 2A^{e}}(a_{i}^{e}+b_{i}^{e}x+c_{i}^{e}y),\\[.3cm]a_{i}^{e}=x_{j}^{e}y_{k}^{e}-x_{k}^{e}y_{j}^{e}\,,\,b_{i}^{e}=y_{jk}\,,\,c_{i}^{e}=x_{kj}\quad {\hbox{with }}y_{jk}=y_{j}^{e}-y_{k}^{e}\,,\,x_{kj}=x_{k}^{e}-x_{j}^{e}\quad i,j,k=1,2,3\end{array}}

(15.c)

and ${\textstyle A^{e}}$ is the element area [28]. In the above expressions and in the following super-index ${\textstyle e}$ denotes element values.

Note that the introduction of the linear approximation for ${\textstyle w}$ into the variational form (13a) will give a vanishing of the virtual bending energy term involving second derivatives of ${\textstyle w}$ . This problem is overcome by “relaxing” further the weak form using an assumed constant value for the curvatures and bending moments over the triangle and computing the curvatures from the integral of the slopes along the element sides as explained next.

3.1 Assumed curvatures and bending moments field

We describe the bending moments ${\textstyle {\boldsymbol {\sigma }}_{b}}$ and the curvatures ${\textstyle {\boldsymbol {\varepsilon }}_{b}}$ (and their virtual values) by constant fields within the triangle, i.e.

(16)

where ${\textstyle (\cdot )^{e}}$ denotes constant values within the triangle.

We define the term “patch of triangular elements”. This is typically formed by four elements: a central triangle and the three adjacent triangles. Central triangles adjacent to a boundary lack the element adjacent to the boundary side (Figure 2).

Eqs.(16) are introduced into the governing integral equations (11)–(13) which are modified as follows.

Constitutive equation for the bending moments

(17)

Recalling that the virtual curvatures are arbitrary we obtain the following constitutive equations for the assumed (constant) bending moments for each element

(18)

where ${\textstyle {\boldsymbol {D}}_{b}^{e}}$ is the average bending constitutive matrix for the triangular element.

Relationship between the bending strains, the deflection and the shear angles

Substituting Eqs.(16) into Eq.(12) gives

(19)


(a)	(b)
Figure 2: Triangular element patches. (a) Standard four element patch; (b) Three element boundary patch. Numbers 1,2,3 in brackets denote local node numbers for the element

As the virtual bending moments are arbitrary and ${\textstyle {\boldsymbol {\varepsilon }}_{b}^{e}}$ is constant within each element we get

(20)

The term ${\textstyle {\boldsymbol {L}}_{w}w}$ in Eq.(20) is integrated by parts to give

(21)

where

(22)

and ${\textstyle n_{x},n_{y}}$ are the components to the outward unit normal ${\textstyle {\boldsymbol {n}}}$ to the boundary of the element ${\textstyle \Gamma ^{e}}$ (Figure 2). The transformation of the integral of the curvature field over the element domain into the integral of the deflection gradient along the element boundary is a distinct feature of the BPT formulation [6,12,15,19,20].

Eq.(21) defines the bending strains as the sum of the integral of the deflection gradient along the boundary of the element and the integral over the element of the term ${\textstyle {\boldsymbol {L}}_{b}{\boldsymbol {\phi }}}$ including the gradients of the shear deformation angles. This term is constant within each element for a linear interpolation of ${\textstyle {\boldsymbol {\phi }}}$ .

Shear forces-shear angles relationship

Making ${\textstyle {\boldsymbol {\varepsilon }}_{s}=-{\boldsymbol {\phi }}}$ (Eq.(8)) and substituting the approximation for ${\textstyle {\boldsymbol {\phi }}}$ of Eq.(14) into Eq.(7.b) gives

(23)

Equilibrium equations

Integrating by parts the first integral in Eq.(13.a) and recalling that the bending moments are constant within each element and the additive property of the element integrals in the FEM, allows us to write Eqs.(13) as

$\sum \limits _{e}\left(\int _{\Gamma ^{e}}{\boldsymbol {T}}({\boldsymbol {\nabla }}\delta w)^{T}\,d\Gamma \right){\boldsymbol {\sigma }}_{b}^{e}-\iint _{A}\left[\delta wq+\left({\boldsymbol {\nabla }}\delta w\right)^{T}{\boldsymbol {m}}\right]\,dA+\sum \limits _{i=1}^{n_{c}}\delta w_{i}P_{i}=0$	(24)
$\sum \limits _{e}\left\{\left(\iint _{A^{e}}({\boldsymbol {L}}_{b}{\boldsymbol {\phi }})^{T}dA\right){\boldsymbol {\sigma }}_{b}^{e}-\iint _{A^{e}}\delta {\boldsymbol {\phi }}^{T}{\boldsymbol {\sigma }}_{s}dA\right\}-\iint _{A}\delta {\boldsymbol {\phi }}^{T}{\boldsymbol {m}}\,dA=0$	(25)

The sum in above expressions extends over all the elements in the mesh.

3.2 Element matrices and vectors

Bending strain matrices

Substituting the linear interpolation for ${\textstyle w}$ and ${\textstyle {\boldsymbol {\phi }}}$ of Eqs.(14) into the integral expressions for the constant bending strain field of Eq.(21) gives

(26.a)

(26.b)

and, therefore

(27)

where ${\textstyle {\boldsymbol {B}}_{w}}$ and ${\textstyle {\boldsymbol {B}}_{\phi }}$ are generalized bending strain matrices and ${\textstyle {\bar {\boldsymbol {w}}}^{e}}$ and ${\textstyle {\bar {\boldsymbol {\phi }}}^{e}}$ are nodal deflection and nodal shear angles vectors for the element. The expression for ${\textstyle {\bar {\boldsymbol {\phi }}}^{e}}$ is given in Eq.(15.b). The expression for ${\textstyle {\bar {\boldsymbol {w}}}^{e}}$ is given below (Eq.(31)).

The computation of ${\textstyle {\boldsymbol {B}}_{\phi }}$ is straightforward from the second integral of Eq.(21) as

(28)

Recalling that ${\textstyle N_{i}}$ are linear shape functions gives (using Eq.(15.c))

{\boldsymbol {B}}_{\phi _{i}}={\begin{bmatrix}\displaystyle {\partial N_{i} \over \partial x}&0\\[.25cm]0&\displaystyle {\partial N_{i} \over \partial x}\\[.45cm]\displaystyle {\partial N_{i} \over \partial y}&\displaystyle {\partial N_{i} \over \partial x}\\\end{bmatrix}}={1 \over 2A^{e}}{\begin{bmatrix}b_{i}^{e}&0\\0&c_{i}^{e}\\c_{i}^{e}&b_{i}^{e}\\\end{bmatrix}}

(29)

The computation of ${\textstyle {\boldsymbol {B}}_{w}^{e}}$ in Eq.(26.a) requires integrating the product of the components of the normal vector and the deflection slopes along the element sides. For the 3-noded triangle this is not so straightforward as ${\textstyle {\boldsymbol {\nabla }}w}$ is discontinuous at the element sides for a linear approximation of ${\textstyle w}$ . A simple way to overcoming this problem is to compute the deflection gradients at the

{\boldsymbol {B}}_{w}={1 \over A^{e}}\left[{\begin{matrix}y_{ij}{\bar {b}}_{i}^{b}+y_{ki}{\bar {b}}_{i}^{d}&y_{ij}{\bar {b}}_{j}^{b}+y_{jk}{\bar {b}}_{j}^{c}&y_{jk}{\bar {b}}_{k}^{c}+y_{ki}{\bar {b}}_{k}^{d}\\-x_{ij}{\bar {c}}_{i}^{b}-x_{ki}{\bar {c}}_{i}^{d}&-x_{ij}{\bar {c}}_{j}^{b}-x_{jk}{\bar {c}}_{j}^{c}&-x_{jk}{\bar {c}}_{k}^{c}-x_{ki}{\bar {c}}_{k}^{d}\\[][y_{ij}{\bar {c}}_{i}^{b}-x_{ij}{\bar {b}}_{i}^{b}&[y_{ij}{\bar {c}}_{j}^{b}-x_{jk}{\bar {b}}_{j}^{b}&[y_{jk}{\bar {c}}_{k}^{c}-x_{jk}{\bar {b}}_{k}^{c}\\+y_{ki}{\bar {c}}_{i}^{d}-x_{ki}{\bar {b}}_{i}^{d}]&+y_{jk}{\bar {c}}_{j}^{c}-x_{jk}{\bar {b}}_{j}^{c}]&+y_{ki}{\bar {c}}_{k}^{d}-x_{ki}{\bar {b}}_{k}^{d}]\end{matrix}}\right.

\left.{\begin{matrix}y_{ij}{\bar {b}}_{l}^{b}&y_{jk}{\bar {b}}_{m}^{c}&y_{ki}{\bar {b}}_{n}^{d}\\-x_{ij}{\bar {c}}_{l}^{b}&-x_{jk}{\bar {c}}_{m}^{c}&-x_{ki}{\bar {c}}_{n}^{d}\\y_{ij}{\bar {c}}_{l}^{b}-x_{ij}{\bar {b}}_{l}^{b}&y_{jk}{\bar {c}}_{m}^{c}-x_{jk}{\bar {b}}_{m}^{c}&y_{ki}{\bar {c}}_{n}^{d}-x_{ki}{\bar {c}}_{n}^{d}\\\end{matrix}}\right]

element sides as the average value of the gradients contributed by the two triangles adjacent to the side [6,12]. This gives

{1 \over A^{e}}\int _{\Gamma ^{e}}{\boldsymbol {T}}{\boldsymbol {\nabla }}wd\Gamma ={1 \over A^{e}}\sum \limits _{j=1}^{3}{l_{j}^{e} \over 2}{\boldsymbol {T}}_{j}\left[\sum \limits _{i=1}^{3}{\boldsymbol {\nabla }}N_{i}^{e}w_{i}^{e}+\sum \limits _{i=1}^{3}{\boldsymbol {\nabla }}N_{i}^{p}w_{i}^{p}\right]_{j}={\boldsymbol {B}}_{w}{\bar {\boldsymbol {w}}}^{e}

(30)

with

(31)

The sum in Eq.(30) extends over the three sides of an element ${\textstyle e}$ , ${\textstyle {\boldsymbol {T}}_{j}}$ is the transformation matrix of Eq.(22) for side ${\textstyle j}$ , ${\textstyle l_{j}^{e}}$ are the lengths of the element sides and superindex ${\textstyle p}$ refers to each of the three triangles ${\textstyle a,b,c}$ adjacent to the central triangle ${\textstyle e}$ with ${\textstyle p=a,b,c}$ for ${\textstyle j=1,2,3}$ (Figure 2). Sides are assigned the number of the opposite node. Thus, side ${\textstyle l_{3}^{e}}$ is opposite to node 3, etc.

Note that ${\textstyle {\boldsymbol {B}}_{w}}$ is a ${\textstyle 3\times 6}$ matrix relating the three bending strains with the deflections at the six nodes of the patch of triangles linked to element ${\textstyle e}$ . The explicit form of matrix ${\textstyle {\boldsymbol {B}}_{w}}$ is given in Box I. This matrix coincides with the curvature matrix of the original BPT element, as presented in [12].

Substituting Eq.(27) into (18) gives the relationship between the bending moments and the nodal values of the deflection and the shear angles for the element as

(32)

Stiffness matrices and nodal force vectors

Substituting Eqs.(14), (23), (26) and (32) into the equilibrium equations (24) and (25) yields

(33.a)

\sum \limits _{e}(\delta {\bar {\boldsymbol {\phi }}}^{e})^{T}\left\{{\boldsymbol {B}}_{\phi }^{T}{\boldsymbol {D}}_{b}^{e}\left[{\boldsymbol {B}}_{w}{\bar {\boldsymbol {w}}}^{e}+{\boldsymbol {B}}_{\phi }{\bar {\boldsymbol {\phi }}}^{e}\right]A^{e}+\left(\iint _{A^{e}}{\boldsymbol {N}}_{\phi }^{T}{\boldsymbol {D}}_{s}{\boldsymbol {N}}_{\phi }dA\right){\bar {\boldsymbol {\phi }}}^{e}\right\}-

(33.b)

After simplification of the virtual nodal deflections and the virtual nodal shear angles we obtain finally the matrix system of equilibrium equations

(34)

where vectors ${\textstyle {\bar {\boldsymbol {w}}}}$ and ${\textstyle {\bar {\boldsymbol {\phi }}}}$ contain the nodal deflections and the nodal shear angles for the whole mesh and the rest of the matrices and vectors are assembled from the element contributions given by

{\begin{array}{l}{\boldsymbol {K}}_{w_{ij}}^{e}=\displaystyle {\boldsymbol {B}}_{w_{i}}^{T}{\boldsymbol {D}}_{b}^{e}{\boldsymbol {B}}_{w_{j}}A^{e}\quad ,\quad {\boldsymbol {K}}_{w\phi _{ij}}^{e}={\boldsymbol {B}}_{w_{i}}^{T}{\boldsymbol {D}}_{b}^{e}{\boldsymbol {B}}_{\phi _{j}}A^{e}\\[.4cm]{\boldsymbol {K}}_{\phi _{ij}}^{e}=\displaystyle {\boldsymbol {B}}_{\phi _{i}}^{T}{\boldsymbol {D}}_{b}^{e}{\boldsymbol {B}}_{\phi _{j}}A^{e}+\iint _{A^{e}}{\boldsymbol {N}}_{i}^{T}{\boldsymbol {D}}_{s}{\boldsymbol {N}}_{j}\,dA\end{array}}

(35)

The element stiffness matrices ${\textstyle {\boldsymbol {K}}_{w}^{e}}$ and ${\textstyle {\boldsymbol {K}}_{w\phi }^{e}}$ can be explicitly computed from the strain and constitutive matrices for the element. The exact computation of the integral in the expression of ${\textstyle {\boldsymbol {K}}_{\phi }}$ requires a 3 Gauss point quadrature. Excellent results have been obtained in all examples analyzed using a reduced one point integration rule for ${\textstyle {\boldsymbol {K}}_{\phi }}$ which indicates no advantage in using a full quadrature for ${\textstyle \mathbf {K} _{\phi }}$ . The simple one point quadrature allows computing all the element stiffness matrices in Eq.(34) explicitly.

The equivalent nodal force vectors for the element are

(36)

For a uniform distribution of ${\textstyle q}$ and ${\textstyle {\boldsymbol {m}}}$

{\begin{array}{l}\displaystyle {f}_{w_{i}}^{e}={qA^{e} \over 3}+{1 \over 2}(b_{i}^{e}m_{x}+c_{i}^{e}m_{y})\\[.4cm]\displaystyle {\boldsymbol {f}}_{\phi _{i}}^{e}={A^{e} \over 3}[m_{x},m_{y}]^{T}\end{array}}

(37)

where ${\textstyle b_{i}^{e}}$ and ${\textstyle c_{i}^{e}}$ are given in Eq.(15).

The vertical load ${\textstyle P_{i}}$ acting at a node ${\textstyle i}$ contributes the following terms to the ${\textstyle i}$ th component of the global force vector ${\textstyle {\boldsymbol {f}}_{w}}$

(38)

The reactions (the vertical force and the bending moment) can be computed at the prescribed nodes from the nodal displacement and rotations field in the usual manner.

Remark 1 . Each node has three DOFs: the deflection ${\textstyle {\bar {w}}_{i}}$ and the two shear angles ${\textstyle {\bar {\phi }}_{x_{i}}}$ and ${\textstyle {\bar {\phi }}_{y_{i}}}$ at the node. The size of all the element stiffness matrices in Eqs.(35) is ${\textstyle 6\times 6}$ . This means that the effective DOFs for each element are in fact the nodal deflections of the patch of four elements assigned to each element (typically six nodal deflections except for boundary elements) and the six nodal shear angles, i.e. two shear angles for each of the three nodes of the element.

3.3 Iterative computation of the deflection and the shear angles at the nodes

Despite that the solution for the ${\textstyle {\bar {\boldsymbol {w}}}}$ and ${\textstyle {\bar {\boldsymbol {\phi }}}}$ variables can be found simultaneously by solving Eqs.(34), the following iterative algorithm is recommended for computing ${\textstyle {\bar {\boldsymbol {w}}}}$ and ${\textstyle {\bar {\boldsymbol {\phi }}}}$

$\displaystyle {\boldsymbol {K}}_{w}{\bar {\boldsymbol {w}}}^{i}\!\!\!=\!\!\!{\boldsymbol {f}}_{w}-{\boldsymbol {K}}_{w\phi }{\bar {\boldsymbol {\phi }}}^{i-1}\rightarrow {\bar {\boldsymbol {w}}}^{i}$	(39.a)
$\displaystyle {\boldsymbol {K}}_{\phi }{\bar {\boldsymbol {\phi }}}^{i}\!\!\!=\!\!\!{\boldsymbol {f}}_{\phi }-{\boldsymbol {K}}_{w\phi }^{T}{\bar {\boldsymbol {w}}}^{i}\rightarrow {\bar {\boldsymbol {\phi }}}^{i}$	(39.b)

where superindex ${\textstyle i}$ denotes the number of iterations. The iterative scheme of Eqs.(39) continues until convergence for the nodal deflection and the nodal shear deformation angles is achieved. Convergence is typically measured by the ${\textstyle L_{2}}$ norm of vectors ${\textstyle {\bar {\boldsymbol {w}}}}$ and ${\textstyle {\bar {\boldsymbol {\phi }}}}$ (Eqs.(58) and (59)). An advantage of the above iterative scheme is that for ${\textstyle i=1}$ and ${\textstyle {\boldsymbol {\phi }}^{0}={\boldsymbol {0}}}$ the value of ${\textstyle {\bar {\boldsymbol {w}}}^{1}}$ corresponds to the thin plate solution of Kirchhoff theory which is accurate enough for many practical cases. The effect of shear deformation is introduced progressively with the number of iterations. Shear deformation effects are negligible for thin plates and hence the ${\textstyle {\bar {\boldsymbol {\phi }}}}$ variables tend rapidly to zero in this case.

3.4 Improved iterative scheme

An enhanced iterative scheme can be devised by performing a smoothing of the shear angles field as follows.

Step 1 Computation of the nodal deflections ̄w¹

(40)

Step 2 Compute ̄ϕⁱ, i ≥ 1

(41)

Step 3 Compute the element shear angles ̄ϕ^e

(42)

where ${\textstyle {\bar {\boldsymbol {\phi }}}^{e}}$ are mean shear angles for element ${\textstyle e}$ for the ${\textstyle i}$ th iteration and ${\textstyle a,b,c}$ are the global numbers of the three nodes of element ${\textstyle e}$ .

Step 4 Compute the smoothed nodal shear angles ̃ϕₖ

The smoothed nodal values ${\textstyle {\tilde {\boldsymbol {\phi }}}_{k}^{i}}$ are computed as

(43)

where ${\textstyle n_{k}}$ is the number of elements sharing the node with global number ${\textstyle k}$ and ${\textstyle N}$ is the total number of nodes in the mesh.

Step 5 Compute ̄wⁱ, i > 1

(44)

Return to step 2.

Convergence of the improved iterative scheme for the deflection field is quite fast (2–4 iterations), even for thick plates (see Section 6).

The accuracy and efficiency of the so called BPT+ triangle is shown in the examples presented in the paper.

4 BOUNDARY CONDITIONS

A BPT element with a side along a boundary edge has one of the triangles belonging to the patch missing (Figure 2). This is taken into account by ignoring the contribution of this element when performing the average of the deflection gradient in Eq.(30) [12,15,19,20].

4.1 Clamped edge (w=0, θ=0)

The condition ${\textstyle w=0}$ is directly imposed at the edge nodes at the solution level when solving the system of equations for ${\textstyle {\bar {\boldsymbol {w}}}}$ in the standard manner.

The condition of zero rotations at a clamped edge introduces additional terms in the system of Eqs.(34).

In order to explain the process, let us consider for example a clamped edge corresponding to side 3 linking nodes 1 and 2 (with global numbers ${\textstyle i}$ and ${\textstyle j}$ ) in the boundary element of Figure 2b. At the clamped edge

(45)

where ${\textstyle {\boldsymbol {\theta }}\vert _{3}^{e}}$ are the rotations at the edge mid point.

The shear angles at the edge mid point are approximated by the average of the nodal values along the edge, i.e.

(46)

Introducing Eqs.(45) and (46) into the definition of matrix ${\textstyle {\boldsymbol {B}}_{w}}$ in Eq.(30) gives

{1 \over A^{e}}\int _{\Gamma ^{e}}\!\!{\boldsymbol {\nabla }}wd\Gamma \!\!=\!\!{1 \over A^{e}}\sum \limits _{j=1}^{2}{l_{j}^{e} \over 2}{\boldsymbol {T}}_{j}\left[\sum \limits _{i=1}^{3}{\boldsymbol {\nabla }}N_{i}^{e}w_{i}^{e}+\sum \limits _{i=1}^{3}{\boldsymbol {\nabla }}N_{i}^{p}w_{i}^{p}\right]_{j}-{1 \over A^{e}}{l_{3}^{e} \over 2}{\boldsymbol {T}}_{3}[{\boldsymbol {\phi }}_{2}+{\boldsymbol {\phi }}_{3}]

(47)

where ${\textstyle {\boldsymbol {B}}_{w}}$ is obtained by disregarding the contributions from the clamped side in the sum along the element sides in the expression of Box I and

(48)

The bending strain field of Eq.(27) is now modified as

(49)

where

(50)

The new matrix ${\textstyle {\bar {\boldsymbol {B}}}_{\phi }}$ substitutes matrix ${\textstyle {\boldsymbol {B}}_{\phi }}$ in the expression for ${\textstyle {\boldsymbol {K}}_{\phi }^{e}}$ of Eq.(35).

We recall that the above modifications are only needed for elements with a clamped edge.

The process is repeated twice if the element has two clamped edges.

4.2 Simply supported edge (w=0, θₛ=0)

The condition ${\textstyle w=0}$ at the nodes laying on a simply supported (SS) edge is prescribed when solving the global system of equations as for the clamped case. Prescribing ${\textstyle w=0}$ at an edge node also automatically implies that ${\textstyle {\partial w \over \partial s}}$ is zero along the SS edge direction ${\textstyle s}$ and this is the option taken for the standard rotation-free thin BPT element. Note that this is equivalent to assuming a “soft” simply support condition in general plate theory (i.e. ${\textstyle w=M_{n}=M_{ns}=0}$ at the SS edge) [29,30].

The “hard” support condition requires prescribing the tangential rotation ${\textstyle \theta _{s}}$ to a zero value. This means

(51)

Prescribing ${\textstyle w}$ at the SS nodes gives automatically ${\textstyle \left({\frac {\partial w}{\partial s}}\right)_{i}=0}$ and hence the nodal variable ${\textstyle {\bar {\phi }}_{s_{i}}}$ must be also prescribed to a zero value at these nodes. This condition can be imposed by transforming the cartesian shear strains to the boundary axes ${\textstyle s,n}$ and making ${\textstyle {\bar {\phi }}_{s_{i}}=0}$ at each simply supported boundary node.

Indeed, prescribing just ${\textstyle w_{i}=0}$ at the support nodes and letting ${\textstyle {\bar {\phi }}_{s_{i}}}$ free reproduces the soft support condition for Reissner-Mindlin theory. This is the approach chosen in the examples presented in the paper ( Figures 5, 6, 9, 10).

It is interesting that the BPT element does not suffer from the difficulties associated to prescribing the tangential rotation at SS nodes which occur in standard thin plate elements [29,32,33]. This is due to the fact that the rotations do not appear explicitly as variables in the BPT formulation and the SS condition is directly imposed by prescribing the deflection at the support nodes.

4.3 Symmetry edge (θₙ=0)

The condition of zero normal rotation ( ${\textstyle \theta _{n}=0}$ ) is imposed by neglecting the contributions from the normal rotation ${\textstyle {\frac {\partial w}{\partial n}}}$ at the symmetry edge when computing Eq.(30). The condition ${\textstyle \theta _{n}={\frac {\partial w}{\partial n}}+\phi _{n}=0}$ for the thick case introduces additional terms in the stiffness equations. The procedure is identical as explained for the clamped edge.

Remark 2 . The condition ${\textstyle {\bar {\phi }}_{n_{i}}=0}$ can also be imposed at free edge nodes. This however has not shown to be improve the quality of the solution or to lead to any computational advantage, other that the reduction in the number of DOFs, in those cases.

5 TWO-NODED ROTATION-FREE BEAM ELEMENT

A two-noded rotation-free beam element with shear deformation effects can be simply derived as a particular case of the formulation for the BPT element previously described. The resulting beam element is termed CCB+ as a reference to the cell-centered approach used to compute the curvature at the element center.

Figure 3 shows the patch of three beam elements needed for computing the constant curvature at the central element using a cell-centered finite-volume type scheme. The starting point in the formulation of the CCB+ element is the standard Timoshenko beam theory [29,30]. The relevant expressions are:

{\begin{array}{l}{\boldsymbol {Deflection}}:\,w\\[.2cm]{\boldsymbol {Rotation}}:\,\theta =\displaystyle {\frac {\partial w}{\partial x}}+\phi \\[.2cm]{\boldsymbol {Displacement~vector}}:\,{\boldsymbol {u}}=[w,\theta ]^{T}\\[.2cm]{\boldsymbol {Curvature}}:\,{\boldsymbol {\varepsilon }}_{b}=[{\kappa }]=\left[\displaystyle {\frac {\partial \theta }{\partial x}}\right]={\boldsymbol {L}}_{b}\theta \quad ,\quad {\boldsymbol {L}}_{b}=\left[\displaystyle {\frac {\partial }{\partial x}}\right]\\[.2cm]{\boldsymbol {Transverse~shear~deformation}}:\,{\boldsymbol {\sigma }}_{s}=[\gamma ]=\left[\displaystyle {\frac {\partial w}{\partial x}}-\theta \right]={\boldsymbol {L}}_{s}{\boldsymbol {u}}\quad ,\quad {\boldsymbol {L}}_{s}=\left[\displaystyle {\frac {\partial }{\partial x}},-1\right]\qquad \qquad \qquad \\\end{array}}

(52)

${\boldsymbol {Constitutive~equations:}}$

(53)

where ${\textstyle I}$ and ${\textstyle A}$ are the inertia modulus and the area of the transverse cross section of the beam and the rest of the terms have been defined previously.

The Hu-Washizu functional for the beam has the same form as in Eq.(10) with ${\textstyle {\boldsymbol {m}}=[m]}$ .

Similarly, the governing equations have identical expressions as for the plate problem Eqs.(11–13).

A linear interpolation is chosen for the deflection ${\textstyle w}$ and the shear angle ${\textstyle \phi }$ as

(54)

Figure 3: (a) Patch of three rotation-free CCB+ elements. (b) Patch of element adjacents to a clamped boundary

{\begin{array}{l}\displaystyle \mathbf {B} _{w}={\frac {1}{2l^{a}l^{e}l^{b}}}\left[l^{b},-l^{b},-l^{a},l^{a}\right]\\[.3cm]\displaystyle \mathbf {B} _{\phi }={\frac {1}{l^{e}}}[1,-1]\end{array}}

Box II. Matrices

and

for the 2-noded rotation-free CCB+ element

where ${\textstyle N_{i}=1/2(1+\xi \xi _{i})}$ with ${\textstyle \xi _{1}=-1}$ and ${\textstyle \xi _{2}=1}$ are the standard linear shape functions for the 2-noded Lagrange element [28,30].

Following an identical process as for the 3-noded BPT+ element, the same system of Eqs.(34) is found. The expressions for the stiffness matrices ${\textstyle {\boldsymbol {K}}_{w}^{e}}$ , ${\textstyle {\boldsymbol {K}}_{\phi }^{e}}$ and ${\textstyle {\boldsymbol {K}}_{w_{\phi }}^{e}}$ coincide with those given in Eq.(35) simply by noting that the integrals change from area to line ones and substituting the area ${\textstyle A^{e}}$ by the element length ${\textstyle l^{e}}$ . The form of matrices ${\textstyle {\boldsymbol {B}}_{w}}$ and ${\textstyle {\boldsymbol {B}}_{\phi }}$ is given in Box II.

The expressions for the equivalent force vector are similar to Eqs.(36)–(38), substituting the area integrals by integrals along the element length. For a uniform distribution of ${\textstyle q}$ and ${\textstyle m}$

(55)

where ${\textstyle \xi _{i}}$ is the value of the natural coordinate ${\textstyle \xi }$ at node ${\textstyle i}$ .

The procedure for imposing the boundary conditions follows precisely the lines described for the BPT+ element in Section 4. The ${\textstyle \Delta {\boldsymbol {B}}_{\phi }}$ matrix for an element with a clamped node is (Eq.(48) and Figure 3b)

(56)

(57)

A similar expression is obtained for an element with a node on the symmetry axis.

The CCB+ element is an extension of the 2-noded rotation-free CCB element based on Euler-Bernouilli beam theory. A description of the CCB element can be found in 11,30.

6 EXAMPLES

6.1 Square and circular plates

The efficiency and accuracy of the BPT+ element has been tested in the analysis of a number of square plates of side ${\textstyle L}$ and circular plates of diameter ${\textstyle 2L}$ under a uniformly distributed loading and a central point load for different thicknesses ranging from to ${\textstyle t/L=10^{-3}}$ (very thin plate) to ${\textstyle t/L=0.1}$ (thick plate). Simply supported and clamped boundary conditions have been considered. For rectangular plates the soft SS condition has been assumed. For circular plates ${\textstyle \phi _{s}}$ is automatically zero at the SS boundary for the problems studied due to symmetry and hence results correspond to the hard SS case. Figure 4 shows some of the different meshes used for the analysis.

Results for the cases studied using the enhanced iterative scheme of Section 3.4 are presented in Figures 5–12. Each figure shows:

The convergence of the vertical deflection values and the shear angles with the number of iterations measured as

(58)

L_{2}^{\boldsymbol {\phi }}=\left[\sum \limits _{j=1}^{N}{\frac {[{\bar {\boldsymbol {\phi }}}_{j}^{i}-{\bar {\boldsymbol {\phi }}}_{j}^{i-1}]^{T}({\bar {\boldsymbol {\phi }}}_{j}^{i}-{\bar {\boldsymbol {\phi }}}_{j}^{i-1})}{[{\boldsymbol {\phi }}_{j}^{i}]^{T}{\boldsymbol {\phi }}_{j}^{i}}}\right]^{1/2}

(59)

where ${\textstyle N}$ is the number of nodes in the mesh and an upper index denotes the iteration number. A value of ${\textstyle {\bar {w}}_{j}^{0}=0}$ and ${\textstyle {\bar {\boldsymbol {\phi }}}_{j}^{0}=\mathbf {0} }$ has been taken. For the examples considered, the iterative scheme stops when ${\textstyle L_{2}^{w}<10^{-3}}$ .

The value of the normalized central deflection ( ${\textstyle w_{c}={\frac {D{\bar {w}}_{c}\cdot 10^{5}}{qL^{4}}}}$ for uniform load and ${\textstyle w_{c}={\frac {D{\bar {w}}_{c}\cdot 10^{4}}{PL^{2}}}}$ for point load with ${\textstyle D={\frac {Et^{3}}{12(1-\nu ^{2})}}}$ ) for the range of thickness ratios ${\textstyle 10^{-3}-10^{-1}}$ for each of the meshes considered. Results are compared with analytical and series values for the thin and thick cases [26,29,30,31] when available, or, alternatively, with FEM results obtained for the deflection at the center of the midle plane using a mesh of ${\textstyle 40\times 40\times 6}$ eight-noded hexahedra in a quarter of plate.
For the point load case the analytical value for the deflection under the load given by thick plate theory is infinity. Hence, results for the deflection for thick plates are compared at the mid-point along a central line in this case.

The distribution of the bending moment ${\textstyle M_{x}}$ and the shear force ${\textstyle Q_{x}}$ along the central line for the thick case ( ${\textstyle t/L=0.10}$ ) for each of the five meshes considered. The isovalues of ${\textstyle M_{x}}$ and ${\textstyle Q_{x}}$ over a quarter of the plate are also shown for the finer mesh.

The following conclusions are drawn from the examples:

Figure 4: Some of meshes of BPT+ elements used for analysis of square and circular plates

Figure 5: Simple supported square plate (soft support) under uniform load. Convergence of central deflection for different thicknesses. Upper curves show convergence of the vertical deflection and the shear angles for a thick plate with the number of iterations. Lower curves show the distribution of

and