4.1. Euler-Bernoulli beam elements

Typically two beam formulations are addressed in literature, the Euler-Bernoulli beam and the Timoshenko beam. The Euler-Bernoulli beam is valid for relatively slender beams, whereas the Timoshenko beam is also applicable where shear deformations are present. In this chapter only in-plane, linear beams are considered. All forces and moments are in-plane. First, we focus on the Euler-Bernoulli beam formulation, which assumes shear deformations are negligible.

Equilibrium of a beam

The bending moment \(M\) in a beam is defined as

(4.1)\[ M = \int_{-h/2}^{h/2} \sigma_{xx} y \, dy \]

and the shear force \(V\) defined by

(4.2)\[ V = \int_{-h/2}^{h/2} \sigma_{xy} \, dy \]

Considering the vertical equilibrium of a beam element, it follows that

(4.3)\[ \int_{d\Omega} V n \, d\Gamma + \int_{\Omega} f_y \, d\Omega = 0 \]

where, for a 1D domain of \(x\in[a,b]\), we have \(\int_{d\Omega} V n \, d\Gamma = V|_{x=b} - V|_{x=a}\). Applying divergence theorem, this time to go from a boundary integral to a volume integral, gives

(4.4)\[ \int_{\Omega} V_{,x} \, d\Omega + \int_{\Omega} f_y \, d\Omega = 0 \]

Since the equilibrium described by Equation (4.4) must hold for an infinitely small segment of a beam, it must hold that

(4.5)\[ V_{,x} + f_y = 0 \]

For rotational equlibrium, it follows that

(4.6)\[ \int_{d\Omega} M n \, d\Gamma - \int_{d\Omega} V n x \, d\Gamma - \int_{\Omega} f_y x \, d\Omega = 0 \]

which can be rearranged such that

(4.7)\[ \int_{\Omega} M_{,x} \, d\Omega - \int_{\Omega} V \, d\Omega - \int_{\Omega} V_{,x} x \, d\Omega - \int_{\Omega} f_y x \, d\Omega = 0 \]

Since satisfaction of the translational equilibrium already ensures that \(\int_{\Omega} (V_{,x} + f_y) x \, d\Omega = 0\), the rotational equilibrium implies that

(4.8)\[ M_{,x} - V = 0 \]

…

Derivation of the strong form

The Euler-Bernoulli beam does not allow for shear deformation. As a result, the rotation of planes perpendicular to the beam axis can be directly related to the displacement \(w\):

(4.9)\[ \theta = \frac{\text{d}w}{\text{d}x} \]

which means that

(4.10)\[ \kappa = \frac{\text{d}^2w}{\text{d}x^2} \]

and the bending moment \(M\) in the beam is determined by:

(4.11)\[ M=-EI \cdot \kappa=-EI \cdot \frac{\text{d}^2w}{\text{d}x^2} \]

where \(EI\) is the bending stiffness of the beam.

Taking the derivative of all terms in Equation (4.8) to \(x\) and substituting in (4.5) yields

(4.12)\[ \frac{\partial^2 M}{\partial x^2} + f_y = 0 \]

Assuming \(EI\) to be constant and substituting in () yields

(4.13)\[ -EI \frac{\partial^4 w}{\partial x^4} + f_y = 0 \]

which is the strong form equation of equilibrium for an Euler-Bernoulli beam. Being a fourth-order equation, two boundary conditions are needed at both ends of the beam. Dirichlet boundary conditions involved the prescription of the displacement or rotation and Neumann involves either the shear force or moment. With appropriate boundary conditions, the boundary value problem is complete and can be solved.

Weak governing equation

Following procedures from (REF to weak form eq.), the weak form of equilibrium for a beam can be developed. Multiplying Equation (4.12) by a weight function \(\bar{w}\), from an appropriately defined space, which is equal to zero where Dirichlet boundary conditions are applied and integration over the beam \(\Omega\) yields:

(4.14)\[ \int_{\Omega} \bar{w} M_{,xx} \, d\Omega + \int_{\Omega} \bar{w} f_y \, d\Omega = 0 \]

Integrating by parts the term involving the moment \(M\) yields:

(4.15)\[ - \int_{\Omega} \bar{w}_{,x} M_{,x} \, d\Omega + \int_{\Gamma} \bar{w} M_{,x} n \, d\Gamma + \int_{\Omega} \bar{w} f_y \, d\Omega = 0 \]

Applying integration by parts again, this time to the term \(\int_{\Omega} \bar{w}_{,x} M_{,x} \, d\Omega\), yields

(4.16)\[ \int_{\Omega} \bar{w}_{,xx} M \, d\Omega - \int_{\Gamma} \bar{w}_{,x} M n \, d\Gamma + \int_{\Gamma} \bar{w} M_{,x} n \, d\Gamma + \int_{\Omega} \bar{w} f_y \, d\Omega = 0 \]

Inserting now the Neumann boundary conditions and the constitutive relation from Equation (), solving the governing weak form equation for a beam involves: find \(w \in S\) such that

(4.17)\[ -\int_{\Omega} \bar{w}_{,xx} EI w_{,xx} \, d\Omega - \int_{\Gamma_M} \bar{w}_{,x} T \, d\Gamma + \int_{\Gamma_V} \bar{w} f_y \, d\Omega + \int_{\Omega} \bar{w} f_y \, d\Omega = 0 \]

where \(S\) and \(V\) are appropriately defined spaces.

Discrete form

A Galerkin problem for a Euler-Bernoulli beam involves: find \(w^h \in S^h\) such that:

(4.18)\[ \int_{\Omega} \bar{w}_{,xx}^{h} EI w_{,xx}^h \, d\Omega = - \int_{\Gamma_M} \bar{w}_{,x}^h T \, d\Gamma + \int_{\Gamma_V} \bar{w}^h F_y \, d\Gamma + \int_{\Omega} \bar{w}^h f_y \, d\Omega = 0 \]

where \(S^h \subset S\) and \(V^h \subset V\) are finite-dimensional spaces. As for continuum elements, we wish to express the displacement field \(w\) in terms of shape functions and nodal degrees of freedom. The problem that arises is that standard \(C^0\) finite element shape functions are not suitable. The above equation requires the evaluation of the second derivative of the interpolated displacement field. However the second derivative of a \(C^0\) continuous function does not exist in a classical sense. The use of \(C^0\) interpolations for fourth-order problems is not mathematically consistent. Crucially, convergence of the solution is not assured for interpolations with insufficient continuity.

The solution is to use \(C^1\) shape functions. Such shape functions can be constructed relatively easily in one dimension (the extension to multiple dimensions is however far from trivial). Hermitian polynomials are \(C^1\) functions, and involve both displacement and rotational degrees of freedom. A Hermitian beam element with two nodes has four degrees of freedom (two displacement degrees of freedom and two rotation degrees of freedom). This results in a cubic interpolation of the displacement along the element. The displacement field \(w^h\) is given by:

(4.19)\[ w^h (x) = \sum_{i}^{2} \left(N_i (x) w_i + M_i (x) \phi_i\right) \]

where \(w_i\) and \(\phi_i\) are degrees of freedom associated with node \(i\), and the shape functions are equal to:

(4.20)\[ N_1 = \frac{-(x-b)^2 \left(-h+2\left(a-x\right)\right)}{h^3} \]

(4.21)\[ N_2 = \frac{(x-a)^2 (h+2(b - x))}{h^3} \]

(4.22)\[ M_1 = \frac{(x-a)(x-b)^2}{h^2} \]

(4.23)\[ M_2 = \frac{(x-a)^2 (x-b)}{h^2} \]

for an element of length \(h\). The displacement at a point in the beam is given by:

(4.24)\[\begin{split} w^h = \mathbf{N} \mathbf{a}_e = \ \begin{bmatrix} N_1 & M_1 & N_2 & M_2 \end{bmatrix} \begin{Bmatrix} w_1 \\ \phi_1 \\ w_2 \\ \phi_2 \end{Bmatrix} \end{split}\]

It is necessary to compute both the first and second derivatives of \(w\) with respect to \(x\). The first derivative is given by:

(4.25)\[\begin{split} w^h_{,x} = \mathbf{N}_{,x} \mathbf{a}_e = \ \begin{bmatrix} N_{1,x} & M_{1,x} & N_{2,x} & M_{2,x} \end{bmatrix} \begin{Bmatrix} w_1 \\ \phi_1 \\ w_2 \\ \phi_2 \end{Bmatrix} \end{split}\]

and the second derivative can be calculated through:

(4.26)\[\begin{split} w^h_{,xx} = \mathbf{N}_{,xx} \mathbf{a}_e = \ \begin{bmatrix} N_{1,xx} & M_{1,xx} & N_{2,xx} & M_{2,xx} \end{bmatrix} \begin{Bmatrix} w_1 \\ \phi_1 \\ w_2 \\ \phi_2 \end{Bmatrix} \end{split}\]

Now that \(w^h\), \(w^h_{,x}\) and \(w^h_{,xx}\) can be computed given \(\mathbf{a}_e\), they can be inserted into the Galerkin problem,

(4.27)\[ \int_{\Omega} (\mathbf{N}_{,xx} \mathbf{c}_e)^T EI \mathbf{N}_{,xx} \mathbf{a}_e \, d\Omega = -\int_{\Gamma_M} (\mathbf{N}_{,x} \mathbf{c}_e)^T T \, d\Gamma + \int_{\Gamma_V} (\mathbf{N} \mathbf{c}_e)^T F_y \, d\Gamma + \int_{\Omega} (\mathbf{N} \mathbf{c}_e)^T f_y \, d\Omega \]

After some rearranging,

(4.28)\[ \int_{\Omega} \mathbf{N}_{,xx}^T EI \mathbf{N}_{,xx} \, d\Omega \mathbf{a}_e = - \int_{\Gamma_M} \mathbf{N}_{,x}^T T \, d\Gamma + \int_{\Gamma_V} \mathbf{N}^T F_y \, d\Gamma + \int_{\Omega} \mathbf{N}^T f_y \, d\Omega \]

The stiffness matrix is then given by:

(4.29)\[ \mathbf{K} = \int_{\Omega} \mathbf{N}_{,xx}^T EI \mathbf{N}_{,xx} \, d\Omega \]

and the RHS vector is given by:

(4.30)\[ \mathbf{f} = \int_{\Gamma_V} \mathbf{N}^T F_y \, d\Gamma - \int_{\Gamma_M} \mathbf{N}^T T \, d\Gamma + \int_{\Omega} \mathbf{N}^T f_y \, d\Omega \]

The operation to form the RHS vector essentially translates the applied loads into equivalent nodal shear forces and moments. Taking derivatives of the Hermitian shape functions in (4.20) until (4.23):

(4.31)\[ N_{1,x} = \frac{2(x-b)(3x+h-2a-b)}{h^3} \]

(4.32)\[ N_{1,xx} = \frac{2(6x+h-2a-4b)}{h^3} \]

(4.33)\[ M_{1,x} = \frac{(x-b)(3x-2a-b)}{h^2} \]

(4.34)\[ M_{1,xx} = \frac{2(3x-a-2b)}{h^2} \]

(4.35)\[ N_{2,x} = - \frac{2(x-a)(3x-h-a-2b)}{h^3} \]

(4.36)\[ N_{2,xx} = - \frac{2(6x-h-4a-2b)}{h^3} \]

(4.37)\[ M_{2,x} = \frac{(x-a)(3x-a-2b)}{h^2} \]

(4.38)\[ M_{2,xx} = \frac{2(3x-2a-b)}{h^2} \]

Assuming the centre of the element is at \(x=0\) (\(-a = b = h/2\)), the above equations can be simplified significantly. Considering just the second derivatives with respect to \(x\),

(4.39)\[ N_{1,xx} = \frac{12x}{h^3} \]

(4.40)\[ M_{1,xx} = \frac{6x}{h^2} - \frac{1}{h} \]

(4.41)\[ N_{2,xx} = - \frac{12x}{h^3} \]

(4.42)\[ M_{2,xx} = \frac{6x}{h^2} + \frac{1}{h} \]

Inserting these terms into equation (4.29), the stiffness matrix is of the form:

(4.43)\[\begin{split} \mathbf{K}_e = \int_{-h/2}^{h/2} \begin{bmatrix} \frac{12x}{h^3} \\ \frac{6x}{h^2} - \frac{1}{h} \\ - \frac{12x}{h^3} \\ \frac{6x}{h^2} + \frac{1}{h} \end{bmatrix} EI \begin{bmatrix} \frac{12x}{h^3} & \frac{6x}{h^2} - \frac{1}{h} & - \frac{12x}{h^3} & \frac{6x}{h^2} + \frac{1}{h} \end{bmatrix} dx. \end{split}\]

Integrating the terms in the stiffness matrix exactly from \(-h/2\) to \(h/2\) (and assuming EI to be constant), the element stiffness matrix is equal to:

\[\begin{split} \mathbf{K}_e = \begin{bmatrix} \frac{12EI}{h^3} & \frac{6EI}{h^2} & -\frac{12EI}{h^3} & \frac{6EI}{h^2}\\ \frac{6EI}{h^2} & \frac{4EI}{h} & -\frac{6EI}{h^2} & \frac{2EI}{h}\\ -\frac{12EI}{h^3} & -\frac{6EI}{h^2} & \frac{12EI}{h^3} & -\frac{6EI}{h^2}\\ \frac{6EI}{h^2} & \frac{2EI}{h} & -\frac{6EI}{h^2} & \frac{4EI}{h} \end{bmatrix} \end{split}\]

Note that this final results is most easily achieved with the assumption that the element is centered around \(x=0\), but that the same stiffness matrix is obtained for any arbitrarily positioned element of length \(h\).