UL formulation with multiplicative update of the deformation gradient

01

From additive strains to multiplicative kinematics

Section 17.3 improved the tangent stiffness approximation and accounted for the different configurations associated with the stored strain and the strain increment. However, the accumulation of Hencky or Biot strain remained an approximation. These strain tensors are not generally additive under successive finite deformations:

The present section changes the quantity transferred between load steps. Instead of adding finite strains, the deformation gradients are composed multiplicatively:

Two implementations are considered. The first, stiff, evaluates the total deformation gradient directly from the initial configuration. The second, stiff_mult, stores the preceding total deformation gradient and updates it with the incremental deformation gradient. The second implementation is the main development of this section.

02

Direct total evaluation: stiff

Let $\boldsymbol{J}_0$ be the Jacobian matrix of the initial element configuration and $\boldsymbol{J}_t$ the Jacobian matrix of the current trial configuration. At each Gauss point, the total deformation gradient is evaluated as:

The right Cauchy–Green tensor and the right stretch tensor then follow:

Using the definitions established in Section 17.1, the total strain is:

No strain from the preceding load step is added. Every trial displacement vector defines a trial geometry, from which the complete strain tensor is recalculated.

J0=H*[xel0 yel0];
Jt=H*[xelt yelt];
FF=(J0\Jt)';
C=FF'*FF;
U=sqrtm(C);

Although the nonlinear solution proceeds through successive load steps, the total strain evaluation always refers directly to the initial configuration. Therefore, this first implementation should not be regarded as a strictly Updated Lagrangian formulation. It works very well numerically, but its reference configuration must be distinguished from that of stiff_mult.

03

Multiplicative update: stiff_mult

The analogy with Section 17.3 is useful. The same load-step and Newton-iteration structure is retained, but a different state quantity is stored and transferred:

At the beginning of the analysis, the total deformation gradient is the identity matrix at every Gauss point:

During load step $i$, $\boldsymbol{J}$ refers to the converged configuration at the beginning of that step, while $\boldsymbol{J}_t$ refers to the current trial configuration. The incremental deformation gradient is:

The total deformation gradient is updated by composition:

where $i$ denotes the current load step, $\boldsymbol{F}_{\rm tot}^{\,i-1}$ is the total deformation gradient stored at the end of the preceding converged load step, and $\boldsymbol{F}_{\rm inc}^{\,i}$ is the incremental deformation gradient evaluated during the current load step.

The order of multiplication follows the sequence of the two mappings: first from the initial configuration to the preceding converged configuration, then from that configuration to the current trial configuration.

In main.m, the storage array has two matrix indices, one Gauss-point index, and one load-step index. Four Gauss points are used per element. The extra load-step slot represents the initial state:

Fprev=zeros(2,2,4*nel,nstep+1);
Fprev(1,1,:,:)=1;
Fprev(2,2,:,:)=1;

For Gauss point ig during load step istep, the preceding converged value is read from slot istep. The current trial value is written to slot istep+1:

Finc=(J\Jt)';
Fp=Fprev(:,:,ig,istep);
FF=Finc*Fp;
Fprev(:,:,ig,istep+1)=FF;

The preceding slot remains unchanged during the Newton iterations. Only the current trial slot is overwritten. When the load step converges, this stored total deformation gradient becomes the preceding value for the next step. The program retains the same convergence tolerance and end-of-iteration storage convention as the other examples.

04

Recalculation of the total strain

After the multiplicative update, the strain is not incremented. It is recalculated from the complete deformation gradient:

The sequence is therefore:

No additive strain update and no reference-configuration correction from Section 17.3 are required. Hencky and Biot remain the same strain measures introduced in Section 17.1; the change concerns how the total deformation is determined.

Both subroutines allow selection of the strain measure through istrain:

istrain=1;   % 1 - Hencky strain
             % 2 - Biot strain

if istrain==1
    ee=logm(U);
    Uic=Ui^2;          % Hencky strain tensor
elseif istrain==2
    ee=U-eye(2);
    Uic=Ui;            % Biot strain tensor
else
    error('Unknown strain tensor option')
end

05

The strain–displacement matrix

In Section 17.3, the strain–displacement matrix $\boldsymbol{B}$ is based on the incremental deformation gradient. In stiff_mult, the total deformation gradient is the product:

During the Newton iterations of the current load step, $\boldsymbol{F}_{\rm prev}$ is fixed. Therefore, a variation of the trial nodal displacements gives:

This additional factor must be included when constructing the elementary matrices used to form $\boldsymbol{B}$. This is why B_matrix is replaced by B_matrix_mult in the multiplicative implementation. One representative line illustrates the change:

B1=FF'*[b(1,1) b(2,1);0 0];

B1=FF'*([b(1,1) b(2,1);0 0]*Fp);

The same postmultiplication by Fp is applied to the other seven elementary matrices, one for each remaining nodal degree of freedom. This is the essential difference between B_matrix and B_matrix_mult associated with the multiplicative update.

In these expressions, FF contains the total deformation gradient and Fp contains the preceding converged total deformation gradient. The eight matrices are then combined to form the three engineering-strain rows of $\boldsymbol{B}$.

06

Updated Lagrangian derivatives

The configuration used to evaluate the displacement derivatives is equally important. In stiff_mult, the derivatives refer to the converged configuration at the beginning of the current load step:

Here $\boldsymbol{H}$ contains the derivatives of the Q4 shape functions with respect to the natural coordinates. The implementation is:

J1=inv(J);
b=J1*H;

By contrast, stiff uses the initial Jacobian:

The first implementation therefore evaluates both the total deformation gradient and these derivatives from the initial configuration. In stiff_mult, the reference for the incremental derivatives is updated after each converged load step. Together with the multiplicative state transfer, this gives the second implementation its intended Updated Lagrangian character.

07

Geometric matrices and tangent approximation

The geometric matrices retain the same form as in Section 17.3:

These are symmetric $8\times8$ matrices formed from the displacement-derivative rows. The element tangent keeps the material contribution $\boldsymbol{B}^T\boldsymbol{D}\boldsymbol{B}$ and the stress-weighted geometric contributions involving $\boldsymbol{G}_x,\boldsymbol{G}_y,\boldsymbol{G}_{xy}$.

The strain–displacement matrix is modified consistently with the multiplicative update of the deformation gradient, while the geometric matrices are retained in the same approximate form used in Section 17.3. The numerical comparison below shows a slightly lower convergence rate than for the direct implementation, but improved consistency of the state update and excellent agreement with the ANSYS solution.

The retained geometric terms must not be interpreted as the exact consistent tangent of the new multiplicative formulation. This remains a tangent approximation, within the small-strain constitutive and thickness assumptions of the present examples.

08

Numerical comparison

The cantilever is the one considered in Section 17.3: length $L=400\ \mathrm{mm}$, rectangular cross-section $5\times20\ \mathrm{mm}$, Young’s modulus $E=1000\ \mathrm{MPa}$, and Poisson’s ratio $\nu=0.3$. The left end is fixed and the free end carries a vertical force of magnitude $100\ \mathrm{N}$.

The same mesh is used for all three MATLAB formulations:

The table compares Hencky-strain calculations with tol=1e-5 and a maximum of 100 Newton iterations per load step. The displacement columns give the maximum absolute nodal displacement components. Timings are local MATLAB R2026a measurements of the solution loop, without plotting; they are machine-dependent.

The ANSYS reference displacement magnitudes are:

Both Section 17.4 implementations converge with only two load steps and show much less sensitivity to load-step size than Section 17.3. The multiplicative implementation may require slightly more iterations after the modification of $\boldsymbol{B}$, but its displacement results remain in excellent agreement with ANSYS.

The slightly smaller displacement error is not the main result. More important is the consistent transfer of the total deformation gradient from one converged state to the next, without treating finite Hencky or Biot strain as additive. These observations concern the present benchmark and are not universal convergence guarantees.

The Biot-strain implementation was also tested and remained stable, requiring approximately the same number of Newton iterations as the corresponding Hencky-strain formulation.

09

What is transferred between load steps?

The change can be summarized by placing the stored quantities and their updates side by side:

The essential improvement is not a different strain measure. It is a different method of transferring the state from one converged load step to the next. Once the total deformation gradient is known, the complete Hencky or Biot strain is evaluated again, rather than assembled from a sequence of strain increments.

10

Q8 finite element

The same multiplicative Updated Lagrangian formulation also applies to the isoparametric Q8 element. The theoretical development remains the same; the main differences are the interpolation functions, sixteen element displacement degrees of freedom instead of eight, the construction and dimensions of the strain–displacement matrix $\boldsymbol B$, the geometric matrices, and $3\times3$ Gauss integration. The linear Q8 formulation is presented in Section 9.6, and the Q8 Total Lagrangian formulation in Section 12.3.

The preceding example is also solved with stiff_mult using Q8 elements. The comparison below uses ten load steps and a Q8 mesh with a number of nodes of the same order as the Q4 mesh. Mesh dimensions are given as height × length:

Each element formulation is compared with its corresponding ANSYS reference:

For these meshes, Q8 uses only one quarter of the number of elements, with a node count of the same order as Q4, and its CPU time is slightly lower. The Q8 displacements are in very good agreement with their ANSYS reference: the differences are approximately 0.029% for $u$ and 0.018% for $v$. These observations apply to the present example and tested meshes, not as general performance claims. The Q4 and Q8 values should each be judged against their respective ANSYS reference, rather than requiring one discretization to reproduce the other.

11

MATLAB programs

main.m calls stiff_mult by default. Replacing this call by stiff selects the direct total implementation. Set istrain=1 for Hencky strain or istrain=2 for Biot strain at the beginning of the selected stiffness subroutine.

The MATLAB package for both implementations can be downloaded below.

12

Conclusions

1. stiff eliminates approximate strain additivity by direct evaluation of the total deformation gradient, but is not strictly Updated Lagrangian because it refers directly to the initial configuration.
2. stiff_mult stores the converged total deformation gradient at every Gauss point. During the next load step, the incremental deformation gradient is composed multiplicatively with this stored value.
3. The total Hencky or Biot strain is recalculated from the resulting right stretch tensor. No additive strain accumulation is required.
4. B_matrix_mult incorporates the variation $\delta\boldsymbol{F}_{\rm tot}=\delta\boldsymbol{F}_{\rm inc}\boldsymbol{F}_{\rm prev}$, with the preceding converged deformation gradient held fixed during the Newton iterations.
5. The geometric matrices remain approximate. The improvement in the state update does not make the complete tangent an exact analytical derivative.
6. For the cantilever considered here, both implementations are robust with few load steps and give displacement results close to ANSYS. The main benefit is conceptual consistency, rather than a small change in the final displacement error.