FEA TipsFinite Element Theory & MATLAB Programs
HomeChapter 12 · Total Lagrangian Formulation (TL)12.1

Section 12.1 of Chapter 12: Total Lagrangian Formulation (TL)

CST Finite Element

A complete Total Lagrangian formulation of the plane-stress constant-strain triangle for small strains and large displacements, including the consistent tangent stiffness matrix and two verified nonlinear examples.

Total LagrangianPlane stressCSTMATLABANSYS

01

Formulation and hypotheses

StrainSmall
DisplacementLarge
MaterialLinear elastic
Stress statePlane stress

In the Total Lagrangian formulation, the initial undeformed geometry is the reference configuration. All derivatives and integrals are evaluated with respect to that initial configuration, and the Green–Lagrange strain measure is used.

Because the strains and stresses are constant within each plane-stress CST finite element, the deformation energy is:

U=12Ael=1nel(∫Aelh εTσ dA)U=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\int_{A_{el}}h\,\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}\,dA\right)
U=12Ael=1nel(εTσ Vel)U=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}\,V_{el}\right)

Here, Ael=1nel\mathcal{A}_{el=1}^{n_{el}} denotes the standard finite element assembly operation from the first through the last element. The signed element area is:

Ael=12[(x2−x1)(y3−y1)−(x3−x1)(y2−y1)]A_{el}=\frac12\left[(x_2-x_1)(y_3-y_1)-(x_3-x_1)(y_2-y_1)\right]

The sign of AelA_{el} depends on whether the element nodes are numbered counterclockwise or clockwise. The physical area is ∣Ael∣|A_{el}|, and the element volume is:

Vel=h∣Ael∣V_{el}=h|A_{el}|

hh is the element thickness, and neln_{el} is the number of finite elements.

02

Strain and stress measures

ε={εxεyγxy}\boldsymbol{\varepsilon}=\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{Bmatrix}
σ={σxσyτxy}\boldsymbol{\sigma}=\begin{Bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{Bmatrix}

The engineering shear component satisfies γxy=2Exy\gamma_{xy}=2E_{xy}. For a linear elastic material in plane stress:

σ=Dε\boldsymbol{\sigma}=\boldsymbol{D}\boldsymbol{\varepsilon}
D=E1−ν2[1ν0ν10001−ν2]\boldsymbol{D}=\frac{E}{1-\nu^2}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&\dfrac{1-\nu}{2}\end{bmatrix}

The two-dimensional Green–Lagrange strain components are:

εx=∂u∂x+12(∂u∂x)2+12(∂v∂x)2\varepsilon_x=\frac{\partial u}{\partial x}+\frac12\left(\frac{\partial u}{\partial x}\right)^2+\frac12\left(\frac{\partial v}{\partial x}\right)^2
εy=∂v∂y+12(∂u∂y)2+12(∂v∂y)2\varepsilon_y=\frac{\partial v}{\partial y}+\frac12\left(\frac{\partial u}{\partial y}\right)^2+\frac12\left(\frac{\partial v}{\partial y}\right)^2
γxy=∂u∂y+∂v∂x+∂u∂x∂u∂y+∂v∂x∂v∂y\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\frac{\partial u}{\partial x}\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\frac{\partial v}{\partial y}

03

CST interpolation

The displacement field is interpolated with the same linear functions used in the CST element of Section 9.1:

{u(x,y)v(x,y)}=N(x,y)uel\begin{Bmatrix}u(x,y)\\v(x,y)\end{Bmatrix}=\boldsymbol{N}(x,y)\boldsymbol{u}_{el}
N(x,y)=[NuNv]=[N10N20N300N10N20N3]\boldsymbol{N}(x,y)=\begin{bmatrix}\boldsymbol{N}_u\\\boldsymbol{N}_v\end{bmatrix}=\begin{bmatrix}N_1&0&N_2&0&N_3&0\\0&N_1&0&N_2&0&N_3\end{bmatrix}
N1(x,y)=12Ael[(y2−y3)x−(x2−x3)y+x2y3−x3y2]N_1(x,y)=\frac{1}{2A_{el}}\left[(y_2-y_3)x-(x_2-x_3)y+x_2y_3-x_3y_2\right]
N2(x,y)=12Ael[(y3−y1)x−(x3−x1)y+x3y1−x1y3]N_2(x,y)=\frac{1}{2A_{el}}\left[(y_3-y_1)x-(x_3-x_1)y+x_3y_1-x_1y_3\right]
N3(x,y)=12Ael[(y1−y2)x−(x1−x2)y+x1y2−x2y1]N_3(x,y)=\frac{1}{2A_{el}}\left[(y_1-y_2)x-(x_1-x_2)y+x_1y_2-x_2y_1\right]

x1,y1,x2,y2,x3,y3x_1,y_1,x_2,y_2,x_3,y_3 are the nodal coordinates. Their derivatives give:

∂u∂x=∂Nu∂xuel=Buxuel,∂u∂y=∂Nu∂yuel=Buyuel\frac{\partial u}{\partial x}=\frac{\partial\boldsymbol{N}_u}{\partial x}\boldsymbol{u}_{el}=\boldsymbol{B}_{u_x}\boldsymbol{u}_{el},\qquad\frac{\partial u}{\partial y}=\frac{\partial\boldsymbol{N}_u}{\partial y}\boldsymbol{u}_{el}=\boldsymbol{B}_{u_y}\boldsymbol{u}_{el}
∂v∂x=∂Nv∂xuel=Bvxuel,∂v∂y=∂Nv∂yuel=Bvyuel\frac{\partial v}{\partial x}=\frac{\partial\boldsymbol{N}_v}{\partial x}\boldsymbol{u}_{el}=\boldsymbol{B}_{v_x}\boldsymbol{u}_{el},\qquad\frac{\partial v}{\partial y}=\frac{\partial\boldsymbol{N}_v}{\partial y}\boldsymbol{u}_{el}=\boldsymbol{B}_{v_y}\boldsymbol{u}_{el}
Bux=∂Nu∂x=12Ael[y2−y30y3−y10y1−y20]\boldsymbol{B}_{u_x}=\frac{\partial\boldsymbol{N}_u}{\partial x}=\frac{1}{2A_{el}}\begin{bmatrix}y_2-y_3&0&y_3-y_1&0&y_1-y_2&0\end{bmatrix}
Buy=∂Nu∂y=12Ael[x3−x20x1−x30x2−x10]\boldsymbol{B}_{u_y}=\frac{\partial\boldsymbol{N}_u}{\partial y}=\frac{1}{2A_{el}}\begin{bmatrix}x_3-x_2&0&x_1-x_3&0&x_2-x_1&0\end{bmatrix}
Bvx=∂Nv∂x=12Ael[0y2−y30y3−y10y1−y2]\boldsymbol{B}_{v_x}=\frac{\partial\boldsymbol{N}_v}{\partial x}=\frac{1}{2A_{el}}\begin{bmatrix}0&y_2-y_3&0&y_3-y_1&0&y_1-y_2\end{bmatrix}
Bvy=∂Nv∂y=12Ael[0x3−x20x1−x30x2−x1]\boldsymbol{B}_{v_y}=\frac{\partial\boldsymbol{N}_v}{\partial y}=\frac{1}{2A_{el}}\begin{bmatrix}0&x_3-x_2&0&x_1-x_3&0&x_2-x_1\end{bmatrix}
uel=[u1v1u2v2u3v3]T\boldsymbol{u}_{el}=\begin{bmatrix}u_1&v_1&u_2&v_2&u_3&v_3\end{bmatrix}^{T}

04

Green–Lagrange strains

After introducing the CST interpolation, the Green–Lagrange strains become:

εx=Buxuel+12uelTGxuel\varepsilon_x=\boldsymbol{B}_{u_x}\boldsymbol{u}_{el}+\frac12\boldsymbol{u}_{el}^{T}\boldsymbol{G}_x\boldsymbol{u}_{el}
εy=Bvyuel+12uelTGyuel\varepsilon_y=\boldsymbol{B}_{v_y}\boldsymbol{u}_{el}+\frac12\boldsymbol{u}_{el}^{T}\boldsymbol{G}_y\boldsymbol{u}_{el}
γxy=(Buy+Bvx)uel+12uelTGxyuel\gamma_{xy}=(\boldsymbol{B}_{u_y}+\boldsymbol{B}_{v_x})\boldsymbol{u}_{el}+\frac12\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{xy}\boldsymbol{u}_{el}

The matrices Gx\boldsymbol{G}_x, Gy\boldsymbol{G}_y, and Gxy\boldsymbol{G}_{xy} are symmetric 6×66\times6 matrices:

Gx=BuxTBux+BvxTBvx\boldsymbol{G}_x=\boldsymbol{B}_{u_x}^{T}\boldsymbol{B}_{u_x}+\boldsymbol{B}_{v_x}^{T}\boldsymbol{B}_{v_x}
Gy=BuyTBuy+BvyTBvy\boldsymbol{G}_y=\boldsymbol{B}_{u_y}^{T}\boldsymbol{B}_{u_y}+\boldsymbol{B}_{v_y}^{T}\boldsymbol{B}_{v_y}
Gxy=BuxTBuy+BuyTBux+BvxTBvy+BvyTBvx\boldsymbol{G}_{xy}=\boldsymbol{B}_{u_x}^{T}\boldsymbol{B}_{u_y}+\boldsymbol{B}_{u_y}^{T}\boldsymbol{B}_{u_x}+\boldsymbol{B}_{v_x}^{T}\boldsymbol{B}_{v_y}+\boldsymbol{B}_{v_y}^{T}\boldsymbol{B}_{v_x}

05

Strain-displacement matrix

The matrix relating the strain increment to the nodal-displacement increment is defined by:

dε=B(uel) duel,B=∂ε∂ueld\boldsymbol{\varepsilon}=\boldsymbol{B}(\boldsymbol{u}_{el})\,d\boldsymbol{u}_{el},\qquad\boldsymbol{B}=\frac{\partial\boldsymbol{\varepsilon}}{\partial\boldsymbol{u}_{el}}
{dεxdεydγxy}=[Bux+uelTGxBvy+uelTGyBuy+Bvx+uelTGxy]⏟Bduel\begin{Bmatrix}d\varepsilon_x\\d\varepsilon_y\\d\gamma_{xy}\end{Bmatrix}=\underbrace{\begin{bmatrix}\boldsymbol{B}_{u_x}+\boldsymbol{u}_{el}^{T}\boldsymbol{G}_x\\\boldsymbol{B}_{v_y}+\boldsymbol{u}_{el}^{T}\boldsymbol{G}_y\\\boldsymbol{B}_{u_y}+\boldsymbol{B}_{v_x}+\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{xy}\end{bmatrix}}_{\boldsymbol{B}}d\boldsymbol{u}_{el}

The matrix is split into a constant linear part and a displacement-dependent nonlinear part:

B(uel)=B0+BL(uel)\boldsymbol{B}(\boldsymbol{u}_{el})=\boldsymbol{B}_0+\boldsymbol{B}_L(\boldsymbol{u}_{el})
[Bux+uelTGxBvy+uelTGyBuy+Bvx+uelTGxy]⏟B=[BuxBvyBuy+Bvx]⏟B0+[uelTGxuelTGyuelTGxy]⏟BL\underbrace{\begin{bmatrix}\boldsymbol{B}_{u_x}+\boldsymbol{u}_{el}^{T}\boldsymbol{G}_x\\\boldsymbol{B}_{v_y}+\boldsymbol{u}_{el}^{T}\boldsymbol{G}_y\\\boldsymbol{B}_{u_y}+\boldsymbol{B}_{v_x}+\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{xy}\end{bmatrix}}_{\boldsymbol{B}}=\underbrace{\begin{bmatrix}\boldsymbol{B}_{u_x}\\\boldsymbol{B}_{v_y}\\\boldsymbol{B}_{u_y}+\boldsymbol{B}_{v_x}\end{bmatrix}}_{\boldsymbol{B}_0}+\underbrace{\begin{bmatrix}\boldsymbol{u}_{el}^{T}\boldsymbol{G}_x\\\boldsymbol{u}_{el}^{T}\boldsymbol{G}_y\\\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{xy}\end{bmatrix}}_{\boldsymbol{B}_L}

06

Potential energy and equilibrium

The total potential energy is:

Π=U+V=12Ael=1nel(εTσVel)−uTF\Pi=U+V=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}V_{el}\right)-\boldsymbol{u}^{T}\boldsymbol{F}

Here, V=−uTFV=-\boldsymbol{u}^{T}\boldsymbol{F} is the potential of the external loads, u\boldsymbol{u} is the global nodal displacement vector, and F\boldsymbol{F} is the global external force vector. The virtual-work equation is:

δΠ=Ael=1nel(δεTσVel)−δuTF=0\delta\Pi=\mathcal{A}_{el=1}^{n_{el}}\left(\delta\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}V_{el}\right)-\delta\boldsymbol{u}^{T}\boldsymbol{F}=0
δΠ=Ael=1nel(δuelTBTσVel)−δuTF=0\delta\Pi=\mathcal{A}_{el=1}^{n_{el}}\left(\delta\boldsymbol{u}_{el}^{T}\boldsymbol{B}^{T}\boldsymbol{\sigma}V_{el}\right)-\delta\boldsymbol{u}^{T}\boldsymbol{F}=0

This relation holds for every geometrically admissible virtual displacement field. The resulting nonlinear equilibrium system is:

Ψ(u)=Ael=1nel(BTσVel)−F=0\boldsymbol{\Psi}(\boldsymbol{u})=\mathcal{A}_{el=1}^{n_{el}}\left(\boldsymbol{B}^{T}\boldsymbol{\sigma}V_{el}\right)-\boldsymbol{F}=\boldsymbol{0}
Ψ(u)=Ael=1nel[(B0T+BLT(uel))σVel]−F=0\boldsymbol{\Psi}(\boldsymbol{u})=\mathcal{A}_{el=1}^{n_{el}}\left[\left(\boldsymbol{B}_0^{T}+\boldsymbol{B}_L^{T}(\boldsymbol{u}_{el})\right)\boldsymbol{\sigma}V_{el}\right]-\boldsymbol{F}=\boldsymbol{0}

07

Tangent stiffness matrix

The first variation of the equilibrium residual is:

δΨ=Ael=1nel(δBTσVel)+Ael=1nel(BTδσVel)=KTδu\delta\boldsymbol{\Psi}=\mathcal{A}_{el=1}^{n_{el}}\left(\delta\boldsymbol{B}^{T}\boldsymbol{\sigma}V_{el}\right)+\mathcal{A}_{el=1}^{n_{el}}\left(\boldsymbol{B}^{T}\delta\boldsymbol{\sigma}V_{el}\right)=\boldsymbol{K}_T\delta\boldsymbol{u}

Because B0\boldsymbol{B}_0 is constant:

δB=δ(B0+BL)=δBL\delta\boldsymbol{B}=\delta(\boldsymbol{B}_0+\boldsymbol{B}_L)=\delta\boldsymbol{B}_L

For one CST finite element, the variation of the nonlinear part gives the geometric or stress-dependent contribution:

δBTσVel=δBLTσVel=δ([[uelTGx]T[uelTGy]T[uelTGxy]T]⏟6×3)σVel=∫V(σxGx+σyGy+τxyGxy)δuelVel\delta\boldsymbol{B}^{T}\boldsymbol{\sigma}V_{el}=\delta\boldsymbol{B}_{L}^{T}\boldsymbol{\sigma}V_{el}=\delta\left(\underbrace{\begin{bmatrix}\left[\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{x}\right]^{T}\\\left[\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{y}\right]^{T}\\\left[\boldsymbol{u}_{el}^{T}\boldsymbol{G}_{xy}\right]^{T}\end{bmatrix}}_{6\times3}\right)\boldsymbol{\sigma}V_{el}=\int_{V}\left(\sigma_x\boldsymbol{G}_x+\sigma_y\boldsymbol{G}_y+\tau_{xy}\boldsymbol{G}_{xy}\right)\delta\boldsymbol{u}_{el}V_{el}
δBTσVel=(σxGx+σyGy+τxyGxy)δuelVel\delta\boldsymbol{B}^{T}\boldsymbol{\sigma}V_{el}=\left(\sigma_x\boldsymbol{G}_x+\sigma_y\boldsymbol{G}_y+\tau_{xy}\boldsymbol{G}_{xy}\right)\delta\boldsymbol{u}_{el}V_{el}

Therefore, the global tangent stiffness matrix is:

KT=Ael=1nel[Vel(σxGx+σyGy+τxyGxy)]+Ael=1nel(VelBTDB)\boldsymbol{K}_T=\mathcal{A}_{el=1}^{n_{el}}\left[V_{el}\left(\sigma_x\boldsymbol{G}_x+\sigma_y\boldsymbol{G}_y+\tau_{xy}\boldsymbol{G}_{xy}\right)\right]+\mathcal{A}_{el=1}^{n_{el}}\left(V_{el}\boldsymbol{B}^{T}\boldsymbol{D}\boldsymbol{B}\right)

The tangent stiffness matrix of one CST finite element is:

kel=Vel(σxGx+σyGy+τxyGxy+BTDB)\boldsymbol{k}_{el}=V_{el}\left(\sigma_x\boldsymbol{G}_x+\sigma_y\boldsymbol{G}_y+\tau_{xy}\boldsymbol{G}_{xy}+\boldsymbol{B}^{T}\boldsymbol{D}\boldsymbol{B}\right)

The first term is the geometric or stress-dependent part; the second is the material part. Equivalently:

Ψ(u)=Ael=1nelfel−F=0,fel=BTσVel\boldsymbol{\Psi}(\boldsymbol{u})=\mathcal{A}_{el=1}^{n_{el}}\boldsymbol{f}_{el}-\boldsymbol{F}=\boldsymbol{0},\qquad\boldsymbol{f}_{el}=\boldsymbol{B}^{T}\boldsymbol{\sigma}V_{el}
kel=∂fel∂uel,KT=Ael=1nelkel=∂Ψ∂u\boldsymbol{k}_{el}=\frac{\partial\boldsymbol{f}_{el}}{\partial\boldsymbol{u}_{el}},\qquad\boldsymbol{K}_T=\mathcal{A}_{el=1}^{n_{el}}\boldsymbol{k}_{el}=\frac{\partial\boldsymbol{\Psi}}{\partial\boldsymbol{u}}

08

Newton–Raphson solution

At Newton–Raphson iteration ii:

Δu(i)=−(KT(i))−1Ψ(i)\Delta\boldsymbol{u}^{(i)}=-\left(\boldsymbol{K}_T^{(i)}\right)^{-1}\boldsymbol{\Psi}^{(i)}
u(i+1)=u(i)+Δu(i)\boldsymbol{u}^{(i+1)}=\boldsymbol{u}^{(i)}+\Delta\boldsymbol{u}^{(i)}
Ψ(i)=Ψ(u(i))\boldsymbol{\Psi}^{(i)}=\boldsymbol{\Psi}(\boldsymbol{u}^{(i)})

The element internal-force vectors and tangent stiffness matrices are assembled in the standard way. A convenient stopping measure is:

ϵ=ΔuTΔun\epsilon=\sqrt{\frac{\Delta\boldsymbol{u}^{T}\Delta\boldsymbol{u}}{n}}

where nn is the total number of equations.

Important remark #1

The Total Lagrangian formulation presented here is not based on an incremental strain approximation. The equilibrium equations are written directly with respect to the initial configuration. Therefore, the accuracy of the formulation itself does not depend on the number of load steps. For problems with moderate nonlinearity, a single load step may be sufficient, provided that the Newton–Raphson iterations converge. Additional load steps may nevertheless be useful to improve convergence in more strongly nonlinear problems.

09

Example 1. Pure Bending of a Cantilever Beam

The cantilever has length L=400 mmL=400\ \mathrm{mm}, a rectangular cross-section with thickness 5 mm5\ \mathrm{mm} and height 20 mm20\ \mathrm{mm}, E=1000 MPaE=1000\ \mathrm{MPa}, and ν=0.3\nu=0.3. It is subjected to the end moment M=19000 N mmM=19000\ \mathrm{N\,mm}. The mesh has 60 subdivisions along the beam and 16 through its height: 1,920 CST elements and 1,037 nodes.

The analytical total rotation is:

θ=MLEI=19000⋅4001000⋅5⋅20312=2.267 rad=129.9∘\theta=\frac{ML}{EI}=\frac{19000\cdot400}{1000\cdot5\cdot\dfrac{20^3}{12}}=2.267\ \mathrm{rad}=129.9^\circ

This analytical value is practically identical to that obtained with the finite-element model.

Figure 1. Deformed configuration and normal stress σx\sigma_x under pure bending.
Figure 2. Green–Lagrange axial strain εx\varepsilon_x.

10

Example 2. Cantilever Beam Subjected to a Tip Force

The same beam geometry and material are used, with a vertical force F=100 NF=100\ \mathrm{N} applied at the free end. The mesh has 120 subdivisions along the beam and 16 through its height: 3,840 CST elements and 2,057 nodes.

ANSYS uses an Updated Lagrangian formulation for this large-displacement analysis and reports Cauchy stresses and logarithmic strains in the rotated element coordinate system. Since the strains in the present example remain small, these results can be meaningfully compared with the Green–Lagrange strains and second Piola–Kirchhoff stresses obtained with the present Total Lagrangian formulation.

The displacement results are:

Figure 3. Vertical displacement vv.
Figure 4. Horizontal displacement uu.
∣vmax⁡∣MATLAB=286.4 mm,∣umax⁡∣MATLAB=156.3 mm|v_{\max}|_{\mathrm{MATLAB}}=286.4\ \mathrm{mm},\qquad|u_{\max}|_{\mathrm{MATLAB}}=156.3\ \mathrm{mm}
∣vmax⁡∣ANSYS=287.1 mm,∣umax⁡∣ANSYS=156.8 mm|v_{\max}|_{\mathrm{ANSYS}}=287.1\ \mathrm{mm},\qquad|u_{\max}|_{\mathrm{ANSYS}}=156.8\ \mathrm{mm}

The normal-stress and Von Mises stress results are:

Figure 5. Second Piola–Kirchhoff normal stress σx\sigma_x.
Figure 6. Von Mises stress σVM\sigma_{\mathrm{VM}}.
σx,min⁡MATLAB=−78.9 MPa,σx,max⁡MATLAB=80.2 MPa\sigma_{x,\min}^{\mathrm{MATLAB}}=-78.9\ \mathrm{MPa},\qquad\sigma_{x,\max}^{\mathrm{MATLAB}}=80.2\ \mathrm{MPa}
σx,min⁡ANSYS=−80.4 MPa,σx,max⁡ANSYS=78.7 MPa\sigma_{x,\min}^{\mathrm{ANSYS}}=-80.4\ \mathrm{MPa},\qquad\sigma_{x,\max}^{\mathrm{ANSYS}}=78.7\ \mathrm{MPa}
σVM,max⁡MATLAB=74.9 MPa,σVM,max⁡ANSYS=74.1 MPa\sigma_{\mathrm{VM},\max}^{\mathrm{MATLAB}}=74.9\ \mathrm{MPa},\qquad\sigma_{\mathrm{VM},\max}^{\mathrm{ANSYS}}=74.1\ \mathrm{MPa}
Figure 7. Green–Lagrange axial strain εx\varepsilon_x.
εx,min⁡MATLAB=−0.0718,εx,max⁡MATLAB=0.0729\varepsilon_{x,\min}^{\mathrm{MATLAB}}=-0.0718,\qquad\varepsilon_{x,\max}^{\mathrm{MATLAB}}=0.0729
εx,min⁡ANSYS=−0.0732,εx,max⁡ANSYS=0.0716\varepsilon_{x,\min}^{\mathrm{ANSYS}}=-0.0732,\qquad\varepsilon_{x,\max}^{\mathrm{ANSYS}}=0.0716

The agreement between the MATLAB and ANSYS results is very good for displacements, stresses, and strains.

Important remark #2

In the Total Lagrangian formulation, strains and stresses are referred to the initial undeformed configuration. Therefore, the mathematically consistent stress map is the one represented on the initial geometry, as shown on the left.

For a more intuitive visualization, the same second Piola–Kirchhoff stresses may also be represented on the deformed geometry, as shown on the right. In this representation, the material directions rotate together with the finite elements; therefore, the x-stress direction becomes approximately tangent to the deformed beam axis.

Cauchy stresses are defined in the current configuration. When strains are small, their values referred to the corresponding local, rotated material directions are very close to the second Piola–Kirchhoff stresses. The distinction becomes important when comparing stress components expressed in different reference frames.

Normal stress represented on the initial and deformed beam geometries, with arrows showing the rotated material direction
Figure 8. The stress direction is fixed with respect to the initial reference frame on the left and rotates with the material directions in the deformed visualization on the right. The explanatory arrows are retained from the author's figure.

11

Programs and downloads

MMATLAB packageFive program files and README · tested in MATLAB R2026aDownload ZIPAANSYS comparison packageMechanical APDL macro and README · inspected, not executedDownload ZIP

The MATLAB package was executed in MATLAB R2026a. All five load steps reached convergence, and the numerical values reproduced the documented results. The ANSYS macro was inspected for consistency with the geometry, thickness, material, plane-stress assumption, loading, constraints, and mesh of Example 2, but it was not executed.