Isoparametric 8-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_8$ and $y_1,\ldots,y_8$ are the nodal coordinates of the finite element in the initial configuration. The Q8 shape functions $N_1,\ldots,N_8$ and node numbering are those of Section 9.6. The physical area transformation is:

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

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 $16\times16$ 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.

For Q8, $\htmlClass{fea-math-bold-italic}{\boldsymbol B}$, $\htmlClass{fea-math-bold-italic}{\boldsymbol B}_0$ and $\htmlClass{fea-math-bold-italic}{\boldsymbol B}_L$ have dimensions $3\times16$.

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 $3\times3$ Gauss numerical quadrature:

The integration uses nine Gauss points, three in each natural-coordinate direction. Their coordinates and one-dimensional weights are:

The two-dimensional weight is the product of the weights in the two natural-coordinate directions. The nine-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 $16\times16$ 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 60 quadrilateral elements along the beam and 8 elements through its height: 480 isoparametric 8-node quadrilateral finite elements, 1,577 nodes and 3,154 displacement degrees of freedom. The load is applied in five steps. The Newton–Raphson solution requires 8, 9, 8, 8 and 7 iterations, respectively, for a tol=1e-5 tolerance.

The comparison with Section 12.2 uses the same geometry, material and applied force. The Q4 and Q8 displacement solutions were calculated with their respective MATLAB programs; the Q4 strain/stress values and the ANSYS values below are those reported in Section 12.2. Displacements are compared by their maximum absolute values:

The Q8 results are close to the Q4 and ANSYS results, with a smaller finite-element discretization.