8.1 2D Euler–Bernoulli Beam Element

Updated Lagrangian formulation

At load step $i$ , the known current configuration becomes the reference configuration for the next load step. The unknown target configuration is reached at load step $i+1$ . Displacements and stresses at $i+1$ are evaluated with respect to the preceding configuration at $i$ , and all derivatives and integrals are evaluated with respect to that preceding configuration [1].

StrainSmall

MaterialLinear elastic

DisplacementLarge

ElementStraight, prismatic, two-node

Beam modelEuler–Bernoulli

ReferenceKnown current configuration

Updated Lagrangian description

The configuration at load step $i$ is known and serves as the reference configuration.
The configuration at load step $i+1$ is the unknown target configuration.
The nodal displacement increments and stress increments for the current step are measured from configuration $i$ .

Configurations and coordinate frames

The global reference frame $(x,y)$ remains fixed for every load step. The known beam configuration at load step $i$ has the local frame $(\bar{x},\bar{y})$ . The next configuration at $i+1$ is described relative to this frame.

The local quantities $\bar u(\bar x)$ and $\bar v(\bar x)$ shown between configurations $i$ and $i+1$ are displacement increments for the current load step. They are not total displacements measured from the initial undeformed configuration. Because these increments are small, the local frames at the beginning and end of one load step are close to one another.

Updated Lagrangian beam element in the known current configuration and the next target configuration — Figure 1. Initial undeformed position, known current configuration at load step $i$ , and target configuration at load step $i+1$ . The barred displacements are increments for the current step.

Deformation energy and strain increments

The deformation energy of the complete structure in the target configuration is:

U_{i+1}=\frac{1}{2}\sum\int_L\left(EA\,\varepsilon_{0,i+1}^{2}+EI\,\kappa_{i+1}^{2}\right)\,d\bar{x}

Using the known generalized strains at load step $i$ and the unknown increments for the current load step,

U_{i+1}=\frac{1}{2}\sum\int_L\left[EA\left(\varepsilon_{0,i}+\Delta\varepsilon_0\right)^2+EI\left(\kappa_i+\Delta\kappa\right)^2\right]d\bar{x}

The summation sign denotes the assembly process, while the generalized strains satisfy

\left\{\begin{aligned}\varepsilon_{0,i+1}&=\varepsilon_{0,i}+\Delta\varepsilon_0\\[5pt]\kappa_{i+1}&=\kappa_i+\Delta\kappa\end{aligned}\right.

During load step $i+1$ , the previously accumulated quantities $\varepsilon_{0,i}$ and $\kappa_i$ are known constants. They therefore do not depend on the element displacement-increment vector $\Delta\boldsymbol u_{el}$ .

\begin{aligned}\frac{\partial\varepsilon_{0,i+1}}{\partial\Delta\boldsymbol u_{el}}&=\frac{\partial\Delta\varepsilon_0}{\partial\Delta\boldsymbol u_{el}},&\qquad \frac{\partial^2\varepsilon_{0,i+1}}{\partial\Delta\boldsymbol u_{el}^{2}}&=\frac{\partial^2\Delta\varepsilon_0}{\partial\Delta\boldsymbol u_{el}^{2}}\\[7pt]\frac{\partial\kappa_{i+1}}{\partial\Delta\boldsymbol u_{el}}&=\frac{\partial\Delta\kappa}{\partial\Delta\boldsymbol u_{el}},&\qquad \frac{\partial^2\kappa_{i+1}}{\partial\Delta\boldsymbol u_{el}^{2}}&=\frac{\partial^2\Delta\kappa}{\partial\Delta\boldsymbol u_{el}^{2}}\end{aligned}

Element residual and tangent stiffness

The resulting expressions have the same form as those developed for the Euler–Bernoulli beam element in Chapter 6, with the total generalized strains and nodal displacements replaced by their current-step increments.

The element internal force vector and tangent stiffness matrix are evaluated with respect to the known reference configuration at load step $i$ . The known reference strain state is represented explicitly by $\varepsilon_0$ and $\kappa$ .

\boldsymbol f_{el}=\frac{\partial U_{el}}{\partial\Delta\boldsymbol u_{el}}=\int_L\left(EA\,\varepsilon_0\frac{\partial\Delta\varepsilon_0}{\partial\Delta\boldsymbol u_{el}}+EI\,\kappa\frac{\partial\Delta\kappa}{\partial\Delta\boldsymbol u_{el}}\right)d\bar{x}

\boldsymbol k_{el}=\frac{\partial\boldsymbol f_{el}}{\partial\Delta\boldsymbol u_{el}}=\frac{\partial^2U_{el}}{\partial\Delta\boldsymbol u_{el}^{2}}=\int_L\left[EA\,\varepsilon_0\frac{\partial^2\Delta\varepsilon_0}{\partial\Delta\boldsymbol u_{el}^{2}}+EA\frac{\partial\Delta\varepsilon_0}{\partial\Delta\boldsymbol u_{el}}\left(\frac{\partial\Delta\varepsilon_0}{\partial\Delta\boldsymbol u_{el}}\right)^T+EI\,\kappa\frac{\partial^2\Delta\kappa}{\partial\Delta\boldsymbol u_{el}^{2}}+EI\frac{\partial\Delta\kappa}{\partial\Delta\boldsymbol u_{el}}\left(\frac{\partial\Delta\kappa}{\partial\Delta\boldsymbol u_{el}}\right)^T\right]d\bar{x}

The nodal unknowns are the displacement increments $\Delta u$ and $\Delta v$ between load steps $i$ and $i+1$ . This distinction is the essential change from the formulation in Chapter 6.

Configuration and strain updates

After convergence of each load step, the nodal geometry is updated:

\left\{\begin{aligned}\boldsymbol{x}_{i+1}&=\boldsymbol{x}_i+\Delta\boldsymbol{u}\\[5pt]\boldsymbol{y}_{i+1}&=\boldsymbol{y}_i+\Delta\boldsymbol{v}\end{aligned}\right.

The generalized strain measures are accumulated consistently:

\varepsilon_0:=\varepsilon_0+\Delta\varepsilon_0

\kappa:=\kappa+\Delta\kappa

This Updated Lagrangian procedure introduces an incremental approximation because the kinematics are written for the displacement increment within each load step. The approximation error becomes small when these increments are small. Increasing the number of load steps therefore reduces the incremental error while also helping the Newton–Raphson iterations follow the equilibrium path.

MATLAB implementation

The nonlinear solution follows the same Newton–Raphson structure as Chapter 6, but the displacement vector is reset at the start of each load step because it represents the current increment. The converged increment is then added to the stored total displacement.

main.mIncremental displacement update

for istep=1:nstep
    S=zeros(neq,1);
    ...
    updt
    St(:,istep+1)=St(:,istep)+S;
end

In stiff.m, the axial strain and curvature from the preceding converged configuration are added to the current increments after evaluating the incremental kinematics:

stiff.mGeneralized-strain accumulation

deriv
ep=ep+strn(ii,1,istep);
ca=ca+strn(ii,2,istep);

The line strn(ii,1:2,istep+1)=[ep ca]; stores the updated generalized strains for configuration $i+1$ . These values become the known reference strain state for the following load step.

The program updt.m is essential to the Updated Lagrangian formulation. After each converged load step it updates the nodal coordinates, element lengths, element orientations, and rotation matrices. The resulting current geometry becomes the reference configuration for the next step.

Example 1 — moment-loaded cantilever

A cantilever of length $L_b=100\ \mathrm{mm}$ has a rectangular cross-section $5\ \mathrm{mm}\times1\ \mathrm{mm}$ (width × thickness), Young's modulus $E=2\times10^5\ \mathrm{MPa}$ , and an applied free-end moment $M=5236\ \mathrm{Nmm}$ . The loading is divided into 100 equidistant load steps.

The analytical maximum rotation is:

\varphi_{\max}=\frac{ML_b}{EI}=\frac{5236\cdot100}{2\times10^5\cdot5\cdot1^3/12}=2\pi\ \mathrm{rad}

The figure shows five of the 100 load steps. A very large total rotation is obtained through many small displacement increments. The verified MATLAB result at the final step is $\varphi_{\max}=6.26691\ \mathrm{rad}$ , close to the analytical value $2\pi=6.28319\ \mathrm{rad}$ .

Moment-loaded cantilever and five deformed configurations approaching a full rotation — Figure 2. Moment-loaded cantilever and five of the 100 equidistant load steps. The plotted coordinates are normalized as $x/L_b$ and $y/L_b$ .

Example 2 — Mattiasson benchmark

This benchmark is taken from Mattiasson [2], where reference solutions for large-deflection beam and frame problems are obtained using elliptic integrals. The two symmetric loading cases and the measured response components are shown below.

Two symmetric large-deflection beam-frame benchmark configurations — Figure 3. Mattiasson beam-frame benchmark: geometry, loading, and displacement measures used in the comparison [2].

Symmetry is used when plotting the deformed structure. The Updated Lagrangian numerical results show very good agreement with the elliptic-integral reference solution from [2].

Symmetric deformed configurations and normalized response curves compared with Mattiasson — Figure 4. Symmetric deformed configurations and normalized load–displacement curves. Solid curves are obtained with the Updated Lagrangian beam model; dashed curves are the elliptic-integral reference values from [2].

Relation to Section 8.2

These two examples will also be considered in Section 8.2 using a two-node isoparametric Timoshenko beam element. Within the adopted assumptions, that formulation does not introduce the same incremental approximation associated with displacement-increment size and may therefore reach the converged solution with very few load steps. The present Euler–Bernoulli Updated Lagrangian algorithm instead requires sufficiently small load-step increments to keep the incremental approximation error small.

MATLAB package

The downloadable package contains the complete program for Example 1, including the nonlinear load-step loop, incremental strain evaluation, geometry update, result printing, and README instructions.

MATLAB verification

Executed successfully in MATLAB 26.1 (R2026a) Update 4. All 100 load steps reached convergence; the final step required four Newton–Raphson iterations, with a convergence measure of $4.61\times10^{-8}$ . The computed maximum nodal rotation was $6.26691\ \mathrm{rad}$ . The code uses standard MATLAB syntax and is expected to be compatible with recent MATLAB versions. No additional toolbox is required.

MUpdated Lagrangian Euler–Bernoulli MATLAB package6 MATLAB files and README · tested in R2026aDownload ZIP

References

E. Chatzi and K. Agathos, Linearized weak form and total Lagrangian formulation of a bar element, Institute of Structural Engineering, ETH Zurich, Lecture 3, 3 October 2019.
K. Mattiasson, Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals, International Journal for Numerical Methods in Engineering, 16, 145–153, 1981.