Plane Stress: Isoparametric 4-Node Quadrilateral Finite Element — A Slightly Different Approach

01

Reference configuration and hypotheses

The isoparametric 4-node quadrilateral finite element is briefly described in Section 9.4.

Updated Lagrangian formulation. The equations for load step $i+1$ are written with respect to the known converged configuration at load step $i$, which becomes the reference configuration for the new increment. Displacements and stresses at the target configuration are therefore determined relative to that preceding configuration, following the reference-configuration discussion of Chatzi [1].

In Section 15.1, the geometry is updated after every iteration. In the present section, the converged configuration at load step $i$ is kept fixed during all iterations used to determine load step $i+1$. Only after convergence is the geometry updated to the new configuration.

02

Strain increment

The displacement and strain increments are measured relative to configuration $i$. All derivatives in the following relations are evaluated with respect to that same configuration. The engineering components of the strain increment are:

The strain increment has the Green–Lagrange form, but it is evaluated with respect to the configuration at the previous load step, not the initial undeformed configuration. Since two successive configurations are close when the load increment is small, the corresponding Green–Lagrange and Euler–Almansi strain increments are also very close.

03

Geometry-update strategy

The essential programming difference between Sections 15.1 and 15.2 is shown below:

while (error > tol) & (iter < itermax)
    iter=iter+1;
    stiff
    dS=K\F;
    S=S-dS;
    error=sqrt(dS'*dS/neq);
    x=x-dS(1:2:neq);
    y=y-dS(2:2:neq);
end

while (error > tol) & (iter < itermax)
    iter=iter+1;
    stiff
    dS=K\F;
    S=S-dS;
    error=sqrt(dS'*dS/neq);
end
x=x+S(1:2:neq);
y=y+S(2:2:neq);

This apparently small programming change has an important consequence: during all iterations of one load step, the geometry remains fixed at the converged configuration of the preceding load step. The accumulated displacement increment $S$ is determined with respect to this fixed configuration. Only after convergence are the nodal coordinates updated.

In the MATLAB program, dS is the correction calculated during one iteration, S is the displacement increment accumulated during the current load step, and St(:,istep+1) is the total displacement accumulated up to that load step.

Apart from the different geometry-update strategy and the strain-update procedure discussed below, the finite-element formulation introduced in Section 15.1 remains applicable here.

04

Simple strain update

The routine stiff uses the simple strain update:

A simple additive strain update would be exact only if the strain measure depended linearly on the deformation. Green–Lagrange and Euler–Almansi strains are nonlinear measures, so such an update is only approximate. For sufficiently small strain increments the error is small and decreases as the load steps are refined.

In the additive expression, $\boldsymbol e_i$ and $\Delta\boldsymbol E$ are associated with configuration $i$, whereas the updated strain is required in configuration $i+1$. The transformation introduced next removes this strain-update kinematic inconsistency.

05

Transformed strain update

Let $\Delta\boldsymbol E$ denote the Green–Lagrange strain increment for the deformation from configuration $i$ to configuration $i+1$, with configuration $i$ used as reference. Let $\boldsymbol e_i$ denote the Euler–Almansi strain already accumulated and expressed in configuration $i$.

The deformation gradient for the increment is:

Here, $\boldsymbol F$ is the deformation gradient for the increment from configuration $i$ to configuration $i+1$, $\boldsymbol I$ is the identity matrix, and $\nabla_i\Delta\boldsymbol u$ is the gradient of the displacement increment with respect to the coordinates of configuration $i$.

In component form:

The Green–Lagrange strain increment is:

The corresponding Euler–Almansi strain increment in configuration $i+1$ is:

The exact transformation between these two representations of the increment is:

The previously accumulated strain must also be transformed to the new current configuration. Therefore:

Both difficulties are removed by transforming the accumulated strain and the new strain increment to the current configuration.

Thus, stiff2 removes the kinematic error associated with the simple additive strain update by transforming the accumulated strain to the new current configuration.

06

MATLAB implementation

The additional transformation in stiff2 is implemented as follows:

FF=[1+ux uy; vx 1+vy];
FF=inv(FF);
ee=[ex exy/2; exy/2 ey];
ee=FF'*ee*FF;
ex =ee(1,1);
ey =ee(2,2);
exy=2*ee(1,2);

After FF=inv(FF);, the MATLAB variable FF contains $\boldsymbol F^{-1}$. Consequently, ee=FF'*ee*FF; performs:

Engineering shear is divided by 2 when the symmetric tensor ee is constructed and multiplied by 2 when exy is recovered.

Both stiff and stiff2 preserve the nonlinear strain–displacement matrix:

The nonlinear matrices retained in the tangent are:

The present implementation retains the nonlinear $\boldsymbol G$-terms in the tangent stiffness. These terms are preserved as part of the numerical procedure used here.

The strain update used in stiff2 is corrected by the tensor transformation described above. The internal-force vector is then evaluated from the resulting current strain and stress state. The tangent stiffness is used to obtain the iterative displacement correction and need not be a fully consistent tangent for the method to converge to equilibrium.

Using an approximate tangent may affect the convergence rate, but not the equilibrium condition itself, provided that the internal-force vector is evaluated consistently with the adopted formulation.

07

Plane-stress thickness update

The small-strain plane-stress thickness convention is:

The current thickness is approximated by:

This thickness relation is a small-strain approximation; it is not a general exact relation for finite elastic strains.

08

Example 1: cantilever comparison

A cantilever of length $L=400\ \mathrm{mm}$ has a rectangular cross-section $5\ \mathrm{mm}\times20\ \mathrm{mm}$, Young's modulus $E=1000\ \mathrm{MPa}$, and Poisson ratio $\nu=0.3$. A vertical force $F=100\ \mathrm{N}$ is applied at the free end. The mesh contains $120\times16=1920$ Q4 elements and 2057 nodes.

The comparison uses five load steps for both Updated Lagrangian variants. The distributed program is configured accordingly with nstep=5.

Green — TLBlue — UL with stiffRed — UL with stiff25 load steps

The table retains the notation $v_{\max}$ used for the maximum downward displacement; its signed value is negative:

In the TL formulation, the number of load steps affects the nonlinear solution procedure but does not introduce an incremental strain-accumulation error.

For the displacement results shown here, the UL formulation using stiff2 with only five load steps is already very close to the TL solution and to the ANSYS result.

09

Programs and download

A MATLAB program package for this example can be downloaded below.

10

Reference

1. Eleni Chatzi, The Finite Element Method for the Analysis of Non-Linear and Dynamic Systems, ETH Zürich, Lecture 3, 15 October 2015.