Euler-Bernoulli beam elements

Euler–Bernoulli beam theory is a simplification of the linear theory of elasticity and covers the case for small deflections of a beam that is subjected to lateral loads only.

The following topics are discussed:

Strain measures
Internal virtual work rate Jacobian
First variations of strains
Second variations of strains
Interpolation
Integration

Strain measures

The internal virtual work rate associated with axial stress is

δ W_{1}^{I} = \int_{L^{f}} \int_{A} σ^{f} δ ε^{f} d A d L^{f},

where $σ^{f}$ and $δ ε^{f}$ are any material stress and strain measures associated with axial deformation at the point $(S, g, h)$ of the beam, since strains are assumed to be small. For this purpose we will use Green's strain so that

ε^{f} = \frac{1}{2} [{(λ^{f})}^{2} - 1],

where ${(λ^{f})}^{2} = {(d l^{f})}^{2} / {(d L^{f})}^{2}$ , the square of the ratio of current configuration length to reference configuration length in the axial direction on the fiber. From Equation 1 and its equivalent in the reference configuration, we have

\begin{matrix} ε^{f} d L^{f} = \frac{1}{2} & [(\frac{d x}{d S} \cdot \frac{d x}{d S} + 2 g \frac{d x}{d S} \cdot \frac{d n_{1}}{d S} + 2 h \frac{d x}{d S} \cdot \frac{d n_{2}}{d S}) \\ \cdot {(\frac{d x}{d S} \cdot \frac{d x}{d S} + 2 g \frac{d X}{d S} \cdot \frac{d N_{1}}{d S} + 2 h \frac{d x}{d S} \cdot \frac{d N_{2}}{d S})}^{- 1} - 1] \\ \cdot {(\frac{d X}{d S} \cdot \frac{d X}{d S} + 2 g \frac{d X}{d S} \cdot \frac{d N_{1}}{d S} + 2 h \frac{d X}{d S} \cdot \frac{d N_{2}}{d S})}^{\frac{1}{2}} d S . \end{matrix}

Again, neglecting all but first-order terms in g and h because of the slenderness assumption, this becomes

\begin{matrix} ε^{f} d L^{f} \approx [ε & + g G^{- 1} (\frac{d x}{d S} \cdot \frac{d n_{1}}{d S} - (1 + ε) \frac{d X}{d S} \cdot \frac{d N_{1}}{d S}) \\ + h G^{- 1} (\frac{d x}{d S} \cdot \frac{d n_{2}}{d S} - (1 + ε) \frac{d X}{d S} \cdot \frac{d N_{2}}{d S})] d L, \end{matrix}

where

G = \frac{d X}{d S} \cdot \frac{d X}{d S} and ε = \frac{1}{2} (λ^{2} - 1) = \frac{1}{2} [{(\frac{d l}{d L})}^{2} - 1] is the Green’s strain of the beam
  axis.

This simplification allows us to write the internal virtual work rate associated with axial stress as

δ W_{1}^{I} = \int_{L^{f}} \int_{A} σ^{f} δ ε^{f} d A d L^{f} = \int_{L} \int_{A} σ^{f} δ {\tilde{ε}}^{f} d A d L,

where

{\tilde{ε}}^{f} = ε - g K_{2} + h K_{1},

with

\begin{matrix} ε & = \frac{1}{2} (λ^{2} - 1) as Green’s strain of the beam
  axis, \\ K_{2} & = - G^{- 1} (\frac{d x}{d S} \cdot \frac{d n_{1}}{d S} - (1 + ε) \frac{d X}{d S} \cdot \frac{d N_{1}}{d S}) \end{matrix}

and

K_{1} = - G^{- 1} (\frac{d x}{d S} \cdot \frac{d n_{2}}{d S} - (1 + ε) \frac{d X}{d S} \cdot \frac{d N_{2}}{d S}) .

Now the cross-sectional base vectors $n_{1}$ and $n_{2}$ are assumed to remain normal to the beam axis, so

\frac{d x}{d S} \cdot n_{1} = \frac{d x}{d S} \cdot n_{2} = 0.

Hence,

\frac{d x}{d S} \cdot \frac{d n_{1}}{d S} = - \frac{d^{2} x}{d S^{2}} \cdot n_{1}, and \frac{d x}{d S} \cdot \frac{d n_{2}}{d S} = - \frac{d^{2} x}{d S^{2}} \cdot n_{2} .

So we have

K_{2} = G^{- 1} (n_{1} \cdot \frac{d^{2} x}{d S^{2}} - (1 + ε) N_{1} \cdot \frac{d^{2} X}{d S^{2}})

and

K_{1} = - G^{- 1} (n_{2} \cdot \frac{d^{2} x}{d S^{2}} - (1 + ε) N_{2} \cdot \frac{d^{2} X}{d S^{2}}) .

This defines the generalized strains associated with axial stretch. For torsional strain the internal virtual work rate is

δ W_{2}^{I} = \int_{L} \int_{A} γ_{α}^{t} τ_{α}^{f} δ e_{1} d A d L,

where

e_{1} = n_{2} \cdot \frac{d n_{1}}{d S} - N_{2} \cdot \frac{d N_{1}}{d S}

and $γ_{α}^{t}$ is the proportionality constant between shear strains and torsional strains (see Beam element formulation for details). For computational simplicity the form of the torsional strains is taken to be

e_{1} = n_{2} \cdot \frac{d {\tilde{n}}_{1}}{d S} - N_{2} \cdot \frac{d N_{1}}{d S},

where

{\tilde{n}}_{1} = C (p_{1} ω^{1} + p_{2} ω^{2}) \cdot N_{1},

and

\begin{matrix} p_{1} & = \frac{1}{2} (1 - g), \\ p_{2} & = \frac{1}{2} (1 + g), \\ - 1 & \leq g \leq 1. \end{matrix}

This assumes a linear interpolation of rotation $ω$ along the beam. Thus, the generalized strain measures for these beams are

$ε$: axial strain,
$K_{1}$ and $K_{2}$: the beam curvature change measures, and
$e_{1}$: the torsional strain.

With these measures, the internal virtual work rate can be written

δ W^{I} = \int_{L} [δ ε \int_{A} σ^{f} d A - δ K_{2} \int_{A} g σ^{f} d A + δ K_{1} \int_{A} h σ^{f} d A + δ e_{1} \int_{A} γ_{α}^{t} τ_{α}^{f} d A] d L .

Internal virtual work rate Jacobian

For the Jacobian matrix of the overall Newton method, Equation 3, the variation of this internal work rate with respect to nodal displacement variations must be formed. The constitutive theory is written as

{\begin{matrix} d σ^{f} \\ d τ^{f} \end{matrix}} = [H] {\begin{matrix} d ε^{f} \\ d γ^{f} \end{matrix}},

and so

d δ W^{I} = \int_{L} [⌊ δ ε δ K_{1} δ K_{2} δ e_{1} ⌋ [B] {\begin{matrix} d ε \\ d K_{1} \\ d K_{2} \\ d e_{1} \end{matrix}} + d δ ε \int_{A} σ^{f} d A + d δ K_{1} \int_{A} h σ^{f} d A

- d δ K_{2} \int_{A} g σ^{f} d A + d δ e_{1} \int_{A} γ_{α}^{f} τ_{α}^{f} d A] d L,

where $[B]$ is the section stiffness matrix obtained by integration over the cross-section. See Beam element formulation for more information on section integration.

First variations of strains

Taking the variations of the above strain definitions gives directly

\begin{matrix} δ ε & = G^{- 1} \frac{d x}{d S} \cdot \frac{d δ u}{d S}, \\ δ K_{2} & = G^{- 1} (δ n_{1} \cdot \frac{d^{2} x}{d S^{2}} + n_{1} \cdot \frac{d^{2} δ u}{d S^{2}} - δ ε N_{1} \cdot \frac{d^{2} X}{d S^{2}}) \\ δ K_{1} & = G^{- 1} (δ n_{2} \cdot \frac{d^{2} x}{d S^{2}} + n_{2} \cdot \frac{d^{2} δ u}{d S^{2}} - δ ε N_{2} \cdot \frac{d^{2} X}{d S^{2}}) \\ and \\ δ e_{1} & = δ n_{2} \cdot \frac{d {\tilde{n}}_{1}}{d S} + n_{2} \cdot \frac{d δ {\tilde{n}}_{1}}{d S} . \end{matrix}

In these expressions we need $δ n_{1}$ and $δ n_{2}$ as well as $δ {\tilde{n}}_{1}$ ; these are now derived. From the expression for ${\tilde{n}}_{1}$ , namely

{\tilde{n}}_{1} = C (p_{1} ω^{1} + p_{2} ω^{2}) \cdot N_{1},

another ancillary vector ${\bar{n}}_{1}$ , normal to the tangent, is defined by

{\bar{n}}_{1} = {\tilde{n}}_{1} - {\tilde{n}}_{1} \cdot tt

so that

n_{1} = {\bar{n}}_{1} / | {\bar{n}}_{1} | .

In addition,

t = g^{- \frac{1}{2}} \frac{d x}{d S}

so that

δ t = g^{- \frac{1}{2}} (I - tt) \cdot \frac{d δ u}{d S} .

Since $(t_{1} n_{1} n_{2})$ form an orthonormal triad,

\begin{matrix} (2) & \begin{matrix} δ n_{1} & = δ n_{1} \cdot (n_{2} n_{2} + tt) \\ = δ n_{1} \cdot n_{2} n_{2} + δ n_{1} \cdot tt \\ = δ n_{1} \cdot n_{2} n_{2} - g^{- \frac{1}{2}} n_{1} \cdot \frac{d δ u}{d S} t, \end{matrix} \end{matrix}

because $n_{1} \cdot t = 0$ . From the definition of $n_{1}$ , it is straightforward to show that

δ n_{1} = \frac{1}{| {\bar{n}}_{1} |} (I - n_{1} n_{1}) \cdot δ \bar{n} .

\begin{matrix} δ n_{1} \cdot n_{2} & = \frac{1}{| {\bar{n}}_{1} |} n_{2} \cdot δ n_{1} \\ = \frac{1}{| {\bar{n}}_{1} |} n_{2} \cdot [δ {\tilde{n}}_{1} - (δ {\tilde{n}}_{1} \cdot t) t - ({\tilde{n}}_{1} \cdot t) δ t] \\ = \frac{1}{| {\bar{n}}_{1} |} [n_{2} \cdot δ {\tilde{n}}_{1} - g^{- \frac{1}{2}} ({\tilde{n}}_{1} \cdot t) n_{2} \cdot \frac{d δ u}{d S}] \end{matrix}

and

δ {\tilde{n}}_{1} = δ ω \times {\tilde{n}}_{1} with ω interpolated
  linearly.

Thus,

δ n_{1} \cdot n_{2} = \frac{1}{| {\bar{n}}_{1} |} (δ ω \cdot {\tilde{n}}_{1} \times n_{2} - g^{- \frac{1}{2}} {\tilde{n}}_{1} \cdot {tn}_{2} \cdot \frac{d δ u}{d S}) .

We can also write

{\tilde{n}}_{1} \times n_{2} = {\tilde{n}}_{1} t - {\tilde{n}}_{1} \cdot t n_{1}

and

| {\bar{n}}_{1} | = {\bar{n}}_{1} \cdot n_{1} = {\tilde{n}}_{1} \cdot n_{1} .

Combining terms appropriately,

δ n_{1} \cdot n_{2} = δ ω \cdot t - \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} [δ ω \cdot n_{1} + g^{- \frac{1}{2}} n_{2} \cdot \frac{d δ u}{d S}] .

Hence, from Equation 2

δ n_{1} = δ ω \cdot {tn}_{2} - \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} [δ ω \cdot n_{1} + g^{- \frac{1}{2}} n_{2} \cdot \frac{d δ u}{d S}] n_{2} - g^{- \frac{1}{2}} n_{1} \cdot \frac{d δ u}{d S} t .

In a similar manner one can show that

δ n_{2} = - δ ω \cdot {tn}_{1} - \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} [δ ω \cdot n_{1} + g^{- \frac{1}{2}} n_{1} \cdot \frac{d δ u}{d S}] n_{1} - g^{- \frac{1}{2}} n_{2} \cdot \frac{d δ u}{d S} t .

The first variations of strain become

\begin{matrix} δ ε & = G^{- 1} \frac{d x}{d S} \cdot \frac{d δ u}{d S} \\ δ K_{2} & = G^{- 1} [n_{2} \cdot \frac{d^{2} x}{d S^{2}} (t - \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} n_{1}) \cdot δ ω \\ - g^{- \frac{1}{2}} (t \cdot \frac{d^{2} x}{d S^{2}} n_{1} + \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} n_{2} \cdot \frac{d^{2} x}{d S^{2}}) n_{2} \cdot \frac{d δ u}{d S} \\ + n_{1} \frac{d^{2} δ u}{d S^{2}} - δ ε N_{1} \cdot \frac{d^{2} X}{d S^{2}}] \\ δ K_{1} & = G^{- 1} [n_{1} \cdot \frac{d^{2} x}{d S^{2}} (t - \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} n_{1}) \cdot δ ω \\ + g^{- \frac{1}{2}} (t \cdot \frac{d^{2} x}{d S^{2}} n_{1} + \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} n_{1} \cdot \frac{d^{2} x}{d S^{2}}) n_{2} \cdot \frac{d δ u}{d S} \\ - n_{2} \frac{d^{2} δ u}{d S^{2}} - δ ε N_{2} \cdot \frac{d^{2} X}{d S^{2}}] \\ and \\ δ e_{1} & = (\frac{d {\tilde{n}}_{1}}{d S} \cdot n_{1} \frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} - \frac{d {\tilde{n}}_{1}}{d S} t) n_{1} \cdot δ ω \\ + g^{- \frac{1}{2}} [\frac{{\tilde{n}}_{1} \cdot t}{{\tilde{n}}_{1} \cdot n_{1}} \frac{d {\tilde{n}}_{1}}{d S} \cdot n_{1} - \frac{d {\tilde{n}}_{1}}{d S} \cdot t] n_{2} \cdot \frac{d δ u}{d S} \\ + {\tilde{n}}_{1} \times n_{2} \cdot \frac{d δ ω}{d S} \end{matrix}

so that

n_{1} = {\bar{n}}_{1} / | {\bar{n}}_{1} | .

In addition,

t = g^{- \frac{1}{2}} \frac{d x}{d S}

so that

δ t = g^{- \frac{1}{2}} (I - tt) \cdot \frac{d δ u}{d S},

since $(t_{1} n_{1} n_{2})$ form an orthonormal triad,

\begin{matrix} δ n_{1} & = δ n_{1} \cdot (n_{2} n_{2} + tt) \\ = δ n_{1} \cdot n_{2} n_{2} + δ n_{1} \cdot tt \\ = δ n_{1} \cdot n_{2} n_{2} - g^{- \frac{1}{2}} n_{1} \cdot \frac{d δ u}{d S} t, \end{matrix}

because

n_{1} \cdot t = 0.

Second variations of strains

The second variation of the axial strain is simply

d δ ε = G^{- 1} \frac{d δ u}{d S} \cdot \frac{d d u}{d S} .

To compute the second variations of bending strain, we need expressions for $d δ n_{1}, d δ n_{2}$ . These are obtained by approximating

\begin{matrix} δ {\tilde{n}}_{1} & \approx n_{2} t \cdot δ ω - t n_{1} \cdot \frac{d δ u}{d S} \\ δ {\tilde{n}}_{2} & \approx - n_{1} t \cdot δ ω - t n_{1} \cdot \frac{d δ u}{d S} \\ δ \tilde{t} & \approx δ ω \cdot n_{2} - δ ω \cdot n_{1} n_{2} \end{matrix}

from which

d δ {\tilde{n}}_{1} = d n_{2} t \cdot δ ω + n_{2} d t \cdot δ ω - d t n_{1} \cdot \frac{d δ u}{d S} - t d n_{1} \cdot \frac{d δ u}{d S} .

Using these expressions, the second variations of bending strains are written as

\begin{matrix} d δ K_{1} = & G^{- 1} [n_{2} \cdot \frac{d^{2} x}{d S^{2}} t \cdot d ω t \cdot δ ω \\ - t \cdot \frac{d^{2} x}{d S^{2}} (n_{1} \cdot \frac{d d u}{d S} t \cdot δ ω + n_{1} \cdot \frac{d δ u}{d S} t \cdot d ω) \\ + n_{1} \cdot \frac{d^{2} d u}{d S^{2}} t \cdot δ ω + t \cdot \frac{d^{2} d u}{d S^{2}} n_{2} \cdot \frac{d δ u}{d S} \\ + n_{1} \cdot \frac{d^{2} δ u}{d S^{2}} t \cdot d ω + t \cdot \frac{d^{2} δ u}{d S^{2}} n_{2} \cdot \frac{d d u}{d S} \\ + d δ ε N_{2} \cdot \frac{d^{2} X}{d S^{2}} \\ + n_{1} \cdot \frac{d^{2} x}{d S^{2}} (n_{2} \cdot d ω n_{1} \cdot δ ω - n_{1} \cdot d ω n_{2} \cdot δ ω + n_{2} \cdot \frac{d δ u}{d S} n_{2} \cdot d ω) \\ - n_{2} \cdot \frac{d^{2} x}{d S^{2}} (n_{2} \cdot \frac{d δ u}{d S} n_{1} \cdot d ω) \\ - t \cdot \frac{d^{2} x}{d S^{2}} (n_{2} \cdot \frac{d d u}{d S} n_{2} t \cdot \frac{d δ u}{d S})] \end{matrix}

and

\begin{matrix} d δ K_{2} = G^{- 1} & [- n_{1} \cdot \frac{d^{2} x}{d S^{2}} t \cdot δ ω t \cdot d ω \\ - t \cdot \frac{d^{2} x}{d S^{2}} (n_{2} \cdot \frac{d δ u}{d S} t \cdot d ω + n_{2} \cdot \frac{d d u}{d S} t \cdot δ ω) \\ + n_{2} \cdot \frac{d^{2} d u}{d S^{2}} t \cdot δ ω - t \cdot \frac{d^{2} d u}{d S^{2}} n_{1} \cdot \frac{d δ u}{d S} \\ + n_{2} \cdot \frac{d^{2} δ u}{d S^{2}} t \cdot d ω - t \cdot \frac{d^{2} δ u}{d S^{2}} n_{1} \cdot \frac{d d u}{d S} \\ - d δ ε N_{1} \cdot \frac{d^{2} X}{d S^{2}} \\ + n_{2} \cdot \frac{d^{2} x}{d S^{2}} (n_{2} \cdot d ω n_{1} \cdot δ ω - n_{1} \cdot d ω n_{2} \cdot δ ω + n_{1} \cdot \frac{d δ u}{d S} n_{1} \cdot d ω) \\ - n_{1} \cdot \frac{d^{2} x}{d S^{2}} (n_{1} \cdot \frac{d δ u}{d S} n_{2} \cdot d ω) \\ - t \cdot \frac{d^{2} x}{d S^{2}} (n_{1} \cdot \frac{d d u}{d S} n_{1} t \cdot \frac{d δ u}{d S})] . \end{matrix}

For the torsional strain contribution to the initial stress matrix, we approximate

\begin{matrix} δ e_{2} & \approx t \cdot \frac{d δ ω}{d S} \\ d δ e_{1} & \approx d t \cdot \frac{d δ ω}{d S} = d ω \times t \cdot \frac{d δ ω}{d S} . \end{matrix}

This matrix is again unsymmetric.

Interpolation

In the virtual work equation the strains include second derivatives of displacement. For this reason continuity of rotation as well as displacement is needed so that the Hermitian polynomial interpolation functions are the minimum interpolation order needed. These are used here. The Hermite cubic is written in terms of the function value and its derivative at the ends of the interval:

u (g) = N_{1} (g) u^{1} + N_{2} (g) u^{2} + N_{3} (g) \frac{d u^{1}}{d g} + N_{4} (g) \frac{d u^{2}}{d g}, - 1 < g < 1,

with node 1 at $g = - 1$ and node 2 at $g = + 1$ .

These functions are used in Abaqus to interpolate the components of displacement $u$ and the initial position vector $x$ , so that the elements are basically isoparametric. In addition, rotation of $(n_{1} n_{2})$ about the beam axis, $ϕ_{a}$ , is interpolated linearly. This interpolation is unsatisfactory for the user, because the nodal variables are

(u_{x} u_{y} u_{z} \frac{d u_{x}}{d S} \frac{d u_{y}}{d S} \frac{d u_{z}}{d S} ϕ_{a}) .

The last four of these variables are difficult to work with; furthermore, making them the same in all elements sharing the same node causes excessive constraint of axial stretch in these elements, especially if the beam axis is not continuous through the node, as in a frame structure or “T” junction.

To avoid this difficulty, the following procedure is adopted. At a node the tangent to the beam axis is

t = \frac{d x}{d l} = \frac{d x}{d S} \frac{d S}{d l},

\frac{d x}{d S} = \frac{d l}{d S} t .

Now suppose we store as degrees of freedom at the node,

(u_{x}, u_{y}, u_{z}, ω_{x}, ω_{y}, ω_{z}, d l / d S),

where $ω$ is the rotation definition introduced above. Since the initial geometry and hence $T$ , the initial direction of the beam axis, is known, $t = C \cdot T,$ where $C$ is the rotation matrix defined by $ω$ , and hence

\begin{matrix} (3) & \frac{d x}{d S} = \frac{d l}{d S} C \cdot T \end{matrix}

is defined by these variables and the initial geometry. Furthermore, $ϕ_{a}$ is directly available from $ω$ and $t$ . Thus, the above set is a satisfactory set of nodal variables. To eliminate the unwanted axial strain constraint, in Abaqus the stretch $d l / d S$ at the node of each such element is taken as an internal variable, local to the element (a third internal node is created for this purpose, and so it is not shared with neighboring elements.)

It should be remarked that the above transformation (Equation 3) is nonlinear. This leads to some complications—for example, the d'Alembert forces no longer have the simple form $M_{M N} {\ddot{\bar{u}}}^{N};$ rather, a matrix ${\bar{M}}_{M N} = Q^{M P} M_{P Q} Q^{Q N}$ replaces $M_{M N}$ where $Q^{M P}$ and $Q^{Q N}$ use the transformation (Equation 3) and its appropriate time derivatives. The resulting Jacobian is nonsymmetric; Abaqus ignores the nonsymmetric terms.

Integration

The cross-section integration has already been discussed in the context of the shear beams—it is the same for these beams. Along the beam axis, the integration schemes are as described below.

Stiffness and internal forces: Three Gauss points are used. Two Gauss points are not sufficient because the torsional strain is not independent.
Mass and distributed loads: Three Gauss points are used. Rotary inertia is not included in the mass, except for rotation of the section about the beam axis, where it is included to avoid singularity in perfectly straight beams.