12.2 Isoparametric 4-Node Quadrilateral Finite Element

Formulation and hypotheses

StrainSmall

DisplacementLarge

MaterialLinear elastic

Stress statePlane stress

In the Total Lagrangian formulation, the initial undeformed geometry is the reference configuration. All coordinates, derivatives, Jacobian terms, and integrals are referred to that initial configuration. Green–Lagrange strains are used.

The deformation energy is:

U=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\int_{A_{el}}h\,\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}\,dA\right)

$\mathcal A_{el=1}^{n_{el}}$ denotes the standard finite-element assembly process, from the first element to the last one. In the natural coordinates $s,t$ :

U=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\int_{-1}^{1}\int_{-1}^{1}h(s,t)\,\boldsymbol{\varepsilon}^{T}(s,t)\boldsymbol{\sigma}(s,t)\,\left|\det\boldsymbol{J}(s,t)\right|\,ds\,dt\right)

Here, $h$ is the element thickness, $A_{el}$ is the finite-element area, and $n_{el}$ is the number of finite elements. Engineering notation is used for the shear component:

\gamma_{xy}=2E_{xy}

\boldsymbol{\varepsilon}=\begin{Bmatrix}\varepsilon_x\\\varepsilon_y\\\gamma_{xy}\end{Bmatrix}

\boldsymbol{\sigma}=\begin{Bmatrix}\sigma_x\\\sigma_y\\\tau_{xy}\end{Bmatrix}

For a linear elastic material in plane stress:

\boldsymbol{\sigma}=\boldsymbol{D}\boldsymbol{\varepsilon}

\boldsymbol{D}=\frac{E}{1-\nu^2}\begin{bmatrix}1&\nu&0\\\nu&1&0\\0&0&\dfrac{1-\nu}{2}\end{bmatrix}

Initial configuration and Jacobian

The Jacobian is:

\boldsymbol{J}(s,t)=\begin{bmatrix}\dfrac{\partial x}{\partial s}&\dfrac{\partial y}{\partial s}\\[4pt]\dfrac{\partial x}{\partial t}&\dfrac{\partial y}{\partial t}\end{bmatrix}=\boldsymbol{H}(s,t)\begin{bmatrix}x_1&y_1\\x_2&y_2\\x_3&y_3\\x_4&y_4\end{bmatrix}

\boldsymbol{H}(s,t)=\frac14\begin{bmatrix}1+t&-(1+t)&-(1-t)&1-t\\1+s&1-s&-(1-s)&-(1+s)\end{bmatrix}

$x_1,\ldots,x_4$ and $y_1,\ldots,y_4$ are the nodal coordinates of the finite element in the initial configuration. The physical area transformation is:

dA=\left|\det\boldsymbol{J}\right|\,ds\,dt

The derivatives of the shape functions with respect to the initial Cartesian coordinates are obtained using the inverse Jacobian; see Section 9.4:

\boldsymbol{b}=\boldsymbol{J}^{-1}\boldsymbol{H}=\begin{bmatrix}b_{11}&b_{12}&b_{13}&b_{14}\\b_{21}&b_{22}&b_{23}&b_{24}\end{bmatrix}

The four displacement-derivative relations are:

\begin{aligned}\frac{\partial u}{\partial x}&=\boldsymbol{B}_{u_x}\boldsymbol{u}_{el}\\[4pt]\boldsymbol{B}_{u_x}&=\begin{pmatrix}b_{11}&0&b_{12}&0&b_{13}&0&b_{14}&0\end{pmatrix}\end{aligned}

\begin{aligned}\frac{\partial u}{\partial y}&=\boldsymbol{B}_{u_y}\boldsymbol{u}_{el}\\[4pt]\boldsymbol{B}_{u_y}&=\begin{pmatrix}b_{21}&0&b_{22}&0&b_{23}&0&b_{24}&0\end{pmatrix}\end{aligned}

\begin{aligned}\frac{\partial v}{\partial x}&=\boldsymbol{B}_{v_x}\boldsymbol{u}_{el}\\[4pt]\boldsymbol{B}_{v_x}&=\begin{pmatrix}0&b_{11}&0&b_{12}&0&b_{13}&0&b_{14}\end{pmatrix}\end{aligned}

\begin{aligned}\frac{\partial v}{\partial y}&=\boldsymbol{B}_{v_y}\boldsymbol{u}_{el}\\[4pt]\boldsymbol{B}_{v_y}&=\begin{pmatrix}0&b_{21}&0&b_{22}&0&b_{23}&0&b_{24}\end{pmatrix}\end{aligned}

The element nodal displacement vector is:

\boldsymbol{u}_{el}=\begin{Bmatrix}u_1\\v_1\\u_2\\v_2\\u_3\\v_3\\u_4\\v_4\end{Bmatrix}

Green–Lagrange strains

The two-dimensional Green–Lagrange strain components are:

\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

\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

\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}

After introducing the displacement interpolation, the strains become:

\varepsilon_x=\boldsymbol{B}_{u_x}\boldsymbol{u}_{el}+\frac12\boldsymbol{u}_{el}^{T}\boldsymbol{G}_x\boldsymbol{u}_{el}

\varepsilon_y=\boldsymbol{B}_{v_y}\boldsymbol{u}_{el}+\frac12\boldsymbol{u}_{el}^{T}\boldsymbol{G}_y\boldsymbol{u}_{el}

\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 symmetric $8\times8$ matrices are:

\boldsymbol{G}_x=\boldsymbol{B}_{u_x}^{T}\boldsymbol{B}_{u_x}+\boldsymbol{B}_{v_x}^{T}\boldsymbol{B}_{v_x}

\boldsymbol{G}_y=\boldsymbol{B}_{u_y}^{T}\boldsymbol{B}_{u_y}+\boldsymbol{B}_{v_y}^{T}\boldsymbol{B}_{v_y}

\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}

Strain-displacement matrix

The matrix $\boldsymbol B$ , relating an increment of the Green–Lagrange strain vector to an increment of the element nodal displacement vector, is defined by:

d\boldsymbol{\varepsilon}_{el}=\boldsymbol{B}\,d\boldsymbol{u}_{el}

\boldsymbol B=\frac{\partial\boldsymbol\varepsilon_{el}}{\partial\boldsymbol u_{el}}

Component by component:

d\varepsilon_x=\left(\boldsymbol B_{u_x}+\boldsymbol u_{el}^{T}\boldsymbol G_x\right)d\boldsymbol u_{el}

d\varepsilon_y=\left(\boldsymbol B_{v_y}+\boldsymbol u_{el}^{T}\boldsymbol G_y\right)d\boldsymbol u_{el}

d\gamma_{xy}=\left(\boldsymbol B_{u_y}+\boldsymbol B_{v_x}+\boldsymbol u_{el}^{T}\boldsymbol G_{xy}\right)d\boldsymbol u_{el}

\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 part, identical to linear theory, and a nonlinear displacement-dependent part:

\boldsymbol{B}(\boldsymbol{u}_{el})=\boldsymbol{B}_0+\boldsymbol{B}_L(\boldsymbol{u}_{el})

\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}

$\boldsymbol B_0$ is the constant part, identical to linear theory, and $\boldsymbol B_L$ is the nonlinear displacement-dependent part.

Potential energy and Gauss quadrature

The total potential energy is:

\Pi=U+V=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\int_{A_{el}}h\,\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}\,dA\right)-\boldsymbol{u}^{T}\boldsymbol{F}

Here, $V=-\boldsymbol{u}^{T}\boldsymbol F$ is the potential of the external loads, $\boldsymbol u$ is the global nodal displacement vector, and $\boldsymbol F$ is the global force vector. In natural coordinates:

\Pi=U+V=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\int_{-1}^{1}\int_{-1}^{1}h(s,t)\,\boldsymbol{\varepsilon}^{T}(s,t)\boldsymbol{\sigma}(s,t)\,\left|\det\boldsymbol J(s,t)\right|\,ds\,dt\right)-\boldsymbol{u}^{T}\boldsymbol F

Using $2\times2$ Gauss numerical quadrature:

\Pi=U+V=\frac12\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h(s_G,t_G)\,\boldsymbol{\varepsilon}^{T}(s_G,t_G)\boldsymbol{\sigma}(s_G,t_G)\,\left|\det\boldsymbol J(s_G,t_G)\right|\right)-\boldsymbol{u}^{T}\boldsymbol F

The integration uses four Gauss points, two in each natural-coordinate direction. For the $2\times2$ Gauss rule, all four weights are equal to 1 and are therefore omitted from the following expressions. The four-point sum is:

\sum_{n_G}E(s_G,t_G)=E\left(\frac{1}{\sqrt3},\frac{1}{\sqrt3}\right)+E\left(\frac{1}{\sqrt3},-\frac{1}{\sqrt3}\right)+E\left(-\frac{1}{\sqrt3},\frac{1}{\sqrt3}\right)+E\left(-\frac{1}{\sqrt3},-\frac{1}{\sqrt3}\right)

where the integrand is:

E(s,t)=h(s,t)\,\boldsymbol{\varepsilon}^{T}(s,t)\boldsymbol{\sigma}(s,t)\,\left|\det\boldsymbol J(s,t)\right|

Virtual work and equilibrium

The virtual-work equation is:

\delta\Pi=\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\,\delta\boldsymbol{\varepsilon}^{T}\boldsymbol{\sigma}\,\left|\det\boldsymbol J\right|\right)-\delta\boldsymbol{u}^{T}\boldsymbol F=0

This relation holds for every geometrically admissible virtual displacement field. Therefore, the nonlinear equilibrium system is:

\boldsymbol{\Psi}(\boldsymbol u)=\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\,\boldsymbol B^{T}\boldsymbol\sigma\,\left|\det\boldsymbol J\right|\right)-\boldsymbol F=\boldsymbol0

Equivalently, using the split of the strain-displacement matrix:

\boldsymbol{\Psi}(\boldsymbol u)=\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\left[\boldsymbol B_0^{T}+\boldsymbol B_L^{T}(\boldsymbol u_{el})\right]\boldsymbol\sigma\,\left|\det\boldsymbol J\right|\right)-\boldsymbol F=\boldsymbol0

Tangent stiffness matrix

The first variation of the equilibrium residual is:

\delta\boldsymbol\Psi=\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\,\delta\boldsymbol B^{T}\boldsymbol\sigma\,\left|\det\boldsymbol J\right|\right)+\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\,\boldsymbol B^{T}\delta\boldsymbol\sigma\,\left|\det\boldsymbol J\right|\right)=\boldsymbol{K}_T\delta\boldsymbol u

Using $\boldsymbol B=\boldsymbol B_0+\boldsymbol B_L(\boldsymbol u_{el})$ , the same variation is:

\delta\boldsymbol\Psi=\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\,\delta\left[\boldsymbol B_0^{T}+\boldsymbol B_L^{T}(\boldsymbol u_{el})\right]\boldsymbol\sigma\,\left|\det\boldsymbol J\right|\right)+\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\left[\boldsymbol B_0^{T}+\boldsymbol B_L^{T}(\boldsymbol u_{el})\right]\delta\boldsymbol\sigma\,\left|\det\boldsymbol J\right|\right)=\boldsymbol{K}_T\delta\boldsymbol u

Because $\boldsymbol B_0$ is constant:

\delta\boldsymbol B=\delta(\boldsymbol B_0+\boldsymbol B_L)=\delta\boldsymbol B_L

Taking into account that $\boldsymbol G_x$ , $\boldsymbol G_y$ , and $\boldsymbol G_{xy}$ are symmetric $8\times8$ matrices, the first contribution becomes:

\delta\boldsymbol B^{T}\boldsymbol\sigma=\delta\boldsymbol B_L^{T}\boldsymbol\sigma=\delta\left(\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}\right)\boldsymbol\sigma

\delta\boldsymbol B^{T}\boldsymbol\sigma=\left(\sigma_x\boldsymbol G_x+\sigma_y\boldsymbol G_y+\tau_{xy}\boldsymbol G_{xy}\right)\delta\boldsymbol u_{el}

The stress variation is:

\delta\boldsymbol\sigma=\boldsymbol D\,\delta\boldsymbol\varepsilon=\boldsymbol D\boldsymbol B\,\delta\boldsymbol u_{el}

Consequently, the global tangent stiffness matrix is:

\boldsymbol{K}_T=\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\left[\sigma_x\boldsymbol G_x+\sigma_y\boldsymbol G_y+\tau_{xy}\boldsymbol G_{xy}\right]\left|\det\boldsymbol J\right|\right)+\mathcal{A}_{el=1}^{n_{el}}\left(\sum_{n_G}h\,\boldsymbol B^{T}\boldsymbol D\boldsymbol B\,\left|\det\boldsymbol J\right|\right)

The tangent stiffness matrix contains two contributions: the material stiffness and the geometric (stress-dependent) stiffness. For one finite element:

\boldsymbol k_{el}=\sum_{n_G}h\left(\sigma_x\boldsymbol G_x+\sigma_y\boldsymbol G_y+\tau_{xy}\boldsymbol G_{xy}+\boldsymbol B^{T}\boldsymbol D\boldsymbol B\right)\left|\det\boldsymbol J\right|

The equilibrium system can also be written in terms of the element internal force vectors:

\boldsymbol\Psi(\boldsymbol u)=\mathcal{A}_{el=1}^{n_{el}}\left(\boldsymbol f_{el}\right)-\boldsymbol F=\boldsymbol0

\boldsymbol f_{el}=\sum_{n_G}h\,\boldsymbol B^{T}\boldsymbol\sigma\,\left|\det\boldsymbol J\right|

The element tangent stiffness matrix is the derivative of the element internal force vector:

\boldsymbol k_{el}=\frac{\partial\boldsymbol f_{el}}{\partial\boldsymbol u_{el}}

Newton–Raphson solution

At Newton–Raphson iteration $i$ :

\Delta\boldsymbol u^{(i)}=-\left(\boldsymbol K_T^{(i)}\right)^{-1}\boldsymbol\Psi^{(i)}

\boldsymbol u^{(i+1)}=\boldsymbol u^{(i)}+\Delta\boldsymbol u^{(i)}

\boldsymbol\Psi^{(i)}=\boldsymbol\Psi\left(\boldsymbol u_i\right)

The element internal-force vectors $\boldsymbol f_{el}$ and tangent stiffness matrices $\boldsymbol k_{el}$ are assembled to obtain $\boldsymbol\Psi^{(i)}$ and $\boldsymbol K_T^{(i)}$ , respectively. A convenient convergence measure is:

\epsilon=\sqrt{\frac{\Delta\boldsymbol u^{T}\Delta\boldsymbol u}{n}}

where $n$ is the 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.

Example and numerical results

The cantilever beam has length $L=400\ \mathrm{mm}$ and a rectangular cross-section with thickness $5\ \mathrm{mm}$ and height $20\ \mathrm{mm}$ . The material properties are $E=1000\ \mathrm{MPa}$ and $\nu=0.3$ . A vertical force $F=100\ \mathrm{N}$ is applied at the free end.

The mesh has 120 quadrilateral elements along the beam and 16 elements through its height: 1,920 isoparametric 4-node quadrilateral finite elements and 2,057 nodes.

Horizontal displacement contour map of the deformed cantilever — **Figure 1.** Horizontal displacement $u$ .

u_{\max}^{\mathrm{MATLAB}}=159.51\ \mathrm{mm}

u_{\max}^{\mathrm{ANSYS}}=159.92\ \mathrm{mm}

Vertical displacement contour map of the deformed cantilever — **Figure 2.** Vertical displacement $v$ .

v_{\max}^{\mathrm{MATLAB}}=288.79\ \mathrm{mm}

v_{\max}^{\mathrm{ANSYS}}=289.48\ \mathrm{mm}

Green–Lagrange axial strain contour map — **Figure 3.** Green–Lagrange axial strain $\varepsilon_x$ .

\varepsilon_{x,\min}^{\mathrm{MATLAB}}=-0.0726

\varepsilon_{x,\min}^{\mathrm{ANSYS}}=-0.0770

\varepsilon_{x,\max}^{\mathrm{MATLAB}}=0.0738

\varepsilon_{x,\max}^{\mathrm{ANSYS}}=0.0750

Second Piola–Kirchhoff normal stress contour map — **Figure 4.** Second Piola–Kirchhoff normal stress $\sigma_x$ .

\sigma_{x,\min}^{\mathrm{MATLAB}}=-82.72\ \mathrm{MPa}

\sigma_{x,\min}^{\mathrm{ANSYS}}=-84.53\ \mathrm{MPa}

\sigma_{x,\max}^{\mathrm{MATLAB}}=84.62\ \mathrm{MPa}

\sigma_{x,\max}^{\mathrm{ANSYS}}=82.63\ \mathrm{MPa}

Von Mises stress contour map — **Figure 5.** Von Mises stress $\sigma_{\mathrm{VM}}$ .

\sigma_{\mathrm{VM},\min}^{\mathrm{MATLAB}}=0.189\ \mathrm{MPa}

\sigma_{\mathrm{VM},\min}^{\mathrm{ANSYS}}=0.152\ \mathrm{MPa}

\sigma_{\mathrm{VM},\max}^{\mathrm{MATLAB}}=79.54\ \mathrm{MPa}

\sigma_{\mathrm{VM},\max}^{\mathrm{ANSYS}}=79.35\ \mathrm{MPa}

The MATLAB and ANSYS results are in very good agreement.

Stress-map interpretation

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.

For a more intuitive visualization, the same second Piola–Kirchhoff stresses may also be represented on the deformed geometry. 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.

This topic will be discussed in a later chapter.

Normal stress represented on the initial and deformed beam geometries with the author-supplied arrows and labels — **Figure 6.** Second Piola–Kirchhoff normal stress represented on the initial configuration and, for intuitive visualization, on the deformed geometry. The original explanatory arrows and labels are preserved.

Programs and downloads

MMATLAB packageSix program files and README · tested in MATLAB R2026aDownload ZIP AANSYS comparison packageMechanical APDL macro and README · historically validated, not re-executedDownload ZIP