SECTION 15.5 OF CHAPTER 15: UPDATED LAGRANGIAN FORMULATION (UL)
Updated Lagrangian Revisited: A Direct and Consistent Tangent Stiffness Matrix for the Q4 Element
This section was developed with the contribution of Chiara.
01
Introduction
StrainSmall
DisplacementLarge
MaterialLinear elastic
Stress statePlane stress
In Section 15.1, the Updated Lagrangian (UL) formulation was developed using a simplified strain–displacement matrix. In particular, the variation of the element geometry through the Jacobian determinant was neglected in the tangent stiffness matrix, and the variation of the Euler–Almansi strains was treated approximately.
Sections 15.2–15.4 introduced alternative incremental formulations combining Euler–Almansi and Green–Lagrange strain measures. The nodal coordinates are updated after each load step, and the subsequent calculations are performed with reference to the updated configuration. These approaches provide accurate results even for relatively few load steps; their accuracy generally improves as the number of load steps increases.
In this section, a direct formulation with a consistent tangent stiffness matrix is developed for the four-node isoparametric Q4 element. Total Euler–Almansi strains are calculated directly from the initial and current configurations, without accumulating strain increments. The tangent stiffness matrix is obtained by consistently linearizing the internal-force vector and accounting for every factor that depends on the nodal displacements.
The converged results for the numerical example are practically independent of the number of load steps, within the prescribed tolerance. The solutions obtained with the approximate formulations of the preceding sections approach these results as the load-step size decreases.
The derivation follows the spatial description of nonlinear equilibrium and its linearization discussed by Bonet and Wood [1], here specialized to the Q4 element and to the adopted plane-stress constitutive model.
02
Element geometry and strain–displacement matrix
The equilibrium equations are written with respect to the current configuration, corresponding to the current load step and Newton iteration. The nodal coordinates are updated after each displacement correction. The interpolation functions remain unchanged, but their spatial derivatives are evaluated using the deformed geometry.
As in Section 9.4, the Jacobian matrix of the Q4 element is:
Although B has the same form as in linear elasticity, its coefficients depend on the current nodal coordinates and therefore change during the Newton iterations.
03
Internal nodal forces and constitutive equations
The internal nodal-force vector is obtained from the Cauchy stresses:
fel=∫Ath1BTσvdA,
where At denotes the current element area and h1 is the current thickness. In natural coordinates,
fel=∫−11∫−11h1BT(s,t)σv(s,t)detJ(s,t)dsdt.
Throughout this section, σ denotes the Cauchy stress tensor, whereas σv denotes its column-vector representation in Voigt notation. The subscript v distinguishes the vector representation from the tensor itself. Likewise, εv is the Voigt strain vector; its shear component is the engineering shear strain γxy=2exy:
For the adopted plane-stress constitutive approximation,
σv=Dεv,D=1−ν2E1ν0ν1000(1−ν)/2.
This relation is applied between Cauchy stresses and total Euler–Almansi strains. It is the constitutive approximation adopted here, not a general finite-strain hyperelastic law.
The total Euler–Almansi strains, introduced in Chapter 13, may be written in terms of the spatial derivatives of the total displacements:
εx=∂x∂u−21(∂x∂u)2−21(∂x∂v)2,
εy=∂y∂v−21(∂y∂u)2−21(∂y∂v)2,
γxy=∂y∂u+∂x∂v−∂x∂u∂y∂u−∂x∂v∂y∂v.
For the following derivation, the tensor representation is more convenient.
04
Total Euler–Almansi strains
The Euler–Almansi strain tensor is:
e=21(I2−Ftot−TFtot−1).
Here Ftot is the total deformation gradient from the initial to the current configuration:
The coordinates x0,y0 refer to the initial configuration, x,y to the current configuration, and u,v to the total displacements.
For an isoparametric element, Ftot can be evaluated directly from the initial and current Jacobian matrices:
Ftot=JTJ0−T=(J0−1J)T,
where:
J0=∂(s,t)∂(x0,y0),J=∂(s,t)∂(x,y).
Thus, the total Euler–Almansi strains are obtained directly from the initial and current nodal coordinates, without accumulating strain increments over the load steps.
05
Variation of the Euler–Almansi strains
The aim is to express the infinitesimal variation of the Euler–Almansi strains in terms of the virtual nodal displacements. Starting from
e=21(I2−Ftot−TFtot−1),
we obtain:
δe=−21[δFtot−TFtot−1+Ftot−TδFtot−1].
The variation of the total deformation gradient is:
where the infinitesimal spatial gradient of the virtual displacement is:
δL=∇δu=[∂x∂(δu)∂x∂(δv)∂y∂(δu)∂y∂(δv)].
All spatial derivatives are evaluated in the current configuration. We use δL to emphasize the infinitesimal nature of the virtual displacement gradient.
The variation of the inverse of a nonsingular matrix satisfies
δFtot−1=−Ftot−1δFtotFtot−1,
δFtot−T=−Ftot−TδFtotTFtot−T.
Using δFtot=δLFtot, these become
δFtot−1=−Ftot−1δL,δFtot−T=−δLTFtot−T.
Introducing the symmetric matrix
Q=(FtotFtotT)−1=Ftot−TFtot−1,
the variation of the Euler–Almansi strains becomes
δe=21(QδL+δLTQ).
This expression is linear in the virtual displacement gradient, although Q depends on the current deformation. For small displacements and rotations, Ftot≈I2 and Q≈I2, giving the familiar infinitesimal strain variation:
δe≈21(δL+δLT).
For the Q4 element, the virtual displacement field is interpolated as
δu=j=1∑4Njδuj,δuj=[δujδvj].
Consequently:
δL=j=1∑4δujbjT=[δu1δu2δu3δu4]bT.
Substituting in the strain variation:
δe=21j=1∑4[QδujbjT+bjδujTQ].
In Voigt notation, recalling that γxy=2exy, the previous relation gives
Thus the blocks of the standard strain–displacement matrix are multiplied by Q. This operation is performed directly at each Gauss point in MATLAB:
BQ=zeros(3,8);
for i=1:4
ii=2*i-1:2*i;
BQ(:,ii)=B(:,ii)*Q;
end
For Q=I2, the standard matrix is recovered: BQ=B.
06
Consistent linearization of the internal forces
The element tangent stiffness matrix is defined as
kel=∂uel∂fel.
The variation is infinitesimal and is evaluated about the current configuration. Using 2×2 Gauss integration, the element internal-force vector is approximated by
fel=p=1∑4hpBpTσv,pdetJp.
All integration weights are equal to one for this rule. Four factors depend on the current nodal displacements: the thickness hp, the matrix Bp, the Cauchy stresses σv,p, and the Jacobian determinant detJp.
The first variation of the internal-force vector is therefore
6.1. Variation of the Euler–Almansi strains: material tangent
From the preceding derivation,
δεv=BQδuel.
As D is constant for the adopted constitutive model,
δσv=DBQδuel.
The first contribution to the tangent stiffness matrix is therefore
kelm=p=1∑4hpBpTDBQ,pdetJp.
This is the material contribution. The two strain–displacement matrices play different roles: BpT comes from the spatial internal-force expression, whereas BQ,p describes the variation of the total Euler–Almansi strains.
6.2. Variation of the strain–displacement matrix
The matrix B depends on the spatial gradients of the interpolation functions. These gradients change with the current element geometry.
For node i,
BiTσv=[σxτxyτxyσy]bi=σbi.
Since bi=J−1Hi, where Hi is the ith column of H and is independent of the nodal displacements,
δbi=δ(J−1)Hi.
The variation of an inverse matrix satisfies
δ(J−1)=−J−1(δJ)J−1.
The variation of the current Jacobian is:
δJ=JδLT.
Consequently:
δbi=−J−1(δJ)bi=−δLTbi.
Since δLT=∑j=14bjδujT,
δbi=−j=1∑4bjbiTδuj.
The contribution to the variation of the internal nodal force associated with node i at one Gauss point is:
δfi(B)=hdetJσδbi.
Therefore:
δfi(B)=−hdetJj=1∑4(σbj)biTδuj,
and the corresponding 2×2 nodal block is:
kelB,ij=−hdetJ(σbj)biT.
This contribution represents the variation of the spatial interpolation derivatives with the current geometry.
6.3. Variation of the Jacobian determinant
For any nonsingular matrix, the determinant variation identity is:
δ(detJ)=detJtr(J−1δJ).
Using δJ=JδLT gives
J−1δJ=δLT,
and hence:
δ(detJ)=detJtr(δL).
The trace can be expressed in terms of virtual nodal displacements:
tr(δL)=j=1∑4bjTδuj=∂x∂(δu)+∂y∂(δv).
Consequently:
δ(detJ)=detJj=1∑4bjTδuj.
For node i, the corresponding internal-force variation is:
δfi(J)=hσbiδ(detJ),
so that:
δfi(J)=hdetJj=1∑4(σbi)bjTδuj.
The associated nodal tangent block is therefore
kelJ,ij=hdetJ(σbi)bjT.
6.4. Geometric contribution to the tangent stiffness matrix
Combining the variations of B and detJ, we obtain the geometric contribution to the tangent stiffness matrix. The 2×2 block associated with nodes i and j is:
Here p indexes the Gauss points, and i,j index the nodes. Each block is 2×2; arranging the sixteen blocks yields the 8×8 matrix kelg. The index g denotes geometric; the index p denotes the Gauss point.
The corresponding MATLAB implementation evaluates the geometric contribution at one Gauss point, before multiplying by the current thickness and Jacobian determinant:
kg=zeros(8);
for i=1:4
ii=2*i-1:2*i;
bi=b(:,i);
sbi=sigma*bi;
for j=1:4
jj=2*j-1:2*j;
bj=b(:,j);
kg(ii,jj)=sbi*bj'-(sigma*bj)*bi';
end
end
The geometric contribution is obtained by direct linearization of the adopted spatial internal-force vector. The resulting tangent matrix is not necessarily symmetric; it should not be artificially symmetrized.
6.5. Variation of the element thickness
The element thickness changes with deformation. For the adopted approximate plane-stress thickness update,
h1=h[1−1−νν(εx+εy)],
where h is the initial thickness. Since h1 depends on the total Euler–Almansi strains, it also depends on the nodal displacements. If the tangent is to remain consistent with the internal-force expression actually implemented, this dependence must be differentiated:
δh1=−h1−νν(δεx+δεy).
Using δεv=BQδuel,
∂uel∂h1=−h1−νν[BQ(1,:)+BQ(2,:)].
This derivative is a 1×8 row matrix. The thickness contribution to the tangent is:
kelh=p=1∑4detJp(BpTσv,p)∂uel∂hp.
The product of an 8×1 column matrix and a 1×8 row matrix gives an 8×8 matrix, as required. In MATLAB,
dth=-th*nu/(1-nu)*(BQ(1,:)+BQ(2,:));
The thickness formula remains approximate. Nevertheless, the term kelh makes the tangent consistent with this particular thickness update.
6.6. Complete tangent stiffness matrix
The complete element tangent stiffness matrix is:
kel=kelm+kelg+kelh.
The material part describes the variation of the Cauchy stresses through the total Euler–Almansi strains. The geometric part describes the variations of the spatial interpolation derivatives and the Jacobian determinant. The thickness part describes the dependence of h1 on the current displacements.
For one Gauss point, the MATLAB assembly statement is
kel=kel+detJ*(th1*(km+kg)+(B'*sig)*dth);
The assembled matrix is the consistent derivative of the adopted discrete internal-force vector, including its approximate constitutive and thickness laws.
07
Numerical example
The cantilever beam considered in Example 2 of Section 15.2 is analyzed again. Its length is L=400mm, and the rectangular cross-section is 20×5mm. The material properties are E=1000MPa and ν=0.3. A vertical force F=100N is applied at the free end. The mesh contains 120×16=1920 Q4 elements and 2057 nodes.
The computation uses total Euler–Almansi strains and the consistent tangent derived above, with the current geometry updated during each Newton iteration. Converged displacement extrema are approximately
vmin=−289.29mm,umax=2.45mm,umin=−159.34mm.
For the tested numbers of load steps, the converged values agree within the numerical tolerance. For comparison, the ANSYS results are:
The agreement is satisfactory, especially for the maximum downward displacement and the minimum horizontal displacement.
Comparison with the preceding formulations
The following results compare the original formulations of Sections 15.1, 15.2 and 15.3. All displacements are in millimeters. The corresponding results for the consistent Section 15.5 formulation are approximately
vmin=−289.29mmumin=−159.34mmumax=2.45mm
Load steps
Section
vmin
umin
umax
Total iterations
3
15.1
−275.9634
−142.0109
2.2672
88
3
15.2
−288.4194
−157.5978
2.4791
169
3
15.3
−288.4280
−157.6271
2.4790
202
5
15.1
−294.5759
−169.2246
2.3310
65
5
15.2
−289.0151
−158.8091
2.4644
100
5
15.3
−289.0151
−158.8091
2.4644
100
10
15.1
−294.0722
−168.1571
2.3112
85
10
15.2
−289.2207
−159.2074
2.4561
104
10
15.3
−289.2207
−159.2074
2.4561
104
20
15.1
−291.5248
−163.8582
2.2918
134
20
15.2
−289.2711
−159.3015
2.4525
150
20
15.3
−289.2711
−159.3014
2.4525
150
50
15.1
−289.9523
−161.2587
2.2807
264
50
15.2
−289.2876
−159.3311
2.4505
300
50
15.3
−289.2876
−159.3310
2.4505
300
The original Section 15.1 formulation is particularly sensitive to the number of load steps when only a few steps are used. Its results approach the consistent solution as the increments become smaller. The formulations in Sections 15.2 and 15.3 give much closer results even for a small number of steps; in this example they are practically identical for five or more steps and approach the same converged displacements.
The direct formulation of the present section produces practically identical converged displacements for different load-step subdivisions within the prescribed numerical tolerance. This observed load-step independence is not, by itself, a proof of tangent consistency; consistency follows from the linearization of the implemented internal-force vector.
08
MATLAB programs
The accompanying MATLAB programs implement the Updated Lagrangian formulation for four-node isoparametric Q4 elements with 2×2 Gauss integration. The main program performs the load-step and Newton-iteration loops and updates the nodal coordinates after each displacement correction.
The script stiff_cons.m calculates the element internal forces and assembles the global tangent stiffness matrix. At each Gauss point, the initial and current Jacobian matrices are used to compute the total deformation gradient, the Euler–Almansi strains and the Cauchy stresses. The material, geometric and thickness terms are then combined.
The files plot_disp.m and plot_stress.m display the nodal displacements and the computed stresses.
[1] Bonet, J., and Wood, R. D. (2008). Nonlinear Continuum Mechanics for Finite Element Analysis. 2nd ed. Cambridge University Press. See especially Chapters 4, 8 and 9 for finite-deformation kinematics, linearized equilibrium and finite-element discretization.