Isoparametric 4-Node Quadrilateral Finite Element

01

Formulation and hypotheses

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:

$\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$:

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:

For a linear elastic material in plane stress:

02

Initial configuration and Jacobian

The Jacobian is:

$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:

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

The four displacement-derivative relations are:

The element nodal displacement vector is:

03

Green–Lagrange strains

The two-dimensional Green–Lagrange strain components are:

After introducing the displacement interpolation, the strains become:

The symmetric $8\times8$ matrices are:

04

Strain-displacement matrix

The matrix $\htmlClass{fea-math-bold-italic}{\boldsymbol B}$, relating an increment of the Green–Lagrange strain vector to an increment of the element nodal displacement vector, is defined by:

Component by component:

The matrix is split into a constant part, identical to linear theory, and a nonlinear displacement-dependent part:

$\htmlClass{fea-math-bold-italic}{\boldsymbol B}_0$ is the constant part, identical to linear theory, and $\htmlClass{fea-math-bold-italic}{\boldsymbol B}_L$ is the nonlinear displacement-dependent part.

05

Potential energy and Gauss quadrature

The total potential energy is:

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

Using $2\times2$ Gauss numerical quadrature:

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:

where the integrand is:

06

Virtual work and equilibrium

The virtual-work equation is:

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

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

07

Tangent stiffness matrix

The first variation of the equilibrium residual is:

Using $\htmlClass{fea-math-bold-italic}{\boldsymbol B}=\htmlClass{fea-math-bold-italic}{\boldsymbol B}_0+\htmlClass{fea-math-bold-italic}{\boldsymbol B}_L(\htmlClass{fea-math-bold-italic}{\boldsymbol u}_{el})$, the same variation is:

Because $\htmlClass{fea-math-bold-italic}{\boldsymbol B}_0$ is constant:

Taking into account that $\htmlClass{fea-math-bold-italic}{\boldsymbol G}_x$, $\htmlClass{fea-math-bold-italic}{\boldsymbol G}_y$, and $\htmlClass{fea-math-bold-italic}{\boldsymbol G}_{xy}$ are symmetric $8\times8$ matrices, the first contribution becomes:

The stress variation is:

Consequently, the global tangent stiffness matrix is:

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

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

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

08

Newton–Raphson solution

At Newton–Raphson iteration $i$:

The element internal-force vectors $\htmlClass{fea-math-bold-italic}{\boldsymbol f}_{el}$ and tangent stiffness matrices $\htmlClass{fea-math-bold-italic}{\boldsymbol k}_{el}$ are assembled to obtain $\htmlClass{fea-math-bold-italic}{\boldsymbol{\mathit{\Psi}}}^{(i)}$ and $\htmlClass{fea-math-bold-italic}{\boldsymbol K}_T^{(i)}$, respectively. A convenient convergence measure is:

where $n$ is the number of equations.

09

Example and numerical results

The cantilever beam has length $L=400\ \mathit{mm}$ and a rectangular cross-section with thickness $5\ \mathit{mm}$ and height $20\ \mathit{mm}$. The material properties are $E=1000\ \mathit{MPa}$ and $\nu=0.3$. A vertical force $F=100\ \mathit{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.

$$u_{\mathop{\mathit{max}}}^{\mathit{MATLAB}}=159.51\ \mathit{mm}$$

$$u_{\mathop{\mathit{max}}}^{\mathit{ANSYS}}=159.92\ \mathit{mm}$$

$$v_{\mathop{\mathit{max}}}^{\mathit{MATLAB}}=288.79\ \mathit{mm}$$

$$v_{\mathop{\mathit{max}}}^{\mathit{ANSYS}}=289.48\ \mathit{mm}$$

The MATLAB and ANSYS results are in very good agreement.

10

Stress-map interpretation

11

Programs and downloads