4. Stabilization

4.1. Principles

In the calculation of the quantities required, integral terms from the variational formulation appear (see § 3). The structure stiffness matrix will be written as follows for a given finite element \(e\):

(4.1)\[ {\ int} _ {{\ mathrm {\ Omega}}} _ {e}}}\ left (\ epsilon\ mathrm {:}\ delta\ sigma\ right) d\ mathrm {\ Omega}\ to {K}\ to {K} _ {e} _ {e}} = {\ e}} = {\ int}} = {\ int} = {\ int}} _ {\ mathrm {\ omega}}\ left ({B} ^ {T} CB\ right) d\ mathrm {\ omega}\]

With \(C\) the Hooke elasticity matrix and \(B\) the linearized deformation matrix. We must therefore use a numerical quadrature scheme to evaluate the integral. This diagram will contain \({n}_{\mathit{pt}}\) \({\xi }_{p=1,{n}_{\mathit{pt}}}\) coordinate points, each with a weight \({\omega }_{p=1,{n}_{\mathit{pt}}}\) such as:

(4.2)\[ {K} _ {e} = {\ int} _ {{\ mathrm {\ omega}}} _ {e}}\ left ({B} ^ {T} CB\ right) d\ mathrm {\ Omega} =\ sum _ {p=1} =\ sum _ {p=1}} ^ {p=1}} ^ {{n}} ^ {n}} _ {p}}} {\ omega} _ {p}) J ({\ xi} _ {p}) {B ({\ xi} _ {p})} ^ {T} CB ({\ xi} _ {p})\]

\(J\) is the Jacobian for transforming the element between parametric space and Cartesian space (see § 2.2). The choice of this diagram, its precision in particular, will depend on the nature of the integrand. For element COQUE_SOLIDE, we will use a 5-point Lobatto-type schema.

This diagram is too « poor » and it will underestimate the rigidity of the element in order to control the locks. But by proceeding in this way, deformation modes that are called « hourglass modes » are also revealed. The element is likely to change shape without expending energy, which does not correspond to physics. To minimize this problem, we will stabilize the element using appropriate technology.

The aim of the stabilization procedure is to correct the rank defect in the stiffness matrix resulting from the reduced integration scheme adopted. The reader is referred to the important works on this subject by [bib18] _, [bib19] _, [bib19] _, [bib14] _, [bib8] _, [bib8] _, [bib9] _, [bib9] _, [bib9] _, [bib9] _, [bib9] _, [bib9] _, [bib9] _, [bib9] _, [bib9] _, [bib9 In some shell-solid formulations, the Jacobian matrix and its inverse are only evaluated at the center of the element, see for example [bib21] _, [bib22] _, [bib23] _, [bib23] _. Such an approximation assumes that the real element can be represented by its equivalent parallelepiped. And it’s pretty accurate for fine meshes and minimally deformed ones.

However, for meshes that are highly deformed in the initial state, this hypothesis shows a lack of precision [bib10] _. In order to take into consideration the real shape of the element and to stabilize the stiffness matrix correctly, a polynomial decomposition of the inverse Jacobian matrix was constructed explicitly.

4.1.1. Stabilization of deformations

Moreover, for the transition between quantities such as deformation between parametric coordinates and real coordinates, given by the matrix \(T\), it is appropriate to use a polynomial decomposition consistent with that chosen for the Jacobian in particular. It was shown in [bib10] _ that very good precision is achieved by simply using the constant, linear terms from \(T\). Therefore, bilinear terms are ignored:

\[\]

: label: eq-81

Tapprox {T} ^ {0} +xi {T} +xi} +eta {T} ^ {eta} +zeta {T} +zeta {T} ^ {zeta} ^ {zeta}

We will use the same form of approximation on \({E}^{u}\) as on \({\widehat{E}}^{u}\):

\[\]

: label: eq-82

{E} ^ {u} =T {widehat {E}}} ^ {u} ^ {u}approx {E} ^ {u}} +zeta {E} ^ {u,zeta} + {zeta} + {zeta}} + {zeta} ^ {2} {2} {E} ^ {u,zeta} +zeta {e} ^ {u,zeta}} +zeta {e} ^ {u,zeta}} +zeta {e} ^ {u,zeta}} +zeta {e} ^ {u,zeta}} +zeta {e} ^ {u,zeta}} +zeta {e} ^ {u,zeta}} +zeta {e}} ^ {u,eta} +xieta {E} ^ {u,xieta} +xizeta {E} ^ {u,xizeta} +etazeta {E} ^ {u,etazeta} +etazeta} +etazeta} +etazeta} +etazeta {E} ^ {u,xizeta} +etazeta {E} ^ {u

Which we will divide into two components:

\[\]

: label: eq-83

{E} ^ {u} = {E} ^ {u,mathit {RI}}} + {E} ^ {u,mathit {STAB}}

{E} ^ {u,mathit {RI}}} = {E}} = {E} ^ {e} ^ {u,zeta}} + {zeta} + {zeta} ^ {2} {2} {E} ^ {u,zetazeta} + {zeta} + {zeta}

{E} ^ {u,mathit {STAB}}} =xi {E}}} =xi {E} ^ {u,eta} +xieta {E} ^ {u,xieta} ^ {u,xieta} +xixeta {eta}} +xizeta {E} ^ {u,xizeta} ^ {u,xizeta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta} ^ {u,xieta}

With the following explicit expressions:

(4.3)\[ \begin{align}\begin{aligned} {E} ^ {u\ mathrm {,0}}} = {T} ^ {0} {\ widehat {E}} ^ {u\ mathrm {,0}}}\\ {E} ^ {u,\ zeta} = {T} ^ {\ zeta} {\ zeta} {\ widehat {E}} ^ {0}} + {T} ^ {0} {\ widehat {E}}} ^ {u,\ zeta}} + {T} ^ {0} {\ widehat {E}}} ^ {u,\ zeta}\\ {E} ^ {u,\ zeta\ zeta} = {T} ^ {\ zeta} {\ widehat {E}} ^ {u,\ zeta} + {T} ^ {0} {\ widehat {E}}} ^ {u,\ zeta\ zeta}} + {T} ^ {0} {0} {\ widehat {E}} {0} {\ widehat {E}}} ^ {0} {\ widehat {E}}\\ {E} ^ {u,\ xi} = {T} ^ {0} {\ widehat {E}} ^ {u,\ xi} + {T} ^ {\ xi} {\ xi} {\ widehat {E}}} ^ {u,0}\\ {E} ^ {u,\ eta} = {T} ^ {0} {\ widehat {E}} ^ {u,\ eta} + {T} ^ {\ eta} {\ eta} {\ eta} {\ widehat {E}}} ^ {u,0}\\ {E} ^ {u,\ xi\ eta} = {T} ^ {0} {T} ^ {0} {\ widehat {E}} ^ {\ xi} + {T} ^ {\ xi} {\ widehat {E}}} {\ widehat {E}}} {\ widehat {E}} {\ widehat {E}} {\ xi} {\ widehat {E}} {\ widehat {E}} {\ widehat {E}} {\ widehat {E}}} {\ widehat {E}} {\ widehat {E}}} {\ widehat {E}} {\ widehat {E}}} {\ widehat {E}} {\\ {E} ^ {u,\ eta\ zeta} = {T} ^ {zeta}} = {T} ^ {0} {\ widehat {E}} ^ {\ zeta} {\ zeta} {\ widehat {E}}} ^ {e}}} ^ {u,\ eta}} ^ {u,\ zeta} + {T} ^ {\ eta} {\ widehat {E}}} ^ {u,\ zeta}} ^ {u,\ zeta}} ^ {u,\ zeta}\\ {E} ^ {u,\ xi\ zeta} = {T} ^ {zeta} = {T} ^ {0} {\ widehat {E}} ^ {\ xi} {\ xi} {\ widehat {E}}} ^ {e}}} ^ {u,\ zeta}} ^ {u,\ zeta} + {T} ^ {\ zeta} {\ widehat {E}}} ^ {u,\ xi}} ^ {u,\ xi}\end{aligned}\end{align} \]

The same decomposition will be used for the tensor of the linearized deformations and therefore for the gradient matrix:

(4.4)\[ \begin{align}\begin{aligned} \ widehat {B} = {\ widehat {B}}} ^ {\ mathit {RI}}} + {\ widehat {B}}} ^ {\ mathit {STAB}}} ^ {\ mathit {}}\\ \ textrm {and}\\ B=T\ widehat {B}\\ {B} ^ {\ mathit {RI}}} = {B} ^ {0} +\ zeta {B} ^ {\ zeta} + {\ zeta} ^ {2} {B} {B} ^ {B} ^ {\ zeta\ zeta}\\ {B} ^ {\ mathit {STAB}}} =\ xi {B}} ^ {\ B} ^ {\ xi} +\ eta {\ eta} +\ xi\ eta {B} ^ {\ xi\ eta} +\ xi\ zeta} +\ xi\ zeta {B} +\ xi\ zeta {B} ^ {\ eta\ zeta} +\ eta\ zeta {B} ^ {\ eta\ eta} +\ xi\ zeta {B} ^ {\ xi\ zeta} +\ xi\ zeta {B} ^ {\ xi\ zeta}\end{aligned}\end{align} \]

With:

(4.5)\[ \begin{align}\begin{aligned} {B} ^ {0} = {T} ^ {0} {\ widehat {B}} ^ {0}\\ {B} ^ {\ zeta} = {T} ^ {\ zeta} {\ zeta} {\ widehat {B}} ^ {0} + {T} ^ {0} {\ widehat {B}}} ^ {\ zeta}\\ {B} ^ {\ zeta\ zeta} = {T} ^ {\ zeta} {\ zeta} {\ widehat {B}} ^ {zeta} + {T} ^ {0} {\ widehat {B}}} ^ {\ zeta\ zeta}\\ {B} ^ {\ xi} = {T} ^ {0} {\ widehat {0}} {\ widehat {B}}} ^ {\ xi} {\ xi} {\ widehat {B}}} ^ {0}\\ {B} ^ {\ eta} = {T} ^ {0} {\ widehat {0}} {\ widehat {B}}} ^ {\ eta} {\ eta} {\ widehat {B}}} ^ {0}\\ {B} ^ {\ xi\ eta} = {T} ^ {0} {\ widehat {B}}} ^ {\ xi\ eta} + {\ xi} {\ xi} {\ widehat {B}}} ^ {\ eta}} ^ {\ eta} + {\ eta} ^ {\ eta} {\ eta} {\ widehat {B}} {\ xi}} ^ {\ xi}\\ {B} ^ {\ eta\ zeta} = {T} ^ {0} {T} ^ {0} {\ widehat {B}}} ^ {\ zeta} {\ zeta} {\ widehat {B}}} ^ {\ eta}} ^ {\ eta} + {\ eta}} ^ {\ eta} {\ eta}} ^ {\ zeta}} ^ {\ zeta}\\ {B} ^ {\ xi\ zeta} = {T} ^ {0} {T} ^ {0} {\ widehat {B}}} ^ {\ xi} {\ xi} {\ widehat {B}}} ^ {\ zeta}} ^ {\ zeta} + {\ zeta}} ^ {\ zeta} {\ zeta} {\ zeta} {\ zeta} {\ widehat {B}} ^ {\ xi} ^ {\ xi}\end{aligned}\end{align} \]

\({E}^{u,\mathit{STAB}}\) and \({B}^{\mathit{STAB}}\) represent the terms cancelled out in the reduced integration. The stabilization process consists in analytically restoring these terms in the stiffness matrix.

Additionally, to eliminate volumetric locking, the B-bar [bib24] _ approach is adopted. Therefore, the hourglass modes of deformation operators and stabilized deformation displacement operators are divided into volumetric and deviatory components:

(4.6)\[ {B} ^ {\ mathit {STAB}}} = {B} _ {\ mathit {dev}}} ^ {\ mathit {STAB}}\ left (\ xi,\ eta,\ zeta\ right) + {B} _ {\ mathit {vol}}}} ^ {\ mathit {vol}}} ^ {\ mathit {vol}}} ^ {\ mathit {STAB}}}\ left (\ mathrm {0,}\ zeta\ right) + {B} _ {\ mathit {vol}}}} {\ mathit {vol}}} ^ {\ mathit {vol}}} ^ {\ mathit {vol}}} ^ {\ mathit {}}}\ left (\ mathrm {0,}\ zeta\ right) + {B} _\]

There is no constant term in the expression for \({B}^{\mathit{STAB}}\), so:

(4.7)\[ {B} _ {\ mathit {vol}}} ^ {\ mathit {STAB}}}\ left (\ mathrm {0,}\ mathrm {0,} 0\ right) =0\]

Only the deviatoric part is retained. It can be shown that the contribution \({B}^{\xi \eta }\) (and \({E}^{u,\xi \eta }\)) is zero in the writing of stabilization.

4.1.2. Stabilization of stresses

In the expression of the energy of deformation of the structure, the separation between the contributions coming from reduced integration and stabilization is repeated:

(4.8)\[ W (E) = {W} ^ {\ mathit {RI}}} ({E} ^ {u,\ mathit {RI}}) + {W} ^ {\ mathit {STAB}}} ({E} {RI}}} ({E} {RI}}}) ({E} ^ {U,\ mathit {STAB}}}) + {W} ^ {}}} ({E} {RI}}} ({E} {RI}}})\]

Since the Piola-Kirchhoff stress tensor of the second kind is the derivative of this energy with respect to the Green-Lagrange deformations, we have:

(4.9)\[ S=\ frac {\ partial W} {\ partial E} =\ frac {\ partial {W}} ^ {\ mathit {RI}}} {\ partial {E} ^ {u,\ mathit {RI}}}}}\ frac {\ partial {RI}}}} {\ partial {W}}}\ frac {\ partial {W}}}}\ frac {\ partial {W}}}}\ frac {\ partial {W}}}}\ frac {\ partial {W}}}\ frac {\ partial {W}}}}\ frac {\ partial {W}}} {\ partial {W}}}\ frac {\ partial {W}}}} {\ partial {W}}} {\ partial {W}}} {\ partial {W}}} {STAB}}} {\ partial {E} ^ {u,\ mathit {STAB}}}}\ frac {\ partial {E} ^ {u,\ mathit {STAB}}}} {\ partial E}}\]

We find the separation between the two contributions (reduced integration and stabilization) for the stress tensor:

\[\]

: label: eq-91

S= {S} ^ {mathit {RI}}} ({E} ^ {u,mathit {RI}}) + {S} ^ {mathit {STAB}}} ({E} ^ {i} ^ {u,mathit {STAB}}})

The tensor \({S}^{\mathit{RI}}\) comes directly from the integration of the law of behavior while the tensor \({S}^{\mathit{STAB}}\) corresponds to the terms of constraints for stabilization.

The key point in stabilization is to find an optimal way to assess this stabilization constraint without integrating it so that the calculation time is minimized. However, as explained above, to avoid volumetric locking, only the deviatoric part of the Piola Kirchhoff stabilization tensor is retained. In doing so, the stabilization of the stress is given as:

\[\]

: label: eq-92

{S} ^ {mathit {STAB}}} = {C}} = {C} ^ {C} ^ {C} ^ {u,mathit {STAB}}}mathrm {:} {E} ^ {u,mathit {}} STAB

\({C}^{\mathit{STAB}}\) being the deviatory part of Saint-Venant Kirchhoff’s behavior:

(4.10)\[\begin{split} {C} ^ {\ mathit {STAB}}} = {\ mu}} ^ {\ mu} ^ {\ mathit {STAB}}\ left [\ begin {array} {cccc}\ frac {4} {3}} & -\ frac {2} {3}} & -\ frac {2} {3} & -\ frac {2} {3} &\ frac {4} {3} & -\ frac {2} {3} {3} & 0& 0& 0& 0\\ 0\ & 0\\ 0} {3} &\ frac {4} {3} &\ frac {4} {3} & 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\ end {array}\ right]\end{split}\]

With the shear modulus \({\mu }^{\mathit{STAB}}\), which, for an isotropic elastic material, is equal to:

(4.11)\[ {\ mu} _ {\ mathit {elas}}} ^ {\ mathit {STAB}} =\ frac {E} {2 (1+\nu)}\]

For non-elastic materials, the use of the elastic shear modulus leads to an overestimation of the stabilization stress. To overcome this problem, the secant module as defined in [bib20] _ is used:

(4.12)\[ 2 {\ mu} ^ {\ mathit {STAB}}} =\ sqrt {\ frac {{\ pi} _ {\ pi}} {{\ pi} _ {E}}}\]

With:

(4.13)\[ \begin{align}\begin{aligned} {\ pi} _ {S} =\ frac {1} {1} {2} {2} {S} _ {\ mathit {dev}}\ mathrm {:} {S} _ {\ mathit {dev}}\\ \ textrm {and}\\ {\ pi} _ {E} =\ frac {1} {1} {2} {2} {E}} _ {\ mathit {dev}} ^ {u}\ mathrm {:} {E} _ {\ mathit {dev}} {\ mathit {dev}}} ^ {u}\end{aligned}\end{align} \]

This way of evaluating the stabilization parameter is very practical because there is no need for coherent linearization of the stabilization matrix, which is very interesting in terms of calculation time. Moreover, the proposed definition of \({C}^{\mathit{STAB}}\) ensures effective and robust stabilization.

4.1.3. Stabilization of discrete operators

Using the previous definitions, we calculate the \({\mu }^{\mathit{STAB}}\) contribution in the following way. At each Gauss point \(g\) we calculate \({{\pi }_{S}}^{g}\) and \({{\pi }_{E}}^{g}\) and then consider the following two situations:

1. If:math: | {{pi} _ {E}}} ^ {e}}} ^ {g} |<mathrm {0.001} `) then:math: {mu} ^ {mathit {STAB} ^ {STAB}, g} = {K}} = {K} _ {K} _ {mathit {vol}}} ^ {g} . That is to say that the stabilization coefficient at the current Guass point:math: `g is equal to the volume part of the material behavior matrix if the deviatoric deformation is less than 0.1%;

2. If not \({\mu }^{\mathit{STAB},g}=\frac{1}{2}\sqrt{\frac{{\pi }_{S}^{g}}{{\pi }_{E}^{g}}}\);

Then the stabilization coefficient is summed on each Gauss point:

(4.14)\[ {\ mu} ^ {\ mathit {STAB}}} = {\ mu}} = {\ mu} ^ {\ mathit {STAB}, g} +\ frac {{\ omega} _ {g}} {g}} {8}\]

The stabilization of the stresses in the case of small deformations consists in calculating the various contributions in the following way (remember that the term \(\xi \eta\) does not contribute):

(4.15)\[ \begin{align}\begin{aligned} {\ epsilon} ^ {\ xi} = {B} ^ {\ xi} u\\ {\ epsilon} ^ {\ eta} = {B} = {B} ^ {\ eta} u` :math:`{\epsilon }^{\xi \zeta }={B}^{\xi \zeta }u` :math: `{\ epsilon} ^ {\ eta\ zeta,\ mathit {STAB}} = {B} ^ {\ eta\ zeta} u\\ {\ sigma} ^ {\ xi,\ mathit {STAB}}} = {C} ^ {\ mathit {STAB}} {\ epsilon} ^ {\ xi}\\ {\ sigma} ^ {\ eta,\ mathit {STAB}}} = {C}} = {C} ^ {\ mathit {STAB}} {\ epsilon} ^ {\ eta} `:math:`{\sigma }^{\xi \zeta ,\mathit{STAB}}={C}^{\mathit{STAB}}{\epsilon }^{\xi \zeta }` :math:` {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma STAB STAB\end{aligned}\end{align} \]

And in large deformations:

(4.16)\[ \begin{align}\begin{aligned} {\ sigma} ^ {\ xi,\ mathit {STAB}}} = {C} ^ {\ mathit {STAB}} {E} ^ {u,\ xi}\\ {\ sigma} ^ {\ eta,\ mathit {STAB}}} = {C}} = {C} ^ {\ mathit {STAB}} {E} ^ {u,\ eta} `:math:`{\sigma }^{\xi \zeta ,\mathit{STAB}}={C}^{\mathit{STAB}}{E}^{\xi \zeta }` :math:` {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {STAB STAB\end{aligned}\end{align} \]

Newton’s linearization process of the equilibrium equation produces two operators: the consistent tangent matrix and the internal force vector. Both will therefore have to integrate stabilization since they depend on constraints. For the tangent matrix, there are two contributions:

• the first concerns stabilization linked to the constraints resulting from the integration of the law of behavior (material matrix), i.e. matrix \({K}_{m}^{\mathit{uu}}\);

• the second corresponds to the large deformation model, it is the geometric part that comes from kinematic nonlinearity, i.e. matrix \({K}_{g}^{\mathit{uu}}\).

For the stabilization of the material tangent matrix, we have four contributions:

\[\]

: label: eq-100

{K} _ {m} ^ {mathit {STAB},xi} = {({B} ^ {xi})} ^ {T} {C} ^ {mathit {STAB}} ^ {mathit {}}} ({B} ^ {xi})

{K} _ {m} ^ {mathit {STAB},eta} = {({B} ^ {eta})} ^ {T} {C} ^ {mathit {STAB}} ^ {mathit {}}} ({B} ^ {eta})

{K} _ {m} ^ {mathit {STAB},xizeta} = {({B} ^ {xizeta})} ^ {T} {C} ^ {mathit {} ^ {mathit {}},xizeta} ^ {xizeta} ^ {xizeta}) STAB

{K} _ {m} ^ {mathit {STAB},etazeta} = {({B} ^ {etazeta})} ^ {T} {C} ^ {mathit {} ^ {mathit {}},etazeta} ^ {etazeta} ^ {etazeta}) STAB

This gives us the global material stabilization matrix with the Jacobian coordinate change ():

\[\]

: label: eq-101

{K} _ {m} ^ {mathit {STAB}}} =frac {8}}} {3} {J} _ {0}left ({K} _ {m} ^ {mathit {STAB} ^ {mathit {},xi},xi} + {K} _ {m}}}} =frac {8} {m}} ^ {mathit {}}} {0} =frac {8} {m} ^ {mathit {}}} ^ {K} _ {m}} ^ {mathit {}}}} =frac {8} {m} ^ {mathit {}}} ^ {K} _ {m}}} =frac {8} {m} ^ {0}left ({K} _ {m} ^ {mathit {STAB},xizeta} + {K} _ {m} ^ {mathit {STAB},etazeta}right) STAB

For the geometric part, as the contribution is the product of a (geometric) matrix with the constraints, it is also necessary to do it by a product for stabilization, for each node \(i\) and \(j\):

\[\]

: label: eq-102

{K} _ {g,mathit {ij}}} ^ {mathit {STAB}}} =frac {8} {3} {J} _ {0}left ({G} _ {mathit {ij}}} {mathit {ij}}}} ^ {j}}} ^ {i}} ^ {u}} ^ {u}} ^ {u}} {sigma} ^ {u}} {sigma} ^ {u}} {sigma} ^ {u}} {sigma} ^ {u}} {sigma} ^ {u,eta} {sigma} ^ {eta,mathit {ij}}}right) +frac {8} {9} {J} _ {0}left ({G} _ {mathit {ij}}} ^ {mathit {ij}}} ^ {u,mathit {ij}}} ^ {u,etazeta}} {mathit {ij}} ^ {mathit {ij}} ^ {mathit {ij}}} ^ {u,xizeta} {sigma}} ^ {xizeta,mathit {STAB}}}right) STAB STAB STAB