8.7 Structural Elements: Two force members, Beams and Plates

 

It would take many millions of three-dimensional continuum elements to mesh a structure such as a bridge or building, even though the stress and strain fields within individual structural members have a simple form.   So, to make computations easier and faster, special ‘structural elements’ are nearly always used to model a complete structure of this kind.   These elements all approximate the internal stress and strain fields in the member in some appropriate way.   In this section, we discuss a few commonly used structural elements to illustrate their features.

 

 

 

8.7.1 Truss Elements.

 

A truss consists of a set of bars with a uniform cross-section, which are connected by pin joints.   The structure may be loaded by applying forces to some of the joints, while other joints may be fixed or subjected to a known displacement.   If the truss is in static equilibrium, its members are all in a state of uniaxial tension or compression.  Each member can be regarded as a single finite element (a ‘truss’ element), and the joints can be regarded as the nodes connecting the elements.  The joint displacements are the unknown degrees of freedom.

 

The figure shows a simple example of a truss. The geometry of the structure is characterized by

 

1. The initial and final positions  $x^a$, $y^a$ $a=\mathrm{1..}N$ of the N joints;

 

2. The joint displacements $u^a=y^a-x^a$;

 

3. The initial and deformed cross-sectional area $A_0^k,A^k$ $k=\mathrm{1..}m$ of the m members;

 

4. The initial and final length $L_0^k,L^k$ of the members;

 

5. The axial and transverse stretch in each member (with joints b and a at its two ends)

$\lambda ={\displaystyle \frac{L}{L_0}=\sqrt{{\displaystyle \frac{\left(y^b-y^a\right)\cdot \left(y^b-y^a\right)}{\left(x^b-x^a\right)\cdot \left(x^b-x^a\right)}}}}$

 

6. The axial component of infinitesimal strain in each member

$\epsilon ={\displaystyle \frac{\left(u^b-u^a\right)\cdot \left(x^b-x^a\right)}{\left(x^b-x^a\right)\cdot \left(x^b-x^a\right)}}$

 

7. The axial component of stretch rate in the member with joints a and b at its ends is given by

$D={\displaystyle \frac{1}{L^2}{\displaystyle \frac{d\left(y^b-y^a\right)}{dt}\cdot \left(y^b-y^a\right)}}$

 

A set of external point forces $P^n$ act at the joints.  The internal forces can be characterized by the axial component of Cauchy stress $\sigma$ in the members.

 

As usual, the finite element method solves the static equilibrium equation by re-writing it as the equivalent principle of virtual work.  For this purpose, we introduce a kinematically admissible virtual velocity $v^a$ for each joint (any displacement component that is prescribed for a joint must have a zero virtual velocity).  The virtual velocity has a corresponding virtual stretch rate

$d^k={\displaystyle \frac{\left(v^b-v^a\right)\cdot \left(y^b-y^a\right)}{\left(y^b-y^a\right)\cdot \left(y^b-y^a\right)}}$

The principal of virtual work is then

${\displaystyle \sum _{k=1}^m{A}^k{L}^k\sigma d^k}={\displaystyle \sum _{j=1}^n{P}^j\cdot v^j}$

where the sums are taken over the members and loaded joints, respectively.   

 

 

Elastic structures with small displacements:  Most practical structures are subjected to modest loads, so their deformations are small and the members remain elastic.   Under these conditions the strain in the members can be approximated by the axial component of the infinitesimal strain $\epsilon$, while the Cauchy stress is given by $\sigma =E\epsilon$, where E is the Young’s modulus of the member.   In addition, changes in the the cross-sectional area and length of the members can be neglected.  The infinitesimal strain in the kth member can then be calculated from the displacements of joints a and b at its two ends using the matrix expression

$\epsilon ={\displaystyle \frac{1}{{\left(L_0^k\right)}^2}\left[\begin{array}{cc}{x}^a-x^b& x^b-x^a\end{array}\right]\left[\begin{array}{c}{u}^a\\ u^b\end{array}\right]}$

Here, the vectors are 6 dimensional (for a 3D structure), with the bold symbols representing either three position vector components or three displacement components.  You can multiply out the vector product to see that this form is equivalent to the dot product in item (6) of the preceding sub-section.  The rate of change of infinitesimal strain can replace the virtual stretch rate.   Consequently, the principle of virtual work can be written as

${\displaystyle \sum _{k=1}^m\left[\begin{array}{cc}{v}^a& v^b\end{array}\right]\left[\begin{array}{c}{x}^a-x^b\\ x^b-x^a\end{array}\right]{\displaystyle \frac{E{A}^k}{{\left(L_0^k\right)}^3}\left[\begin{array}{cc}{x}^a-x^b& x^b-x^a\end{array}\right]}}\left[\begin{array}{c}{u}^a\\ u^b\end{array}\right]={\displaystyle \sum _{j=1}^n{P}^j\cdot v^j}$

The virtual work equation must be satisfied for all admissible joint velocities.  By choosing a set of virtual displacements in which only one displacement component at one joint is nonzero in turn for all the joints yields a system of linear equations

$Ku=PK={\displaystyle \sum _{k=1}^m{\displaystyle \frac{E{A}^k}{{\left(L_0^k\right)}^3}\left[\begin{array}{c}{x}^a-x^b\\ x^b-x^a\end{array}\right]\otimes \left[\begin{array}{cc}{x}^a-x^b& x^b-x^a\end{array}\right]}}$

where the symbol $\otimes$ denotes an ‘outer’ or dyadic product of the two vectors, P represents a column vector containing all the external forces, and the sum in the expression for the stiffness represents the assembly of element stiffness matrices into a global stiffness. After modifying the equations to account for any prescribed displacements, they may be solved to determine the unknown joint displacements.

We can illustrate the procedure in more detail using the simple 3-noded planar truss shown in the figure.   The joints have coordinates

$x^1=(0,0)x^2=(H,0)x^3=(0,H)$

The two element stiffnesses are therefore

$\begin{array}{l}{k}^1={\displaystyle \frac{EA}{H^3}\left[\begin{array}{c}-H\\ 0\\ H\\ 0\end{array}\right]\otimes \left[\begin{array}{cccc}-H& 0& H& 0\end{array}\right]={\displaystyle \frac{EA}{H}\left[\begin{array}{cccc}1& 0& -1& 0\\ 0& 0& 0& 0\\ -1& 0& 1& 0\\ 0& 0& 0& 0\end{array}\right]}}\\ k^2={\displaystyle \frac{EA}{{\left(\sqrt{2}H\right)}^3}\left[\begin{array}{c}H\\ -H\\ -H\\ H\end{array}\right]\otimes \left[\begin{array}{cccc}H& -H& -H& H\end{array}\right]={\displaystyle \frac{EA}{2\sqrt{2}H}\left[\begin{array}{cccc}1& -1& -1& 1\\ -1& 1& 1& -1\\ -1& 1& 1& -1\\ 1& -1& -1& 1\end{array}\right]}}\end{array}$

The element stiffnesses may be combined into a 6x6 global stiffness to yield the (singular) equation system

${\displaystyle \frac{EA}{2\sqrt{2}H}\left[\begin{array}{cccccc}2\sqrt{2}& 0& -2\sqrt{2}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ -2\sqrt{2}& 0& (1+2\sqrt{2})& -1& -1& 1\\ 0& 0& -1& 1& 1& -1\\ 0& 0& -1& 1& 1& -1\\ 0& 0& 1& -1& -1& 1\end{array}\right]\left[\begin{array}{c}{u}_1^1\\ u_2^1\\ u_1^2\\ u_2^2\\ u_1^3\\ u_2^3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ P_1\\ P_2\\ 0\\ 0\end{array}\right]}$

Finally, the equation system must be modified to ensure that the displacements at joints 1 and 3 are zero.   Usually this is done so that the global stiffness remains symmetric, in which case the final system of equations is

${\displaystyle \frac{EA}{2\sqrt{2}H}\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& (1+2\sqrt{2})& -1& 0& 0\\ 0& 0& -1& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{u}_1^1\\ u_2^1\\ u_1^2\\ u_2^2\\ u_1^3\\ u_2^3\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ P_1\\ P_2\\ 0\\ 0\end{array}\right]}$

 The system can be solved for the unknown displacement at node 2 with the result

$\left[\begin{array}{c}{u}_1^2\\ u_2^2\end{array}\right]={\displaystyle \frac{H}{EA}\left[\begin{array}{cc}1& 1\\ 1& 1+2\sqrt{2}\end{array}\right]\left[\begin{array}{c}{P}_1\\ P_2\end{array}\right]}$

 

 

Trusses made from nonlinear elastic materials that experience large displacements: More generally, the members of the truss may be made from a rubbery material that can be stretched significantly.   It is more difficult to calculate the deformation and internal forces in this kind of structure, because the stress-strain relation and the equilibrium equation are nonlinear.

 

As an example, consider a truss that has its members made from an incompressible nonlinear elastic material with uniaxial Cauchy stress-v-stretch relation $\sigma (\lambda )$  Some specific stress-stretch relations for rubber elasticity models are listed in Section 3.5.5: for example, the ‘neo-hookean’ material has a Cauchy stress

$\sigma =\mu (\lambda -{\lambda }^{-2})$

 

where $\mu$ is the shear modulus.

 

 In a nonlinear structure, it is usually helpful to calculate the joint displacements as the load is applied in a series of increments.   So, suppose that the joint positions $y_t^a$ at a time instant t are known, and we wish to determine the new displacements $y_{t+\Delta t}^a$ at time $t+\Delta t$.   The virtual work equation at $t+\Delta t$ is

${\displaystyle \sum _{k=1}^m\left[\begin{array}{cc}{v}^a& v^b\end{array}\right]\left[\begin{array}{c}{y}_{t+\Delta t}^a-y_{t+\Delta t}^b\\ y_{t+\Delta t}^b-y_{t+\Delta t}^a\end{array}\right]}{\displaystyle \frac{A_0^k}{L_0^k{\lambda }^2}\sigma (\lambda )={\displaystyle \sum _{j=1}^n{P}^j\cdot v^j}}$

where we have noted that $L=\lambda L_0$ ( $\lambda$ is the stretch of the kth member), and that $AL=A_0{L}_0$ because the material is incompressible.  The virtual equation must be satisfied for all admissible virtual velocities, which yields a set of nonlinear equations

$R={\displaystyle \sum _{j=1}^n{P}^j}-{\displaystyle \sum _{k=1}^m\left[\begin{array}{c}{y}_{t+\Delta t}^a-y_{t+\Delta t}^b\\ y_{t+\Delta t}^b-y_{t+\Delta t}^a\end{array}\right]}{\displaystyle \frac{A_0^k}{L_0^k}{\displaystyle \frac{\sigma (\lambda )}{{\lambda }^2}=0}}$

that must be solved numerically.

 

Nonlinear equilibrium equations in finite element analysis are usually solved using Newton-Raphson iteration.   For a truss, this means that we repeatedly correct a current approximation for the joint displacements $w$ by an amount $\Delta w$ until the equilibrium equation is satisfied.   Ideally the correction should satisfy $R^a(w+\Delta w)=0$, but of course it is not possible to solve this equation, so instead we take a Taylor expansion and rearrange

$K\Delta w=R(w)K=-{\displaystyle \frac{dR}{dy}}$

This is now a set of linear equations which, together with equations for joints with prescribed displacements, can be solved for $\Delta w$.   The stiffness matrix (prior to applying constraints) is

$K={\displaystyle \sum _{k=1}^m{\displaystyle \frac{A_0^k}{L_0^k}\left({\displaystyle \frac{d\sigma }{d\lambda }-{\displaystyle \frac{2\sigma }{\lambda }}}\right){\displaystyle \frac{1}{L_0^{k2}{\lambda }^3}\left[\begin{array}{c}{w}^a-w^b\\ w^b-w^a\end{array}\right]}}}\otimes \left[\begin{array}{cc}{w}^a-w^b& w^b-w^a\end{array}\right]+{\displaystyle \frac{A_0^k}{L_0^k}{\displaystyle \frac{\sigma (\lambda )}{{\lambda }^2}\left[\begin{array}{cc}I& -I\\ -I& I\end{array}\right]}}$

where I represents a (3x3) identity matrix (for a 3D structure), and the summation represents the assembly of the stiffness matrix for the k th element (which has nodes (a,b) at its two ends) into a global stiffness.

 

The Newton-Raphson equation solver usually works well, but if the joint deflections are so large that the structure ‘snaps through’ an unstable equilibrium configuration, it may fail to converge.   In this case a simple iterative relaxation solver may work better.  This approach involves the following steps.  Start with an initial guess for the joint coordinates w (which is usually the positions at the end of the preceding load-step).  Then:

 

1. Update the coordinates of any joints with prescribed displacements to their deformed positions;

 

2. Compute the force vector R;

 

3. Set rows of R corresponding to degrees of freedom with known displacements to zero

 

4. Check for convergence $\left|R\right|<TOL$ (where TOL is an appropriate small tolerance).   If the solution has converged update the joint coordinates to y=w and proceed to the next load-step.

 

5. If the solution has not yet converged, updated the approximate joint positions to $w=w+\gamma R$, where $\gamma$ is a numerical relaxation factor (for the neo-hookean material $\gamma \approx 1/\max (A\mu /L)$, where the max is taken over all members) and proceed to step 2.

 

As an example, the figure above illustrates the predicted deformation of a 3 noded truss structure (with members made from a hyperelastic material with modulus $\mu$ ) subjected to horizontal and vertical forces with ratio $P_2/P_1=6$ .   Fig. (a) shows the shape of the deformed structure for various magnitudes of the applied force; while the graph (b) shows the displacement of the loaded joint as a function of the magnitude of the force.

 

 

Trusses made from elastic-plastic materials that experience large displacements: Finally, we show how to analyze a truss that has its members made from a metallic material, which displays elastic behavior at low stresses, and experiences permanent plastic deformation if the stresses exceed yield. 

 

Specifically, consider a truss with members made from an elastic $–$ rate independent plastic material that displays power law hardening.  For this type of material, the stretch $\lambda$ in the members can be separated into a reversible elastic part ${\lambda }^e$, and a permanent plastic part ${\lambda }^p$ , so that $\lambda ={\lambda }^e{\lambda }^p$.   The yield stress of the members is related to the accumulated plastic strain ${\epsilon }^p$ by

$Y=Y_0{(1+{\epsilon }^p)}^m$

where $Y_0$ is the initial yield stress, m<1 is the strain hardening exponent, and the accumulated plastic strain is related to the plastic stretch by

${\displaystyle \frac{d{\epsilon }^p}{dt}={\displaystyle \frac{1}{{\lambda }^p}\left|{\displaystyle \frac{d{\lambda }^p}{dt}}\right|}}$

The uniaxial Cauchy (true) stress in the members is related to the elastic part of the stretch by

$\sigma =E\log ({\lambda }^e)$

where E is the Young’s modulus.  The plastic stretch rate is zero if the stress is below the yield stress, or if the stress magnitude is decreasing.   Otherwise, the plastic stretch rate is given by

${\displaystyle \frac{d{\lambda }^p}{dt}={\displaystyle \frac{(1+{\epsilon }^p){\lambda }^p}{mY}{\displaystyle \frac{d\sigma }{dt}}}}$

 

Changes in volume of the members can be neglected. 

 

The virtual work equation and the Newton-Raphson iterations for an elastic-plastic truss have the same form as those listed for a hyperelastic truss in the preceding section.   Only the expressions for the stress and its derivative with respect to stretch need to be changed.   

 

The stress-strain relation for a plastically deforming material depends on the history of loading, so the finite element solution must calculate the gradual change in shape of the structure as the load is applied in a series of increments.  Suppose that at the start of the nth time increment the plastic stretch and the total accumulated plastic strain have values ${\lambda }_n^p,{\epsilon }_n^p$.   Given the total stretch $\lambda$ at the end of the time step we must now calculate $\sigma ,{\lambda }_{n+1}^p,{\epsilon }_{n+1}^p$ at the end of the time-step.   This can be accomplished using the following (implicit) stress update:

 

1.       Calculate the stress that would occur if the bar remains elastic, which is given by

$\sigma =E\left(\log \lambda -\log {\lambda }_n^p\right)$.

 

2.       Check whether the stress is below the yield stress: $\left|\sigma \right|\le Y_0{(1+{\epsilon }_n^p)}^m$.  If this condition is satisfied, the bar remains elastic and the stress is given by the formula in step 1.   The tangent is $d\sigma /d\lambda =E/\lambda$.

 

3.       If the estimate for stress in step 1 exceeds the yield stress, the plastic stretch and the accumulated plastic strain follow from the condition that the stress must equal the yield stress.  This requires

$\begin{array}{c}E\log (\lambda /{\lambda }_{n+1}^p)=\text{sign}(\sigma )Y_0{\left(1+{\epsilon }_{n+1}^p\right)}^m\\ {\epsilon }_{n+1}^p={\epsilon }_n^p+\text{sign}(\sigma )\log ({\lambda }_{n+1}^p/{\lambda }_n^p)\end{array}$

These two equations can be combined to eliminate ${\epsilon }_{n+1}^p$ and then solved for ${\lambda }_{n+1}^p$ using Newton-Raphson iteration (or any convenient numerical method).  

 

4.       The stress follows as $\sigma =E\left(\log \lambda -\log {\lambda }_{n+1}^p\right)$

 

5.       Finally, a straightforward but tedious calculation shows that the derivative of Cauchy stress with respect to stretch is

${\displaystyle \frac{d\sigma }{d\lambda }={\displaystyle \frac{E}{\lambda }{\displaystyle \frac{m{Y}_0{(1+{\epsilon }_{n+1}^p)}^{m-1}}{E+m{Y}_0{(1+{\epsilon }_{n+1}^p)}^{m-1}}}}}$

 

A truss with elastic-plastic members may experience a sudden loss of stability (particularly if the truss is statically determinate).   Because of the instability, the Newton-Raphson solution to the virtual work equation may fail to converge.    The relaxation method described in the preceding section can sometimes produce a better result.

 

 

As an example, the figure illustrates the predicted deformation of a 3 noded truss structure (with members made from an elastic-plastic material with yield stress $Y_0$, Young’s modulus $E=200Y_0$ and hardening exponent $m=0.25$ ) subjected to horizontal and vertical forces with ratio $P_2/P_1=10$.    Fig. (a) shows the shape of the deformed structure for various magnitudes of the applied force; while graph (b) shows the displacement of the loaded joint as a function of the magnitude of the force.

 

 

 

8.7.2 Beam Elements.

 

A beam is a slender member, with uniform (or slowly varying) cross-section, and length much greater than any cross-sectional dimension.   It may be subjected to forces and moments at its ends, or transverse forces along its length.   These cause the beam to bend, twist and stretch.  But because it has a slender cross-section, the stress and strain fields in a beam have a simple form.   There are two approaches to analyzing deformation in an assembly of beams using the finite element method.  The first is to set up and solve a finite element approximation to the classical differential equations governing deformation of a beam.  These equations were derived by approximating the strain field inside the beam in some appropriate way, so these approximations are built into the finite element solution from the beginning.  This idea has two limitations: (i) it works only if the beam remains elastic; and (ii) the most common version of classical beam theory (Euler-Bernoulli theory) assumes that planes transverse to the axis of the undeformed beam remain transverse to the deformed beam axis.   This is accurate for beams with length exceeding about 10 times the beam thickness or width, but not for shorter beams $–$ for short beams, ‘shear flexible’ (Timoshenko) beam theory gives better results.   It is also possible to get around these limitations by meshing a beam with standard 3D (or 2D) isoparametric elements, but the mesh would need to include at least 5-10 elements through the thickness of the beam, and since the elements would need to remain roughly cubic or square to avoid locking, a huge number of elements would needed even to analyze a single beam.  To avoid this, a better approach is to modify isoparametric elements so that they produce the correct strain field in a slender beam.   These are called ‘continuum beam’ elements.   Continuum beam theory in 3D is too complicated to discuss here, but we will show how to implement finite element versions of both Euler-Bernoulli and Timoshenko beam theory.

 

 

Summary of beam theory: Beam theory is discussed in (tedious) detail in Chapter 10, but to understand how to solve the equations using the finite element method, only a short summary is needed. 

 

The figure shows a beam.   We will assume that the beam is straight, with a uniform cross-section, and is made from an isotropic elastic material with Young’s modulus E and shear modulus $\mu$.  It experiences only small deformations.  It is often assumed that the beam does not stretch or twist, but in the version described here we will account (approximately) for both, but we will assume that sections of the beam that are transverse to its axis remain plane.  In practice, twisting the rod will cause its cross-section to warp.   The procedure necessary to correct for this effect in an elastic beam is discussed in more detail in Chapter 10.

 

The geometry of the beam is usually described using a basis $\{e_1,e_2,e_3\}$ shown in the figure, with $e_3$ is parallel to the beam’s axis, and $\{e_1,e_2\}$ orientated in some convenient direction with respect to its cross section.   We then define the following geometrical quantities:

 

1. The cross-sectional area $A={\displaystyle \underset{A}{\int }dA}$ 

 

2. The position of the neutral line in the cross section (this is a fiber in the beam that does not stretch as the beam bends)

$\overline{r}={\overline{x}}_1{e}_1+{\overline{x}}_2{e}_2={\displaystyle \frac{1}{A}{\displaystyle \underset{A}{\int }\left(x_1{e}_1+x_2{e}_2\right)dA}}$

 

3. The components of the area moment of inertia tensor for the cross-section

$\begin{array}{l}I=\left[\begin{array}{ccc}{I}_{11}& -I_{12}& 0\\ -I_{12}& I_{22}& 0\\ 0& 0& I_{33}\end{array}\right]\\ I_{11}={\displaystyle \underset{A}{\int }{\left(x_2-{\overline{x}}_2^{}\right)}^{2}dA}{I}_{22}={\displaystyle \underset{A}{\int }{\left(x_1-{\overline{x}}_1^{}\right)}^{2}dA}{I}_{33}=I_{11}+I_{22}\\ I_{12}={\displaystyle \underset{A}{\int }\left(x_1-{\overline{x}}_1^{}\right)\left(x_2-{\overline{x}}_2^{}\right)dA}\end{array}$

 

The effects of warping can be included by replacing $I_{33}$ by a modified polar moment of area.   The procedure is discussed in detail in Chapter 10.    

 

 

The deformed shape of a straight beam is described by specifying the displacement $u(x_3)$ of each point on the neutral line of the cross-section, as a function of position along the length of the beam. The axial stretch of the neutral section is related to the displacement by

${\epsilon }_{33}=e_3\cdot {\displaystyle \frac{du}{d{x}_3}={\displaystyle \frac{d{u}_3}{d{x}_3}}}$

It is also necessary to calculate the (small) angle of rotation of the beam’s cross-section.   The rotation caused by twisting the beam is quantified by a small rotation ${\theta }_3$ about the $e_3$ axis (with right hand screw convention).  Two different approaches are used to calculate the rotation of the cross-section about the transverse axes:

 

1. Euler-Bernoulli theory (which is valid for long, slender beams with length greater than about 10 times its cross-sectional dimensions) assumes that planes transverse to the beam’s neutral line before deformation remain transverse as the beam bends, so they rotate through small angles

${\theta }_1=-{\displaystyle \frac{d{u}_2}{d{x}_3}{\theta }_2={\displaystyle \frac{d{u}_1}{d{x}_3}}}$

The sign conventions used here are confusing $–$ the angles represent small rotations, in radians, about the $e_1,e_2$ axes.  The small rotation is a vector, and is related to the displacement by $\theta =e_3\times du/d{x}_3+{\theta }_3{e}_3$.    Taking cross products of both sides of this expression with $e_3$ shows also that $du/d{x}_3=-e_3\times \theta +\left(e_3\cdot du/d{x}_3\right)e_3$

 

2. Timoshenko beam theory (which is valid for shorter beams, but will also work for long ones) allows the cross-section to rotate relative to the neutral line. 

 

In both theories, the curvature and twist per unit length of the beam are quantified by a vector $\kappa =d\theta /d{x}_3$ 

 

A beam may be loaded by subjecting its ends to forces $P^{(0)},P^{(L)}$ or moments $Q^{(0)},Q^{(L)}$.  Alternatively, one or more components of displacement or rotation may be prescribed at the ends of the beam.  In addition, a distributed force per unit length $p(x_3)$ may act along the beam’s length.

 

The internal forces in a beam are quantified by internal force and moment vectors acting on each cross-section.  To make this precise, introduce an imaginary cut perpendicular to the axis of the beam, as shown in the figure. The stresses acting on an interior face with normal parallel to the $e_3$ direction exert resultant forces T, and bending moments M. In an elastic beam with Youngs modulus E and shear modulus $\mu$, the bending moment vector is related to the curvature vector by

$\left[\begin{array}{c}{M}_1\\ M_2\\ M_3\end{array}\right]=\left[\begin{array}{ccc}E{I}_{11}& -E{I}_{12}& 0\\ -E{I}_{12}& E{I}_{22}& 0\\ 0& 0& \mu I_{33}\end{array}\right]\left[\begin{array}{c}{\kappa }_1\\ {\kappa }_2\\ {\kappa }_2\end{array}\right]$

Here, $I_{33}$ may be replaced by a modified polar moment of area to correct for warping of the cross-section.

 

For an Euler-Bernoulli beam the force vector T is a constraint (reaction) force and must be calculated using the equilibrium equation.  For a stretchable Timoshenko beam the forces are related to the displacement of the neutral line and the rotation vector by

$T=\beta A\mu \left({\displaystyle \frac{du}{d{x}_3}+e_3\times \theta }\right)+\left((E-\beta \mu )A{e}_3\cdot {\displaystyle \frac{du}{d{x}_3}}\right)e_3$

where $\beta$ is a constant known as the ‘Timoshenko shear coefficient’ for the beam. Since Timoshenko beam theory is (like all beam theories) an approximation, the value of $\beta$ depends on how it is calculated.  It depends on the geometry of the cross-section and Poisson’s ratio, but for most practical applications one can assume $\beta \approx 5/6$.

 

The static equilibrium equations for the internal forces are

${\displaystyle \frac{dT}{d{x}_3}+p=0{\displaystyle \frac{dM}{d{x}_3}+e_3\times T=0}}$

As usual the finite element method replaces these by the equivalent principle of virtual work.  For a Timoshenko beam these are

$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}T\cdot {\displaystyle \frac{d\delta u}{d{x}_3}}}d{x}_3={\displaystyle \underset{0}{\overset{L}{\int }}p}\cdot \delta ud{x}_3+P^{(L)}\cdot \delta u(L)+P^{(0)}\cdot \delta u(0)\\ {\displaystyle \underset{0}{\overset{L}{\int }}M\cdot {\displaystyle \frac{d\delta \theta }{d{x}_3}}}d{x}_3={\displaystyle \underset{0}{\overset{L}{\int }}{e}_3\times T\cdot \delta \theta d{x}_3}+Q^{(L)}\cdot \delta \theta (L)+Q^{(0)}\cdot \delta \theta (0)\end{array}$

For the Euler-Bernoulli beam the components of T acting transverse to the beam must be eliminated from the virtual work equation, which gives

$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}\left(M\cdot {\displaystyle \frac{d\delta \theta }{d{x}_3}+T_3\delta {\epsilon }_{33}}\right)}d{x}_3={\displaystyle \underset{0}{\overset{L}{\int }}p\cdot \delta u}\\ +Q^{(0)}\cdot \delta \theta (0)+Q^{(L)}\cdot \delta \theta (L)+P^{(0)}\cdot \delta u(0)+P^{(L)}\cdot \delta u(L)\end{array}$

Here, $\delta u$ and $\delta \theta$ are kinematically admissible variations in displacement and rotation, which must be related by the constraint equation $\delta \theta =-e_3\times d\delta u/d{x}_3+\delta {\theta }_3{e}_3$

 

Finite element method for elastic Euler-Bernoulli beams: There are three new features in the virtual work equation for an Euler-Bernoulli beam:

 

1. The equation includes moments.  These are related to the second derivative of the displacements, or equivalently, the first derivative of the components of the rotation vector.   The appearance of this second derivative means that the displacement and virtual displacement field must be continuous for the finite element method to converge as the mesh size is reduced.

 

2. The virtual work equation includes the rotation vector;

 

3. The virtual work equation requires that the rotation is correctly related to the displacement.

 

These new features mean that we must use a special scheme to interpolate the displacement field.   In fact, we solve not only for the displacements, but also for the rotation, at discrete points along the beam.

 

To this end, the beam is divided into a series of elements bounded by nodes, as shown in the figure.   The unknowns in the finite element solution are the three components of displacement, and the three components of rotation, at each node.  We will set up a scheme that allows these components to be expressed in an arbitrary $\{i,j,k\}$ basis (which allows the beam to have any orientation in space), but for now, suppose that we wish to solve for their components in the $\{e_1,e_2,e_3\}$ basis aligned with the beam.

 

The axial displacement, and the axial component of rotation (i.e. the twist) can then be interpolated linearly, since these are differentiated only once in the virtual work equation.   The transverse components of displacement and rotation must be interpolated using cubic Hermitian polynomials, which ensure that the rotation and displacement are continuous across neighboring elements, and also that the rotation is correctly related to the displacement.    The interpolations can be conveniently written as a matrix expression

 

where $u_i^1,u_i^2$ and ${\theta }_i^1,{\theta }_i^2$ represent the displacements and rotations at the two nodes of the element, and the interpolation functions are

$\begin{array}{l}{N}^1=(1-\xi )/2N^2=(1+\xi )/2\\ M_u^1={\left(\xi -1\right)}^2\left(\xi +2\right)/4M_u^2={\left(\xi +1\right)}^2\left(2-\xi \right)/4\\ M_{\theta }^1={\left(\xi -1\right)}^2\left(\xi +1\right)/8M_{\theta }^2={\left(\xi +1\right)}^2\left(\xi -1\right)/8\\ {M^{\prime }}_u^1=3\left(\xi -1\right)\left(\xi +1\right)/(2L_e){M^{\prime }}_u^2=3\left(\xi +1\right)\left(1-\xi \right)/(2L_e)\\ {M^{\prime }}_{\theta }^1=\left(3{\xi }^2-2\xi -1\right)/(4L_e){M^{\prime }}_{\theta }^2=\left(3{\xi }^2+2\xi -1\right)/(4L_e)\\ L_e=\left|x^2-x^1\right|\end{array}$

with $-1<\xi <1$ representing the normalized distance along the element. The curvature and axial strain can then be computed as follows

 

where

$\begin{array}{l}{N^{\prime }}^1=-1/L_e{N^{\prime }}^2=1/L_e\\ {M^{″}}_u^1=6\xi /L_e^2{M^{″}}_u^2=-6\xi /L_e^2\\ {M^{″}}_{\theta }^1=\left(3\xi -1\right)/L_e^2{M^{″}}_{\theta }^2=\left(3\xi +1\right)/L_e^2\end{array}$

The bending moment, twisting moment and axial force can then be calculated as

$\begin{array}{l}{\left[\begin{array}{cccc}{M}_1& M_2& M_3& T_3\end{array}\right]}^T=D{\left[{\kappa }_1,{\kappa }_2,{\kappa }_3,{\epsilon }_{33}\right]}^T\\ D=\left[\begin{array}{cccc}E{I}_{11}& -E{I}_{12}& 0& 0\\ -E{I}_{12}& E{I}_{22}& 0& 0\\ 0& 0& \mu I_{33}& 0\\ 0& 0& 0& EA\end{array}\right]\end{array}$

Finally, in many applications (particularly those involving several beams welded together to create a structure) it is convenient to solve for displacement and rotation components in a global $\{i,j,k\}$ basis, instead of using the components in the $\{e_1,e_2,e_3\}$ basis, which would differ for each beam.   This can be accomplished by defining a basis change matrix

$\Omega =\left[\begin{array}{cccc}R& 0& 0& 0\\ 0& R& 0& 0\\ 0& 0& R& 0\\ 0& 0& 0& R\end{array}\right]R=\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\end{array}\right]$

where 0 represents a 3x3 null matrix, and $e_1,e_2,e_3$ represents the components of these vectors in the global $\{i,j,k\}$ basis.

 

With these definitions, the principle of virtual work reduces to a system of linear equations

$KU=FK={\displaystyle \sum _{elements}{k}_{el}}F={\displaystyle \sum _{nodes}{f}_n}+{\displaystyle \sum _{elements}{r}_{el}}$

where U is a vector of unknown nodal displacements and rotations (as components in the  $\{i,j,k\}$ basis), F is a vector of nodal forces and moments, and

$\begin{array}{l}{k}_{el}={\displaystyle \underset{0}{\overset{L_e}{\int }}{\Omega }^T{B}^{T}DB\Omega d{x}_3}\approx {\displaystyle \frac{L_e}{2}{\displaystyle \sum _i{\Omega }^T{B}^T({\xi }_i)DB({\xi }_i)\Omega w_i}}\\ r_{el}={\displaystyle \underset{0}{\overset{L_e}{\int }}{\Omega }^T{N}^T\Omega pd{x}_3}\approx {\displaystyle \frac{L_e}{2}{\displaystyle \sum _i{\Omega }^T{N}^T({\xi }_i)\Omega p{w}_i}}\end{array}$

are the element stiffness matrix and force vector, with ${\xi }_i,w_i$ the coordinates and weights for a 2 point Gauss integration scheme ( ${\xi }_1=-1/\sqrt{3},{\xi }_2=1/\sqrt{3}$ and $w_1=w_2=1$ ).

 

 

 

Finite element method for elastic Timoshenko beams: The principle of virtual work for the Timoshenko beam consist of two coupled virtual work equations for the unknown displacement and rotation of transverse cross-sections of the beam.    The two are independent variables, so

 

1. There is no need to ensure that the rotation are related by a compatibility constraint;

 

2. The displacement and rotation must be continuous across element boundaries, but there is no need for continuity of the derivative of the displacement.

 

This means that standard linear or quadratic interpolation functions can be used to approximate their variations within an element.  Higher order interpolation schemes can also be found in the literature.

 

To proceed with the finite element formulation, it is helpful to combine the two virtual work equations for the Timoshenko beam into a single, symmetric form.   For this purpose we define a shear/axial strain vector

$\gamma ={\displaystyle \frac{du}{d{x}_3}+e_3\times \theta }$

Adding the two virtual work equations then yields

$\begin{array}{l}{\displaystyle \underset{0}{\overset{L}{\int }}M\cdot {\displaystyle \frac{d\delta \theta }{d{x}_3}}}d{x}_3+{\displaystyle \underset{0}{\overset{L}{\int }}T\cdot \delta \gamma d{x}_3}={\displaystyle \underset{0}{\overset{L}{\int }}p}\cdot \delta ud{x}_3+P^{(L)}\cdot \delta u(L)\\ +P^{(0)}\cdot \delta u(0)+Q^{(L)}\cdot \delta \theta (L)+Q^{(0)}\cdot \delta \theta (0)\end{array}$

This must hold for all admissible variations of $\delta u,\delta \theta$.

 

As before, the finite element version of the virtual work principle is derived by dividing the beam into a series of elements bounded by nodes, as shown in the figure.   The unknowns in the finite element solution are the three components of displacement, and the three components of rotation, at each node that is adjacent to two neighboring elements. We will again set up a scheme that allows the displacement and rotation components to be expressed in an arbitrary $\{i,j,k\}$ basis (which allows the beam to have any orientation in space), but begin by assuming that we wish to solve for their components in the $\{e_1,e_2,e_3\}$ basis aligned with the beam.

 

The displacement and rotation of the beam may be interpolated linearly between the nodes as follows

$\begin{array}{l}{\left[\begin{array}{cccccc}{u}_1& u_2& u_3& {\theta }_1& {\theta }_2& {\theta }_3\end{array}\right]}^T=\\ N{\left[\begin{array}{cccccccccccc}{u}_1^1& u_2^1& u_3^1& {\theta }_1^1& {\theta }_2^1& {\theta }_3^1& u_1^2& u_2^2& u_3^2& {\theta }_1^2& {\theta }_2^2& {\theta }_3^2\end{array}\right]}^T\\ N=\left[\begin{array}{cccccccccccc}{N}_1& 0& 0& 0& 0& 0& N_2& 0& 0& 0& 0& 0\\ 0& N_1& 0& 0& 0& 0& 0& N_2& 0& 0& 0& 0\\ 0& 0& N_1& 0& 0& 0& 0& 0& N_2& 0& 0& 0\\ 0& 0& 0& N_1& 0& 0& 0& 0& 0& N_2& 0& 0\\ 0& 0& 0& 0& N_1& 0& 0& 0& 0& 0& N_2& 0\\ 0& 0& 0& 0& 0& N_1& 0& 0& 0& 0& 0& N_2\end{array}\right]\end{array}$

Where

$N_1=(1-\xi )/2N_2=(1+\xi )/2$

The curvature and shear/axial strain vectors can then be calculated using

 

The internal forces and moments follow as

$\begin{array}{l}{\left[\begin{array}{cccccc}{T}_1& T_2& T_3& M_1& M_2& M_3\end{array}\right]}^T=D{\left[\begin{array}{cccccc}{\gamma }_1& {\gamma }_2& {\gamma }_3& {\kappa }_1& {\kappa }_2& {\kappa }_3\end{array}\right]}^T\\ D=\left[\begin{array}{cccccc}\beta A\mu & 0& 0& 0& 0& 0\\ 0& \beta A\mu & 0& 0& 0& 0\\ 0& 0& EA& 0& 0& 0\\ 0& 0& 0& E{I}_{11}& -E{I}_{12}& 0\\ 0& 0& 0& -E{I}_{12}& E{I}_{22}& 0\\ 0& 0& 0& 0& 0& \mu I_{33}\end{array}\right]\end{array}$

 

Finally, to re-write the finite element equations in terms of displacement and rotation components in a global $\{i,j,k\}$ basis we once again define thebasis change matrix

$\Omega =\left[\begin{array}{cccc}R& 0& 0& 0\\ 0& R& 0& 0\\ 0& 0& R& 0\\ 0& 0& 0& R\end{array}\right]R=\left[\begin{array}{c}{e}_1\\ e_2\\ e_3\end{array}\right]$

where 0 represents a 3x3 null matrix, and $e_1,e_2,e_3$ represents the components of these vectors in the global $\{i,j,k\}$ basis.

 

With these definitions, the principle of virtual work reduces to a system of linear equations

$KU=FK={\displaystyle \sum _{elements}{k}_{el}}F={\displaystyle \sum _{nodes}{f}_n}+{\displaystyle \sum _{elements}{r}_{el}}$

where U is a vector of unknown nodal displacements and rotations (as components in the  $\{i,j,k\}$ basis), F is a vector of nodal forces and moments, and

$\begin{array}{l}{k}_{el}={\displaystyle \underset{0}{\overset{L_e}{\int }}{\Omega }^T{B}^{T}DB\Omega d{x}_3}\approx L_e{\Omega }^T{B}^T(0)DB(0)\Omega \\ r_{el}={\displaystyle \underset{0}{\overset{L_e}{\int }}{\Omega }^T{N}^T\Omega pd{x}_3}\approx L_e{\Omega }^T{N}^T(0)\Omega p\end{array}$

are the element stiffness matrix and force vector.  The integrals have been evaluated with a 1-point Gaussian quadrature scheme, in which the integrand is calculated at $\xi =0$, and the integration weight is $w=1$.   A higher order integration scheme cannot be used, because it causes the element to ‘lock,’ causing a spuriously high element stiffness matrix.  This is particularly a problem in long, slender beams with little shear deformation.

 

Example: A simple demonstration of both the Euler-Bernoulli and Timoshenko beam elements are shown in the figure.  Fig. (a) compares the predictions of the Euler-Bernoulli beam element for the deflection and rotation of an end-loaded cantilever beam with the exact solution.  The FEA predictions are within 0.1% of the exact result with only 10 elements along the length of the beam.   Fig. (b) shows the predicted deflection of an end loaded Timoshenko beam, for several values of the normalized ratio of bending to shear stiffness $EI/(\beta \mu A{L}^2)$.   For $EI/(\beta \mu A{L}^2)\to 0$ the Timoshenko elements converge to the limit of the Euler-Bernoulli beam.   As $EI/(\beta \mu A{L}^2)>0.02$ shear deformation in the beam becomes significant, and the Timoshenko elements predict a significantly larger deflection than Euler-Bernoulli theory.

 

The codes that plot these graphs are provided in the files named FEM_beam_Euler.m and FEM_beam_Timoshenko.m at

https://github.com/albower/Applied_Mechanics_of_Solids/

 

 

 

8.7.3 Plate Elements.

 

We next consider simplified solid mechanics theories that describe motion of flat plates.   The figure shows the problem to be solved.  Our goal is to calculate the deflections and forces in a plate that is subjected to transverse forces on its surface, and either supported or clamped on its edge, or subjected to in-plane or transverse forces around its circumference.   To keep the discussion (reasonably) simple here, we will assume

 

· The solid is a flat plate, with a uniform thickness h (in more sophistication FEA analysis, plates can be curved - then they are called shells -  and the thickness can vary)

 

· Deflections are small, and the plate remains elastic.

 

· The plate is oriented so that the $e_3$ direction is normal to the undeformed mid-section of the plate.   In more complex FEA simulations the plate may have an arbitrary orientation: in this case the element stiffness matrices for the plate are first formulated using a local coordinate system with $e_3$ normal to the plate, and then transformed to a global coordinate frame.   The procedure was discussed in Section 8.7.2 for beams, but to keep things brief will not be repeated in this section.

 

 

Summary of plate theory: Plate theory is discussed in (tedious) detail in Chapter 10, but to understand how to solve the equations using the finite element method, only a short summary is needed. 

 

The deformed shape of a flat plate is described by specifying the displacement $u(x_1,x_2)$ of each point on the mid-section of the plate as a function of position. 

 

It is also necessary to calculate the (small) rotation of the planes normal to the $e_1$ and $e_2$ directions in the undeformed plate.   The rotation is quantified by a small rotation vector $\theta$, using one of two different approaches:

 

1. Kirchhoff plate theory is an approximate theory (analogous to the Euler-Bernoulli beam), which is intended to model plates that are thinner than 1/10^th of their in-plane dimensions.  It assumes that planes transverse to the mid-plane of the plate remain transverse after deformation.   In this case, the rotation is related to the out-of-plane deflection of the plate by

${\theta }_1={\displaystyle \frac{d{u}_3}{d{x}_2}{\theta }_2=-{\displaystyle \frac{d{u}_3}{d{x}_1}}}$

These angles represent the small rotation of a line that is transverse to the mid-plane of the undeformed plate.  They are components of a vector, and are related to the transverse displacement by $\theta =-e_3\times \nabla u_3$.

 

2. Reissner-Mindlin plate theory is a more complicated version, which allows the transverse planes to rotate relative to the mid-plane of the plate.   The rotation is then an independent variable, which must be calculated as part of the solution.

 

In both theories, the curvature of the plate is defined as

${\kappa }_{11}={\displaystyle \frac{d{\theta }_2}{d{x}_1}{\kappa }_{12}=-{\displaystyle \frac{d{\theta }_1}{d{x}_1}{\kappa }_{21}={\displaystyle \frac{d{\theta }_2}{d{x}_2}{\kappa }_{22}=-{\displaystyle \frac{d{\theta }_1}{d{x}_2}}}}}$

The curvature is a (2D) tensor, and is related to the rotation vector by $\kappa =-\nabla \left(e_3\times \theta \right)$ 

 

The strain distribution in the plate is then

${\epsilon }_{\alpha \beta }={\displaystyle \frac{1}{2}\left({\displaystyle \frac{\partial u_{\alpha }}{\partial x_{\beta }}+{\displaystyle \frac{\partial u_{\beta }}{\partial x_{\alpha }}+x_3\left({\kappa }_{\alpha \beta }+{\kappa }_{\beta \alpha }\right)}}\right)}$

where the Greek indices have values 1 or 2.

 

The plate is in a state of plane stress, with internal stresses

$$ $\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{12}\end{array}\right]={\displaystyle \frac{E}{(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ 2{\epsilon }_{12}\end{array}\right]}$

Reissner-Mindlin plate theory also accounts approximately for the effects of the transverse shear stress components ${\sigma }_{31},{\sigma }_{32}$, but these cannot be calculated directly from the strain field

 

A plate may be loaded by subjecting its edge to a force $P$ per unit length or moment $Q$ per unit length.  Alternatively, one or more components of displacement or rotation may be prescribed.  In addition, a distributed force per unit area $p(x)e_3$ may act perpendicular to surface of the plate.

 

The forces in a plate are quantified by internal force and moment tensors acting on each cross-section.  To make this precise, cut a square element from the plate with planes normal to the coordinate axes as shown in the figure.  A set of resultant forces and moments act on the exposed faces:

 

· The moments $M_{ij}$ are defined by

$M_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{11}{x}_{3}d{x}_3}{M}_{12}=M_{21}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{12}{x}_{3}d{x}_3}{M}_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{22}{x}_{3}d{x}_3}$ 

The physical significance of the components $M_{ij}$ is illustrated in the figure: $M_{1j}$ characterizes the moment per unit length acting on planes inside the shell that are normal to the $e_1$ direction, while $M_{2j}$ characterizes the moment per unit length acting on planes that are normal to $e_2$.  Note that $M_{i1}$ represents a moment about the $e_2$ axis, while $M_{i2}$ is a moment acting about the $-e_1$ axis. This can be expressed mathematically as $M=M_{\alpha \beta }{e}_{\alpha }\otimes (e_3\times e_{\beta })$ where the repeated indices $\alpha ,\beta$ are summed over 1 and 2, but this expression is only helpful if you are really good at visualizing dyadic and cross products.

 

· The in-plane forces $T_{ij}$ are defined by

$T_{11}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{11}d{x}_3}{T}_{12}=T_{21}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{12}d{x}_3}{T}_{22}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{22}d{x}_3}$

They represent resultant forces exerted by stresses on the internal planes: $T_{1j}$ are the components of force acting on the plane perpendicular to the $e_1$ direction, while $T_{2j}$ are those acting on the plane perpendicular to the $e_2$ direction. 

  

· The transverse forces $V_1,V_2$ represent forces acting in the $e_3$ direction on planes perpendicular to the $e_1$ and $e_2$ direction, respectively.   In Kirchhoff theory they are constraint forces and must be calculated from the equilibrium equation.   In an exact version of Reissner-Mindlin theory they would be related to the shear stresses acting on transverse planes by

 

$V_1={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{31}d{x}_3}{V}_2={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{32}d{x}_3}$

but the theory is approximate, so the transverse forces are calculated instead through the rotation of the cross-section, using the expressions given below.

 

 

The internal forces are related to the displacement and curvature of the mid-plane of the plate by

$\begin{array}{l}\left[\begin{array}{c}{T}_{11}\\ T_{22}\\ T_{12}\end{array}\right]={\displaystyle \frac{Eh}{(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]\left[\begin{array}{c}\partial u_1/\partial x_1\\ \partial u_2/\partial x_2\\ (\partial u_1/\partial x_2+\partial u_2/\partial x_1)\end{array}\right]}\\ \left[\begin{array}{c}{M}_{11}\\ M_{22}\\ M_{12}\end{array}\right]={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]\left[\begin{array}{c}{\kappa }_{11}\\ {\kappa }_{22}\\ {\kappa }_{12}+{\kappa }_{21}\end{array}\right]}\end{array}$

 

In Reissner-Mindlin theory, the transverse shear forces are related to the displacement and rotation by

$\left[\begin{array}{c}{V}_1\\ V_2\end{array}\right]=2\beta h\mu \left[\begin{array}{c}\partial u_3/\partial x_1+{\theta }_2\\ \partial u_3/\partial x_2-{\theta }_1\end{array}\right]$

 

where $\beta$ is a constant known as the ‘shear factor, ’ analogous to the shear coefficient for a beam.  Since the theory is approximate, the value of $\beta$ depends on how it is calculated.  It depends on Poisson’s ratio, but for most practical applications one can assume $\beta \approx 5/6$.

 

The static equilibrium equations for a plate are

${\displaystyle \frac{\partial V_1}{\partial x_1}+{\displaystyle \frac{\partial V_2}{\partial x_2}+p_3=0}}$

    $\begin{array}{l}{\displaystyle \frac{\partial M_{11}}{\partial x_1}+{\displaystyle \frac{\partial M_{12}}{\partial x_2}-T_{11}{\displaystyle \frac{\partial u_3}{\partial x_1}-T_{12}{\displaystyle \frac{\partial u_3}{\partial x_2}-V_1=0}}}}\\ {\displaystyle \frac{\partial M_{21}}{\partial x_1}+{\displaystyle \frac{\partial M_{22}}{\partial x_2}-T_{22}{\displaystyle \frac{\partial u_3}{\partial x_2}-T_{21}{\displaystyle \frac{\partial u_3}{\partial x_1}-V_2=0}}}}\end{array}$

${\displaystyle \frac{\partial T_{11}}{\partial x_1}+{\displaystyle \frac{\partial T_{21}}{\partial x_2}=0{\displaystyle \frac{\partial T_{12}}{\partial x_1}+{\displaystyle \frac{\partial T_{22}}{\partial x_2}=0}}}}$

This version of plate theory accounts approximately for the effects of large in-plane loads on the out-of-plane motion of the plate (so the theory can model a thin stretched membrane as well as a plate).   In simpler versions of plate theory (which are more commonly used in practice) the terms involving products of $T_{\alpha \beta }$ and $\partial u_3/\partial x_{\beta }$ in the moment equilibrium equation are neglected.

 

For a finite element solution the equilibrium equations must be re-cast as the equivalent principle of virtual work:

${\displaystyle \underset{A}{\int }V\cdot \nabla \delta u_{3}dA}={\displaystyle \underset{C}{\int }{P}_3}\delta u_{3}ds+{\displaystyle \underset{A}{\int }{p}_3\delta u_3}dA$

${\displaystyle \underset{A}{\int }M\cdot \delta \kappa dA}-{\displaystyle \underset{A}{\int }V\cdot }\left(e_3\times \delta \theta \right)dA+{\displaystyle \underset{A}{\int }\left(\nabla \delta u_3\right)\cdot T\cdot (\nabla u_3)dA}={\displaystyle \underset{C}{\int }Q\cdot }\delta \theta ds$

${\displaystyle \underset{A}{\int }T\cdot \nabla \delta wdA}={\displaystyle \underset{A}{\int }p\cdot \delta wdA}+{\displaystyle \underset{C}{\int }P}\cdot \delta wds$

For the Kirchhoff plate, V must be eliminated from the first two virtual work equations with the result

${\displaystyle \underset{A}{\int }M\cdot \delta \kappa dA}+{\displaystyle \underset{A}{\int }\left(\nabla \delta u_3\right)\cdot T\cdot (\nabla u_3)dA}={\displaystyle \underset{C}{\int }Q\cdot }\delta \theta ds+{\displaystyle \underset{C}{\int }{P}_3\delta u_{3}ds}+{\displaystyle \underset{A}{\int }{p}_3\delta u_{3}dA}$

 

 

Finite element method for in-plane displacements and forces: In many practical applications the in-plane displacements and forces in the plate may be neglected, but if needed, they can be calculated before the out-of-plane displacements.  Notice that the virtual work equation for the in-plane displacements w and forces T decouple from those for the out-of-plane displacements and rotations of the plate.   Moreover, the equation is identical (aside from the distinction between forces and stress) to that for plane stress deformation of a plate.   Consequently, the equation can be solved using the procedure described in Chapter 7, in which the in-plane displacement field is approximated by a linear variation between the three nodes on the element.   In this section we review the procedure briefly, for convenience.

 

A representative triangular plate element is illustrated in the figure.  It has three nodes with coordinates $x_k=(x_k{e}_1+y_k{e}_2)$, where k=1,2,3.  Each node has two in-plane displacement components $w^{(k)}=u_{\alpha }^{(k)}{e}_{\alpha }$, (the Greek indices have values 1 or 2), the out of plane displacement $u_3^{(k)}$ , and two rotation components ${\theta }^{(k)}={\theta }_{\alpha }^{(k)}{e}_{\alpha }$. 

 

The in-plane displacements and forces are then calculated as follows:

 

1.  Assemble the global stiffness matrix and force vector

$K={\displaystyle \sum _{elements}{k}_{el}}f={\displaystyle \sum _{elements}{f}_{el}}$

where the element stiffness and force vector are

$k_{el}=A_{el}{B}_0^T{D}_0{B}_0{f}_{el}={\displaystyle \frac{l_{el}}{2}\left[\begin{array}{c}{P}_1\\ P_2\end{array}\right]}$

Here,  $A_{el}$ is the area of the element, $l_{el}$ is the length of a loaded element side, and $P_{\alpha }$ the in-plane components of the (constant) force per unit length acting on the side of the element, while

 

where

$\begin{array}{c}{b}_1=y_2-y_3{b}_2=y_3-y_1{b}_3=y_1-y_2\\ c_1=x_3-x_2{c}_2=x_1-x_3{c}_3=x_2-x_1\end{array}$

and

$D_0={\displaystyle \frac{Eh}{(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]}$

 

2.  Modify the global stiffness to impose any constraints on the in-plane displacements at the boundary

 

3.  Solve the equation system $Kw=f$ for the in-plane displacements

 

4.  Calculate the in-plane forces (which are constant within each element) using $T=D_0{B}_0{w}_{el}$, where $w_{el}$ is the 6 dimensional vector of the in-plane displacements at the corners of the element.

 

 

Finite element method for out-of-plane deformation of elastic Kirchhoff plates: We next describe a method for calculating the out-of-plane deformation of a Kirchhoff plate.  The virtual work equation for the out-of-plane motion resembles that for a beam, in that the equation involves the curvature of the plate.   This means that ideally the out-of-plane displacements and its derivative must both be continuous across the boundaries of an element. For a beam, it was possible to find an interpolation scheme (using Hermitian polynomial shape functions) that satisfied this condition exactly.    There is unfortunately no equivalent interpolation scheme for a plate.   Consequently, a large number of different approaches to setting up interpolation schemes for plates have been developed over the years, and new ones continue to appear.  We will not attempt to review these here, and instead will summarize how to set up a simple triangular plate element called the ‘Specht triangle’ in honor of its author.   

 

 

The out-of-plane motion of the plate can be analyzed using the same mesh of triangular elements that is used to model its in-plane motion.   A representative triangular plate element is illustrated in Fig 8.54.   It has three nodes with coordinates $x_k=(x_k{e}_1+y_k{e}_2)$, where k=1,2,3.  Each node now has an unknown out of plane displacement $u_3^{(k)}$ , and two rotation components ${\theta }^{(k)}={\theta }_{\alpha }^{(k)}{e}_{\alpha }$. 

 

As always, the goal is to calculate the displacement $u_3$ and rotation ${\theta }_1,{\theta }_2$ at an arbitrary point $(x,y)$ within the element, given their values at the three corner nodes with coordinates $(x_i,y_i)$.  We define

$u_3=\left[\begin{array}{ccccccc}{N}_u^1& N_{\theta 1}^1& N_{\theta 2}^1& \cdots & N_u^3& N_{\theta 1}^3& N_{\theta 2}^3\end{array}\right]\left[\begin{array}{c}{u}^1\\ {\theta }_1^1\\ {\theta }_3^1\\ ⋮\\ u^3\\ {\theta }_1^3\\ {\theta }_3^3\end{array}\right]$

where the interpolation functions are calculated as follows:

 

1.  Calculate the ‘areal coordinates’ $L_i=A_i/A_{el}$, where $A_i$ is the area of the triangle formed by a point within the element and the two vertices opposite to node i, and $A_{el}$ is the area of the element:

$\begin{array}{l}{L}_1=\left[(x_2-x)(y_3-y)-(x_3-x)(y_2-y)\right]/(2A_{el})\\ L_2=\left[(x_3-x)(y_1-y)-(x_1-x)(y_3-y)\right]/(2A_{el})\\ L_3=\left[(x_1-x)(y_2-y)-(x_2-x)(y_1-y)\right]/(2A_{el})\\ 2A_{el}=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\end{array}$

 

2.  Assemble a vector of products of $L_i$

 

where

${\mu }_1={\displaystyle \frac{d_3-d_2}{d_1}{\mu }_2={\displaystyle \frac{d_1-d_3}{d_2}{\mu }_1={\displaystyle \frac{d_2-d_1}{d_3}}}}$

$d_1={\left(x_3-x_2\right)}^2+{\left(y_3-y_2\right)}^2{d}_2={\left(x_1-x_3\right)}^2+{\left(y_1-y_3\right)}^2{d}_3={\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2$

 

3.  The shape functions can then be calculated as weighted sums of the components of g:

$\left[\begin{array}{c}{N}_u^1\\ N_{\theta 1}^1\\ N_{\theta 2}^1\\ N_u^2\\ N_{\theta 1}^2\\ N_{\theta 1}^2\\ N_u^3\\ N_{\theta 1}^3\\ N_{\theta 2}^3\end{array}\right]=\left[\begin{array}{c}{g}_1-g_4+g_6+2(g_7-g_9)\\ -b_2(g_9-g_6)-b_3{g}_7\\ -c_2(g_9-g_6)-c_3{g}_7\\ g_2-g_5+g_4+2(g_8-g_7)\\ -b_3(g_7-g_4)-b_1{g}_8\\ -c_3(g_7-g_4)-c_1{g}_8\\ g_3-g_6+g_5+2(g_9-g_8)\\ -b_1(g_8-g_5)-b_1{g}_9\\ -c_1(g_8-g_5)-c_1{g}_9\end{array}\right]$

where (as before)

$\begin{array}{c}{b}_1=y_2-y_3{b}_2=y_3-y_1{b}_3=y_1-y_2\\ c_1=x_3-x_2{c}_2=x_1-x_3{c}_3=x_2-x_1\end{array}$

 

4.  The derivatives of the shape functions are also needed: the first derivative can be calculated as

${\displaystyle \frac{dg}{dx}=\left[\begin{array}{cc}d{g}_1/dx& d{g}_1/dy\\ ⋮& \\ d{g}_9/dx& d{g}_9/dy\end{array}\right]={\displaystyle \frac{dg}{dL}{\displaystyle \frac{dL}{dx}}}}$

where

$\begin{array}{l}{\displaystyle \frac{dg}{dL}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ L_2& L_1& 0\\ 0& L_3& L_2\\ L_3& 0& L_1\\ 2L_1{L}_2& L_1^2& 0\\ 0& 2L_2{L}_3& L_2^2\\ L_3^2& 0& 2L_3{L}_1\end{array}\right]+{\displaystyle \frac{3}{2}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 2(1-{\mu }_3)L_1{L}_2{L}_3& (1-{\mu }_3)L_1^2{L}_3& (1-{\mu }_3)L_1^2{L}_2\\ (1-{\mu }_1)L_2^2{L}_3& 2(1-{\mu }_1)L_1{L}_2{L}_3& (1-{\mu }_1)L_2^2{L}_1\\ (1-{\mu }_2)L_2{L}_3^2& (1-{\mu }_2)L_1{L}_3^2& 2(1-{\mu }_2)L_1{L}_2{L}_3\end{array}\right]}}\\ +{\displaystyle \frac{1}{2}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ (1+3{\mu }_3)(L_3-L_2)L_2{L}_3& (1+3{\mu }_3)(L_3-2L_2)L_1{L}_3& (1+3{\mu }_3)(2L_3-L_2)L_1{L}_2\\ (1+3{\mu }_1)(2L_1-L_3)L_2{L}_3& (1+3{\mu }_1)(L_1-L_3)L_1{L}_3& (1+3{\mu }_1)(L_1-2L_3)L_1{L}_2\\ (1+3{\mu }_2)(L_2-2L_1)L_2{L}_3& (1+3{\mu }_2)(2L_2-L_1)L_1{L}_3& (1+3{\mu }_2)(L_2-L_1)L_1{L}_2\end{array}\right]}\end{array}$

 

${\displaystyle \frac{dL}{dx}={\displaystyle \frac{1}{2A_{el}}\left[\begin{array}{cc}{b}_1& c_1\\ b_2& c_2\\ b_3& c_3\end{array}\right]2A_{el}=c_3{b}_2-b_3{c}_2}}$

 

5.  Similarly, the second derivatives are

${\displaystyle \frac{d^{2}g}{d{x}^2}=\left[\begin{array}{ccc}{d}^2{g}_1/d{x}^2& d^2{g}_1/d{y}^2& d^2{g}_1/dxdy\\ & ⋮& \\ d^2{g}_9/d{x}^2& d^2{g}_9/d{y}^2& d^2{g}_9/dxdy\end{array}\right]={\displaystyle \frac{d^{2}g}{d{L}^2}{\displaystyle \frac{d^{2}L}{d{x}^2}}}}$

Where

 

${\displaystyle \frac{d^{2}L}{d{x}^2}={\displaystyle \frac{1}{{\left(2A_{el}\right)}^2}\left[\begin{array}{ccc}{b}_1^2& c_1^2& b_1{c}_1\\ b_2^2& c_2^2& b_2{c}_2\\ b_3^2& c_3^2& b_3{c}_3\\ 2b_1{b}_2& 2c_1{c}_2& b_{1}c{}_2+b_2{c}_1\\ 2b_1{b}_3& 2c_1{c}_3& b_1{c}_3+b_3{c}_1\\ 2b_2{b}_3& 2c_2{c}_3& b_2{c}_3+b_3{c}_2\end{array}\right]2A_{el}=c_3{b}_2-b_3{c}_2}}$

The curvature and rotation of the plate can then be calculated as

$\begin{array}{l}{\left[-{\theta }_2,{\theta }_1\right]}^T=B_{\theta }[\begin{array}{ccccccc}{u}_3^1& {\theta }_1^1& {\theta }_2^1& \cdots & u_3^3& {\theta }_1^3& {\theta }_2^3{]}^T\end{array}\\ {\left[{\kappa }_{11},{\kappa }_{22},2{\kappa }_{12}\right]}^T=B_{\kappa }[\begin{array}{ccccccc}{u}_3^1& {\theta }_1^1& {\theta }_2^1& \cdots & u_3^3& {\theta }_1^3& {\theta }_2^3{]}^T\end{array}\end{array}$

where

$B_{\theta }=\left[\begin{array}{ccccccc}\partial N_u^1/\partial x_1& \partial N_{\theta 1}^1/\partial x_1& \partial N_{\theta 2}^1/\partial x_1& \cdots & N_u^3/\partial x_1& \partial N_{\theta 1}^3/\partial x_1& \partial N_{\theta 2}^3/\partial x_1\\ \partial N_u^1/\partial x_2& \partial N_{\theta 1}^1/\partial x_2& \partial N_{\theta 2}^1/\partial x_2& \cdots & \partial N_u^3/\partial x_2& \partial N_{\theta 1}^3/\partial x_2& \partial N_{\theta 2}^3/\partial x_2\end{array}\right]$

 

The bending moments follow as

${\left[M_{11},M_{22},M_{12}\right]}^T=D_{\kappa }{\left[{\kappa }_{11},{\kappa }_{22},2{\kappa }_{12}\right]}^T$

where

$D_{\kappa }={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]}$

 

 

With these definitions, the principle of virtual work reduces to a system of linear equations

$KU=FK={\displaystyle \sum _{elements}{k}_{el}}F={\displaystyle \sum _{nodes}{f}_n}+{\displaystyle \sum _{elements}{r}_{el}}$

where U is a vector of unknown nodal displacements and rotations (as components in the  $\{i,j,k\}$ basis), F is a vector of nodal forces and moments, and

$\begin{array}{l}{k}_{el}={\displaystyle \underset{A_{el}}{\overset{}{\int }}{B}_{\kappa }^T{D}_{\kappa }{B}_{\kappa }+B_{\theta }^{T}T{B}_{\theta }dA}\approx {\displaystyle \frac{A_{el}}{3}{\displaystyle \sum _i{B}_{\kappa }^T({\xi }_i)D_{\kappa }{B}_{\kappa }({\xi }_i)}+B_{\theta }^T({\xi }_i)T{B}_{\theta }({\xi }_i)}\\ r_{el}={\displaystyle \underset{A_{el}}{\overset{}{\int }}{N}^T{p}_{3}dA}\approx {\displaystyle \frac{A_{el}}{3}{\displaystyle \sum _i{N}^T({\xi }_i)p_3}}\end{array}$

are the element stiffness matrix and force vector, with ${\xi }_i$ the coordinates for a 3 point Gauss integration scheme

$\begin{array}{l}{\xi }_1=0.5(x_1{e}_1+y_1{e}_2)+0.5(x_2{e}_1+y_2{e}_2)\\ {\xi }_2=0.5(x_2{e}_1+y_2{e}_2)+0.5(x_3{e}_1+y_3{e}_2)\\ {\xi }_3=0.5(x_1{e}_1+y_1{e}_2)+0.5(x_3{e}_1+y_3{e}_2)\end{array}$

 

 

Finite element method for out-of-plane deformation of elastic Reissner-Mindlin plates: The principle of virtual work for the Reissner-Mindlin plate consists of two coupled virtual work equations for the unknown displacement and rotation of transverse cross-sections of the plate.    The out of plane displacement and the rotation of the plate are now two are independent variables,

 

1. There is no need to ensure that the rotation are related by a compatibility constraint;

 

2. The displacement and rotation must be continuous across element boundaries, but there is no need for continuity of the derivative of the displacement.

 

At first sight this suggests that displacements and rotations could be interpolated using the standard 2D shape functions, but unfortunately these all suffer from shear locking if the plate thickness is small.   There have been many attempts to design locking resistant plate elements.   We will not attempt to review these here, and instead describe the simple ‘Min3’ triangular element developed by Tessler and Hughes as a representative example.

 

To proceed with the finite element formulation, it is helpful to combine the two virtual work equations for the plate beam into a single, symmetric form.   For this purpose we define a shear/axial strain vector

$\gamma =\nabla u_3-e_3\times \theta$

Adding the two virtual work equations then yields

${\displaystyle \underset{A}{\overset{}{\int }}M\cdot \delta \kappa }dA+{\displaystyle \underset{A}{\overset{}{\int }}V\cdot \delta \gamma dA}+{\displaystyle \underset{A}{\int }\left(\nabla \delta u_3\right)\cdot T\cdot (\nabla u_3)dA}={\displaystyle \underset{C}{\int }{P}_3}\delta u_{3}ds+{\displaystyle \underset{A}{\int }{p}_3\delta u_3}dA+{\displaystyle \underset{C}{\int }Q\cdot }\delta \theta ds$

 

As for the Kirchhoff plate, the elements are triangular, such as the example illustrated in Fig 8.54.   A generic element has three nodes with coordinates $x_k=(x_k{e}_1+y_k{e}_2)$, where k=1,2,3.  Each node has an unknown out of plane displacement $u_3^{(k)}$ , and two rotation components ${\theta }^{(k)}={\theta }_{\alpha }^{(k)}{e}_{\alpha }$. 

 

The displacement and rotation of the plate are now interpolated separately as

${\left[u_3,{\theta }_1,{\theta }_2\right]}^T=N{\left[\begin{array}{ccccccc}{u}_3^1& {\theta }_1^1& {\theta }_2^1& \cdots & u_3^3& {\theta }_1^3& {\theta }_2^3\end{array}\right]}^T$

where

$N=\left[\begin{array}{ccccccc}{L}_1& N_{\theta 1}^1& N_{\theta 2}^1& & L_3& N_{\theta 1}^3& N_{\theta 2}^3\\ 0& L_1& 0& \cdots & 0& L_3& 0\\ 0& 0& L_1& & 0& 0& L_3\end{array}\right]$

with $L_i$ the areal coordinates

$\begin{array}{l}{L}_1=\left[(x_2-x)(y_3-y)-(x_3-x)(y_2-y)\right]/(2A_{el})\\ L_2=\left[(x_3-x)(y_1-y)-(x_1-x)(y_3-y)\right]/(2A_{el})\\ L_3=\left[(x_1-x)(y_2-y)-(x_2-x)(y_1-y)\right]/(2A_{el})\\ 2A_{el}=(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\end{array}$

and

$\begin{array}{c}{N}_{\theta 1}^1={\displaystyle \frac{1}{2}\left(b_2{L}_3{L}_1-b_3{L}_1{L}_2\right)}\\ N_{\theta 1}^2={\displaystyle \frac{1}{2}\left(b_3{L}_1{L}_2-b_1{L}_2{L}_3\right)}\\ N_{\theta 1}^3={\displaystyle \frac{1}{2}\left(b_1{L}_2{L}_3-b_2{L}_3{L}_1\right)}\end{array}\begin{array}{c}{N}_{\theta 2}^1={\displaystyle \frac{1}{2}\left(c_2{L}_3{L}_1-c_3{L}_1{L}_2\right)}\\ N_{\theta 2}^2={\displaystyle \frac{1}{2}\left(c_3{L}_1{L}_2-c_1{L}_2{L}_3\right)}\\ N_{\theta 2}^3={\displaystyle \frac{1}{2}\left(c_1{L}_2{L}_3-c_2{L}_3{L}_1\right)}\end{array}$

with

$\begin{array}{c}{b}_1=y_2-y_3{b}_2=y_3-y_1{b}_3=y_1-y_2\\ c_1=x_3-x_2{c}_2=x_1-x_3{c}_3=x_2-x_1\end{array}$

The derivatives of these shape functions can be calculated as

${\displaystyle \frac{d{N}_{\theta }}{dx}={\displaystyle \frac{d{N}_{\theta }}{dL}{\displaystyle \frac{dL}{dx}}}}$

where

$\begin{array}{l}{\displaystyle \frac{\partial N_{\theta }}{\partial L}=\left[\begin{array}{ccc}\partial N_{\theta 1}^1/\partial L_1& \partial N_{\theta 1}^1/\partial L_2& \partial N_{\theta 1}^1/\partial L_3\\ \partial N_{\theta 2}^1/\partial L_2& N_{\theta 2}^1/\partial L_2& \\ ⋮& & \end{array}\right]}\\ ={\displaystyle \frac{1}{2}\left[\begin{array}{ccc}{b}_2{L}_3-b_3{L}_2& -b_3{L}_1& b_2{L}_1\\ c_2{L}_3-c_3{L}_2& -c_3{L}_1& c_2{L}_1\\ b_3{L}_2& b_3{L}_1-b_1{L}_3& -b_1{L}_2\\ c_3{L}_2& c_3{L}_1-c_1{L}_3& -c_1{L}_2\\ -b_2{L}_3& b_1{L}_3& b_1{L}_2-b_2{L}_1\\ -c_2{L}_3& c_1{L}_3& c_1{L}_2-c_2{L}_1\end{array}\right]}\end{array}$

$\begin{array}{l}{\displaystyle \frac{dL}{dx}=\left[\begin{array}{cc}\partial L_1/\partial x_1& \partial L_1/\partial x_2\\ \partial L_2/\partial x_1& \partial L_2/\partial x_2\\ \partial L_3/\partial x_1& \partial L_3/\partial x_2\end{array}\right]={\displaystyle \frac{1}{2A_{el}}\left[\begin{array}{cc}{b}_1& c_1\\ b_2& c_2\\ b_3& c_3\end{array}\right]}}\\ 2A_{el}=c_3{b}_2-b_3{c}_2\end{array}$

 

The displacement gradient, shear strain and curvatures are related to the nodal displacements and rotations by

 

$\begin{array}{l}{\left[\begin{array}{cc}\partial u_3/\partial x_1& \partial u_3/\partial x_2\end{array}\right]}^T=B_{\nabla u}\left[\begin{array}{ccccccccc}{u}_3^1& {\theta }_1^1& {\theta }_2^1& u_3^2& {\theta }_1^2& {\theta }_2^2& u_3^3& {\theta }_1^3& {\theta }_2^3\end{array}\right]\\ {\left[\begin{array}{cc}{\gamma }_1& {\gamma }_2\end{array}\right]}^T=B_{\gamma }{\left[\begin{array}{ccccccccc}{u}_3^1& {\theta }_1^1& {\theta }_2^1& u_3^2& {\theta }_1^2& {\theta }_2^2& u_3^3& {\theta }_1^3& {\theta }_2^3\end{array}\right]}^T\\ {\left[\begin{array}{ccc}{\kappa }_{11}& {\kappa }_{22}& {\kappa }_{12}+{\kappa }_{21}\end{array}\right]}^T=B_{\kappa }{\left[\begin{array}{ccccccccc}{u}_3^1& {\theta }_1^1& {\theta }_2^1& u_3^2& {\theta }_1^2& {\theta }_2^2& u_3^3& {\theta }_1^3& {\theta }_2^3\end{array}\right]}^T\end{array}$

$B_{\nabla u}=\left[\begin{array}{ccccccccc}{\displaystyle \frac{\partial L_1}{\partial x_1}}& {\displaystyle \frac{\partial N_{\theta 1}^1}{\partial x_1}}& {\displaystyle \frac{\partial N_{\theta 2}^1}{\partial x_1}}& {\displaystyle \frac{\partial L_2}{\partial x_1}}& {\displaystyle \frac{\partial N_{\theta 1}^2}{\partial x_1}}& {\displaystyle \frac{\partial N_{\theta 2}^2}{\partial x_1}}& {\displaystyle \frac{\partial L_3}{\partial x_1}}& {\displaystyle \frac{\partial N_{\theta 1}^3}{\partial x_1}}& {\displaystyle \frac{\partial N_{\theta 2}^3}{\partial x_1}}\\ {\displaystyle \frac{\partial L_1}{\partial x_2}}& {\displaystyle \frac{\partial N_{\theta 1}^1}{\partial x_2}}& {\displaystyle \frac{\partial N_{\theta 2}^1}{\partial x_2}}& {\displaystyle \frac{\partial L_2}{\partial x_2}}& {\displaystyle \frac{\partial N_{\theta 1}^2}{\partial x_2}}& {\displaystyle \frac{\partial N_{\theta 2}^2}{\partial x_2}}& {\displaystyle \frac{\partial L_3}{\partial x_2}}& {\displaystyle \frac{\partial N_{\theta 1}^3}{\partial x_2}}& {\displaystyle \frac{\partial N_{\theta 2}^3}{\partial x_2}}\end{array}\right]$

$B_{\gamma }=B_{\nabla u}+\left[\begin{array}{ccccccccc}0& 0& L_1& 0& 0& L_2& 0& 0& L_3\\ 0& -L_1& 0& 0& -L_2& 0& 0& -L_3& 0\end{array}\right]$

 

The internal forces follow as

${\left[\begin{array}{ccc}{M}_{11}& M_{22}& M_{12}\end{array}\right]}^T=D_{\kappa }{\left[\begin{array}{ccc}{\kappa }_{11}& {\kappa }_{22}& {\kappa }_{12}+{\kappa }_{21}\end{array}\right]}^T$

${\left[\begin{array}{cc}{V}_1& V_2\end{array}\right]}^T=D_{\gamma }{\left[\begin{array}{cc}{\gamma }_1& {\gamma }_2\end{array}\right]}^T$

$D_{\kappa }={\displaystyle \frac{E{h}^3}{12(1-{\nu }^2)}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& (1-\nu )/2\end{array}\right]D_{\gamma }=2\beta h\mu \left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]}$

 

With these definitions, the principle of virtual work reduces to a system of linear equations

$KU=FK={\displaystyle \sum _{elements}{k}_{\kappa }+{\beta }_{*}{k}_{\gamma }+k_T}F={\displaystyle \sum _{nodes}{f}_n}+{\displaystyle \sum _{elements}{r}_{el}}$

where U is a vector of unknown nodal displacements and rotations, F is a vector of nodal forces and moments, and

$\begin{array}{l}{k}_{\kappa }={\displaystyle \underset{A_{el}}{\overset{}{\int }}{B}_{\kappa }^T{D}_{\kappa }{B}_{\kappa }dA}\approx {\displaystyle \frac{A_{el}}{3}{\displaystyle \sum _i{B}_{\kappa }^T({\xi }_i)D_{\kappa }{B}_{\kappa }({\xi }_i)}}\\ k_{\gamma }={\displaystyle \underset{A_{el}}{\overset{}{\int }}{B}_{\gamma }^T{D}_{\gamma }{B}_{\gamma }dA}\approx {\displaystyle \frac{A_{el}}{3}{\displaystyle \sum _i{B}_{\gamma }^T({\xi }_i)D_{\gamma }{B}_{\gamma }({\xi }_i)}}\\ k_T={\displaystyle \underset{A_{el}}{\overset{}{\int }}{B}_{\nabla u}^{T}T{B}_{\nabla u}dA}\approx {\displaystyle \frac{A_{el}}{3}{\displaystyle \sum _i{B}_{\nabla u}^T({\xi }_i)T{B}_{\nabla u}({\xi }_i)}}\\ r_{el}={\displaystyle \underset{A_{el}}{\overset{}{\int }}{N}^{T}pdA}\approx {\displaystyle \frac{A_{el}}{3}{\displaystyle \sum _i{N}^T}({\xi }_i)p}\end{array}$

are the element stiffness matrix and force vector.  The integrals have been evaluated with a 3-point Gaussian quadrature scheme, with

$\begin{array}{l}{\xi }_1=0.5(x_1{e}_1+y_1{e}_2)+0.5(x_2{e}_1+y_2{e}_2)\\ {\xi }_2=0.5(x_2{e}_1+y_2{e}_2)+0.5(x_3{e}_1+y_3{e}_2)\\ {\xi }_3=0.5(x_1{e}_1+y_1{e}_2)+0.5(x_3{e}_1+y_3{e}_2)\end{array}$

A 1 point scheme will work as well.  The coefficient ${\beta }_{*}$ is a numerical correction that improves the performance of the element.   Tessler and Hughes report that using

${\beta }_{*}={\displaystyle \frac{2{\displaystyle {\sum }_{i=2,3,5,6,8,9}{k}_{ii}^{(\kappa )}}}{2{\displaystyle {\sum }_{i=2,3,5,6,8,9}{k}_{ii}^{(\kappa )}}+{\displaystyle {\sum }_{i=2,3,5,6,8,9}{k}_{ii}^{(\gamma )}}}}$

 

gives good results.

 

Example: A simple demonstration of both the Kirchhoff and Reissner-Mindlin elements are shown in the figure.  Fig. (a) compares the predictions of the Kirchhoff element for the deflection and rotation of a simply supported circular plate with radius R, thickness h, Young’s modulus E, and Poisson’s ratio $\nu$ subjected to uniform pressure p with the exact solution (the mesh is shown inset).  The FEA predictions are within 0.1% of the exact result.   Fig. (b) shows the predicted deflection of circular Reissner-Mindlin plate, for several values of the normalized ratio of bending to shear stiffness $E{h}^2/(\beta \mu R^2)$.   The shear deformation increases the deflection slightly for $E{h}^2/(\beta \mu R^2)$ >0.1 but the effect of shearing for the circular plate is less pronounced than the corresponding behavior of a cantilever beam discussed in the preceding section.

 

The codes that plot these graphs are provided in the files named FEM_plate_Kirchhoff.m and FEM_plate_Mindlin.m at

https://github.com/albower/Applied_Mechanics_of_Solids/