16.3. Numerical evaluation of the tangent stiffness matrix

Main idea

This chapter was written by me, Chiara, based on an idea by Mircea — an idea I liked so much that I simply would not let him write it himself.

StrainSmall

DisplacementLarge

MaterialLinear elastic

Stress statePlane stress

Sections 16.1 and 16.2 use different approximations of the element tangent stiffness matrix. In the present section, the tangent matrix is obtained directly by central numerical differentiation of the element internal nodal force vector:

\boldsymbol{k}_{el}=\frac{\partial\boldsymbol{f}_{el}}{\partial\boldsymbol{u}_{el}}

The nodal displacement vector of the 4-node quadrilateral element is:

\boldsymbol{u}_{el}=\left(u_1\;v_1\;u_2\;v_2\;u_3\;v_3\;u_4\;v_4\right)^T

Central finite differences

Let $\delta$ denote the numerical differentiation perturbation. For each element degree of freedom $i$ , two perturbed displacement vectors are introduced:

\boldsymbol{u}_{el}^{+}=\boldsymbol{u}_{el}+\delta\boldsymbol{e}_i

\boldsymbol{u}_{el}^{-}=\boldsymbol{u}_{el}-\delta\boldsymbol{e}_i

Here, $\boldsymbol{e}_i$ is the eight-component unit vector associated with degree of freedom $i$ : its $i$ th component is one and all other components are zero. Column $i$ of the element tangent matrix is approximated by:

\boxed{\boldsymbol{k}_{el}(:,i)\approx\frac{\boldsymbol{f}_{el}\!\left(\boldsymbol{u}_{el}+\delta\boldsymbol{e}_i\right)-\boldsymbol{f}_{el}\!\left(\boldsymbol{u}_{el}-\delta\boldsymbol{e}_i\right)}{2\delta}}

For example, for $i=3$ :

\boldsymbol{e}_3=\left(0\;0\;1\;0\;0\;0\;0\;0\right)^T

The MATLAB variable corresponding to $\delta$ is dd.

Element internal-force reevaluation

For each $+\delta$ or $-\delta$ perturbation, the $8\times1$ element internal nodal force vector $\boldsymbol{f}_{el}$ is recomputed.

The rot routine is called again for every perturbation. It computes the rigid-body rotation and the corotated nodal-displacement vector:

\widetilde{\boldsymbol{u}}_{el}=(\boldsymbol{I}_2-\boldsymbol{R})\boldsymbol{x}_{el0}+\boldsymbol{u}_{el}

where $\boldsymbol{x}_{el0}$ contains the nodal coordinates of the initial element configuration and $\boldsymbol{u}_{el}$ contains the current nodal displacements.

Consequently, each perturbed state includes the updated rigid-body rotation, corotational configuration, corotated nodal displacements, Jacobian matrix, strain-displacement matrix, strains, stresses, current thickness, and internal nodal force vector. The calculation follows the sequence:

\boldsymbol{u}_{el}\longrightarrow\boldsymbol{R}\longrightarrow\widetilde{\boldsymbol{u}}_{el}\longrightarrow\boldsymbol{J},\boldsymbol{B}\longrightarrow\boldsymbol{\varepsilon}\longrightarrow\boldsymbol{\sigma}\longrightarrow\boldsymbol{f}_{el}

MATLAB implementation

The calculation of the element internal nodal force vector is placed in the separate subprogram fel_el. Therefore, the same sequence is used for the actual state and for the $+\delta$ and $-\delta$ states without repeating the code.

No analytical differentiation subprogram such as deltaB, used in Section 16.2, is required.

In the following excerpt, nod contains the four node numbers of the current element; u and v are the global nodal-displacement vectors in the $x$ - and $y$ -directions; and inod is the local node position associated with the perturbed degree of freedom:

The relevant part of stiff.m is shown below.

    % Internal nodal force vector in the current configuration:
    istore=1;                % store strains and stresses
    fel_el
    fel=felt;

    % Central numerical differentiation of the element internal force:
    istore=0;                % do not store perturbed strains and stresses
    for ipert=1:8
        inod=ceil(ipert/2);

        % +dd perturbation:
        if rem(ipert,2)==1
            u(nod(inod))=u(nod(inod))+dd;
        else
            v(nod(inod))=v(nod(inod))+dd;
        end
        fel_el
        felp=felt;

        % Change from +dd to -dd:
        if rem(ipert,2)==1
            u(nod(inod))=u(nod(inod))-2*dd;
        else
            v(nod(inod))=v(nod(inod))-2*dd;
        end
        fel_el
        felm=felt;

        % Restore the unperturbed displacement:
        if rem(ipert,2)==1
            u(nod(inod))=u(nod(inod))+dd;
        else
            v(nod(inod))=v(nod(inod))+dd;
        end

        % Column ipert of the tangent stiffness matrix:
        kel(:,ipert)=(felp-felm)/(2*dd);
    end

The relevant part of fel_el.m is shown below.

    strn=B*sel;
    sg=DHooke*strn;
    th1=(1-nu/(1-nu)*(strn(1)+strn(2)))*th;
    felt=felt+th1*(B'*sg)*detJ;

Sensitivity to the perturbation δ

Example 1 of Section 16.2 was used to investigate the influence of the numerical differentiation perturbation $\delta$ .

The computations used:

n_{\mathrm{step}}=10\qquad \mathrm{tol}=10^{-5}

$\delta$	CPU time (s)	Total iterations	Mean $\eta_{\mathrm{sym}}$
$10^{-3}$	44.78	67	$1.3096\times10^{-2}$
$10^{-5}$	44.396	67	$1.3096\times10^{-2}$
$10^{-6}$	45.55	67	$1.3096\times10^{-2}$
$10^{-7}$	46.53	67	$1.3096\times10^{-2}$
$10^{-8}$	44.866	67	$1.3096\times10^{-2}$
$10^{-9}$	43.774	67	$1.3096\times10^{-2}$
$10^{-10}$	58.13	69	$1.3100\times10^{-2}$

Convergence is practically unchanged for $10^{-3}$ through $10^{-9}$ . At $\delta=10^{-10}$ , a slight deterioration appears. For $\delta=10^{-11}$ , the iterative process no longer converges. The value adopted in the downloadable MATLAB program is:

\boxed{\delta=10^{-5}}

Symmetry criterion

The relative asymmetry of an element tangent matrix is measured by:

\boxed{\eta_{\mathrm{sym}}=\frac{\left\|\boldsymbol{k}_{el}-\boldsymbol{k}_{el}^{T}\right\|_F}{\left\|\boldsymbol{k}_{el}\right\|_F}}

Here, the Frobenius norm of a matrix $\boldsymbol{A}$ is:

\left\|\boldsymbol{A}\right\|_F=\sqrt{\sum_i\sum_j a_{ij}^{2}}

The symmetry indicator therefore measures the relative magnitude of the antisymmetric part of $\boldsymbol{k}_{el}$ . A value close to zero indicates an approximately symmetric matrix.

For Section 16.3, the mean value is:

\overline{\eta}_{\mathrm{sym}}\approx1.31\times10^{-2}

No artificial symmetrization is used. Artificial symmetrization was tested. It made convergence somewhat slower, and for $\delta=10^{-10}$ the symmetrized tangent no longer converged.

Comparison with Sections 16.1 and 16.2

All three formulations use the same element internal nodal force vector $\boldsymbol{f}_{el}$ and therefore give the same final equilibrium results. The final stresses, displacements, and displacement derivatives are identical in Sections 16.1, 16.2, and 16.3.

For Example 1, the post-processed nodal axial stresses are:

\sigma_{\max}=50.519\,\mathrm{MPa}

\sigma_{\min}=-50.711\,\mathrm{MPa}

Section	Tangent	Total iterations	CPU time (s)	Mean $\eta_{\mathrm{sym}}$
16.1	Simple approximation; $s_r=0.5$	384	≈ 40.18	$6.87\times10^{-17}$
16.2	Improved analytical approximation	68	≈ 14	$1.9060\times10^{-2}$
16.3	Central numerical differentiation, $\delta=10^{-5}$	67	≈ 44.4	$1.310\times10^{-2}$

In Section 16.1, symmetry follows directly from the form:

\boldsymbol{B}^{T}\boldsymbol{D}\boldsymbol{B}

because $\boldsymbol{D}$ is symmetric. The measured asymmetry is therefore essentially roundoff error. Section 16.2 is the fastest of the three implementations. Section 16.3 gives essentially the same convergence rate without analytical differentiation of the tangent matrix.

Main conclusion and MATLAB package

\boxed{\begin{gathered}\text{The internal nodal force vector}\\\text{determines the equilibrium solution.}\end{gathered}}

\boxed{\begin{gathered}\text{The tangent stiffness matrix}\\\text{determines mainly how efficiently}\\\text{that solution is reached.}\end{gathered}}

The identical final results are expected because the same $\boldsymbol{f}_{el}$ is used in all three formulations.

The MATLAB package containing the central numerical-difference tangent and the two examples can be downloaded below.

MMATLAB packageDOWNLOAD ZIP