17.5. Direct numerical evaluation of the strain–displacement and tangent stiffness matrices

Motivation

Large displacements

Small strains

Linear elastic material

Section 17.4 introduced a multiplicative update of the deformation gradient to avoid the approximate additive accumulation of Hencky or Biot strains. Section 17.5 retains this update:

\boxed{\boldsymbol F_{\rm tot}^{\,i}=\boldsymbol F_{\rm inc}^{\,i}\boldsymbol F_{\rm prev}^{\,i-1}}

At the beginning of each load step, the preceding converged configuration provides the reference for the increment. The Jacobian matrices $\boldsymbol J$ and $\boldsymbol J_t$ correspond to this configuration and the current trial configuration, respectively. The incremental and total deformation gradients are:

\boldsymbol F_{\rm inc}=(\boldsymbol J^{-1}\boldsymbol J_t)^T,\qquad \boldsymbol F_{\rm tot}=\boldsymbol F_{\rm inc}\boldsymbol F_{\rm prev}

The total Hencky or Biot strain is recalculated from $\boldsymbol F_{\rm tot}$ at each trial evaluation. Finite strain tensors are therefore not added from one load step to the next. The preceding converged $\boldsymbol F_{\rm prev}$ remains fixed during the numerical perturbations.

Can the strain–displacement matrix and the element tangent stiffness matrix be obtained directly by numerical differentiation, without deriving their analytical expressions? The computational sequence considered here is:

\boxed{\boldsymbol F_{\rm prev}\rightarrow\boldsymbol F_{\rm inc}\rightarrow\boldsymbol F_{\rm tot}\rightarrow\text{total strain}\rightarrow\boldsymbol B_{\rm num}\rightarrow\boldsymbol f_{\rm el}\rightarrow\boldsymbol k_{\rm el,num}}

The total strain is obtained from the right stretch tensor, using the definitions established earlier:

\boldsymbol C=\boldsymbol F_{\rm tot}^T\boldsymbol F_{\rm tot},\qquad \boldsymbol U=\sqrt{\boldsymbol C},\qquad \boldsymbol E_H=\ln\boldsymbol U,\qquad \boldsymbol E_B=\boldsymbol U-\boldsymbol I_2

Numerical evaluation of the strain–displacement matrix

A Q4 element has eight nodal displacement degrees of freedom. Let $\boldsymbol u_{\rm el}$ be the element displacement vector and $\boldsymbol e_j$ the unit vector associated with degree of freedom $j$ . For example, for the degrees of freedom $u_1$ and $v_2$ :

\boldsymbol e_1=\left\{\begin{array}{c}1\\0\\0\\0\\0\\0\\0\\0\end{array}\right\},\qquad\boldsymbol e_4=\left\{\begin{array}{c}0\\0\\0\\1\\0\\0\\0\\0\end{array}\right\}

The unit entry occupies the position of the perturbed degree of freedom.

\boldsymbol u_{\rm el}=\{u_1,v_1,u_2,v_2,u_3,v_3,u_4,v_4\}^{T}

\boxed{\boldsymbol B=\frac{\partial\boldsymbol e}{\partial\boldsymbol u_{\rm el}}}

Each column of $\boldsymbol B$ is evaluated by central differences:

\boxed{\boldsymbol B_j\simeq\frac{\boldsymbol e(\boldsymbol u_{\rm el}+ddB\,\boldsymbol e_j)-\boldsymbol e(\boldsymbol u_{\rm el}-ddB\,\boldsymbol e_j)}{2\,ddB}}

Here $\boldsymbol e$ contains the engineering strain components returned by the total-strain algorithm at the Gauss point. The parameter ddB is the displacement perturbation used for this first differentiation.

\boxed{\boldsymbol u_{\rm el}\in\mathbb R^{8\times1},\qquad\boldsymbol e\in\mathbb R^{3\times1},\qquad\boldsymbol B\in\mathbb R^{3\times8},\qquad\boldsymbol k_{\rm el}\in\mathbb R^{8\times8}}

Accordingly, $\boldsymbol B$ maps the eight element displacement degrees of freedom to the three in-plane strain components, while $\boldsymbol k_{\rm el}$ relates the eight element force and displacement components.

For each perturbed trial configuration, the multiplicative update includes the MATLAB sequence:

Multiplicative deformation-gradient update

Finc=(J\Jt)';
Ftot=Finc*Fp;

The entire strain algorithm is differentiated numerically. Each positive and negative perturbation passes through the trial Jacobian, the incremental and total deformation gradients, and the stretch and strain calculations. Thus, the intermediate operations do not need to be differentiated analytically. Both evaluations use the same stored $\boldsymbol F_{\rm prev}$ ; the perturbations must not alter the converged history.

Numerical evaluation of the tangent stiffness matrix

Once $\boldsymbol B_{\rm num}$ is available, the element internal-force vector $\boldsymbol f_{\rm el}$ is evaluated by the implemented Gauss-point calculation. A second central differentiation then gives each column of the element tangent stiffness matrix:

\boxed{\boldsymbol k_{{\rm el},j}\simeq\frac{\boldsymbol f_{\rm el}(\boldsymbol u_{\rm el}+ddK\,\boldsymbol e_j)-\boldsymbol f_{\rm el}(\boldsymbol u_{\rm el}-ddK\,\boldsymbol e_j)}{2\,ddK}}

The parameter ddK is the displacement perturbation for the second differentiation. Every perturbed internal-force evaluation includes a fresh numerical calculation of $\boldsymbol B$ . Consequently:

\boxed{\boldsymbol k_{\rm el}=\frac{\partial\boldsymbol f_{\rm el}}{\partial\boldsymbol u_{\rm el}}}

The geometric matrices $\boldsymbol G_x$ , $\boldsymbol G_y$ , and $\boldsymbol G_{xy}$ are no longer required explicitly. The tangent follows from differentiation of the complete implemented internal-force calculation, including the displacement dependence contained in that calculation.

Computational cost

For each Gauss point, calculating the numerical $\boldsymbol B$ requires one unperturbed strain evaluation, eight positive perturbations and eight negative perturbations:

1+8+8=17

For the four Gauss points of a Q4 element, this gives the following number of strain evaluations per internal-force evaluation:

4\times17=68

The numerical tangent requires one unperturbed internal-force evaluation and sixteen perturbed evaluations. Because every one of these calculations repeats the numerical strain differentiation, the approximate number of strain evaluations for one element in one Newton iteration is:

17\times4\times17=1156

For the reduced mesh of 150 elements, the corresponding count is:

\boxed{150\times1156=173400}

strain evaluations per Newton iteration. Each evaluation repeats the kinematic and total-strain calculations. The two nested differentiations therefore replace analytical work with a very large amount of numerical work, even for this relatively small mesh. This explains the execution times reported below.

Influence of perturbation size

Central differences require a balance between two errors. Excessively small perturbations cause cancellation and round-off when nearly equal values are subtracted. Excessively large perturbations increase truncation error and make the derivative less local.

The following tests use $nstep=3$ (three load steps) and equal values of ddB and ddK:

$ddB=ddK$	Newton iterations	$u\ (\mathrm{mm})$	$v\ (\mathrm{mm})$	Time (s)
$10^{-6}$	No convergence	—	—	—
$10^{-5}$	34	142.8110	276.1843	106.5687
$10^{-4}$	30	142.8110	276.1843	98.7424
$10^{-3}$	30	142.8110	276.1844	97.5643
$10^{-2}$	31	142.8158	276.1888	86.3961
$10^{-1}$	No convergence	—	—	—

For $10^{-5}\le ddB=ddK\le10^{-3}$ , the displacement solutions are practically identical. The calculation with $10^{-2}$ still converges, but the difference in the displacement values becomes visible. The tests with $10^{-6}$ and $10^{-1}$ do not converge.

This robust interval applies to the present example. The perturbations are dimensional, so an appropriate value also depends on the displacement units and the scale of the problem. No single perturbation value is universally optimal.

Influence of the number of load steps

For the present example, one or two load steps are insufficient, whereas three load steps already provide convergence. For three or more load steps, the displacement solution varies smoothly and approaches a stable value as the load-step subdivision is refined.

The calculations with 3, 5, 10 and 20 load steps converge. Within the robust perturbation range, the load-step size has a stronger influence on nonlinear robustness than the perturbation size itself.

Comparison between Sections 17.4 and 17.5

Both formulations use the same reduced mesh of 186 nodes and 150 Q4 elements. The ANSYS reference displacement magnitudes are:

u_{\rm ANSYS}=142.795\ \mathrm{mm},\qquad v_{\rm ANSYS}=276.148\ \mathrm{mm}

Load steps	Method	$ddB=ddK$	Total Newton iterations	$u\ (\mathrm{mm})$	$v\ (\mathrm{mm})$	$t\ (\mathrm{s})$
1	17.4	—	18	142.8849	276.3104	0.791
1	17.5	$10^{-5}\ \&\ 10^{-3}$	No convergence	—	—	—
2	17.4	—	28	142.8083	276.1525	0.947
2	17.5	$10^{-5}\ \&\ 10^{-3}$	No convergence	—	—	—
3	17.4	—	37	142.8051	276.1334	1.23
3	17.5	$10^{-5}$	34	142.8110	276.1843	106.57
3	17.5	$10^{-3}$	30	142.8110	276.1844	97.56
5	17.4	—	54	142.8025	276.1211	1.78
5	17.5	$10^{-5}$	43	142.8081	276.1717	152.17
5	17.5	$10^{-3}$	43	142.8082	276.1718	171.94
10	17.4	—	104	142.8000	276.1126	3.49
10	17.5	$10^{-5}$	72	142.8054	276.1631	274.72
10	17.5	$10^{-3}$	66	142.8055	276.1631	265.15
20	17.4	—	197	142.7986	276.1086	6.080
20	17.5	$10^{-5}$	120	142.8039	276.1590	556.35
20	17.5	$10^{-3}$	106	142.8040	276.1590	429.46

The results obtained with the double numerical differentiation are very close to those of Section 17.4 and to the ANSYS reference solution. As the number of load steps increases, the Section 17.5 displacement solution approaches approximately $u=142.804$ mm and $v=276.159$ mm. The purpose of Section 17.5 is nevertheless not to improve the displacement accuracy, but to demonstrate that both the strain–displacement matrix and the tangent stiffness matrix can be constructed directly by numerical differentiation.

Newton convergence and execution time

The numerical tangent can reduce the number of Newton iterations, confirming its good consistency with the internal-force evaluation. However, each iteration is extremely expensive because of the two nested numerical differentiations.

For example, with 20 load steps, Section 17.4 requires 197 iterations and 6.080 s. Section 17.5 with $ddB=ddK=10^{-3}$ requires 106 iterations and 429.46 s. The reduction in iteration count is outweighed by the cost of each iteration:

\boxed{\text{fewer Newton iterations}\not\Rightarrow\text{shorter execution time}}

Practical assessment

The double numerical differentiation provides a tangent that is highly consistent with the internal-force evaluation and produces displacement results very close to both Section 17.4 and ANSYS. Its principal disadvantage is the very large computational cost. Therefore, the formulation remains mainly useful as a didactic and verification tool rather than as an efficient practical alternative to the analytical formulation.

MATLAB programs

The MATLAB package provides Hencky and Biot strain options. The strain_el.m routine evaluates the total strain using the multiplicative kinematics. The fel_el.m routine constructs the numerical strain–displacement matrix and the internal-force vector, and stiff.m evaluates the numerical tangent and assembles the global equations.

The MATLAB package can be downloaded below.

MMATLAB packageMultiplicative update with numerical B and numerical tangent stiffnessDownload ZIP

Conclusions

Section 17.5 retains the multiplicative deformation-gradient update of Section 17.4.
The total Hencky or Biot strain is recalculated from $\boldsymbol F_{\rm tot}$ .
The strain–displacement matrix $\boldsymbol B$ is obtained by numerical differentiation of the total strain.
The tangent stiffness matrix $\boldsymbol k_{\rm el}$ is obtained by a second numerical differentiation of the internal-force vector.
The matrices $\boldsymbol G_x$ , $\boldsymbol G_y$ , and $\boldsymbol G_{xy}$ are no longer required explicitly.
For the present example, three load steps already provide convergence; one or two load steps are insufficient.
For $nstep=3$ , the solution is practically unchanged for $10^{-5}\le ddB=ddK\le10^{-3}$ , while excessively small or excessively large perturbations may lead to loss of convergence.
The numerical tangent can reduce Newton iterations, but the large computational cost makes this primarily a didactic and verification tool.

\boxed{\text{Greater analytical simplicity can be purchased at a very high numerical cost.}}

The double numerical differentiation strategy used in this section was developed with the assistance of Chiara.