3. Stabilization of QUAD4 at a Gauss point#
The stabilization of QUAS4 takes place in two steps:
Enrich the discretized gradient operator :math:`{B}_{c}` and thus make it possible to calculate the deformation energy related to hourglass modes of travel (*hourglass modes);
Interpolate a deformation/stress field to account for the deformation/stresses on the entire element while integrating the QUAD4 in the center (*assumed strain stabilization).
3.1. Variational principle of the problem#
This one is taken from the weak form of the Hu-Washizu variational principle:
\(\delta \pi (\upsilon ,\stackrel{ˉ}{\dot{\varepsilon }},\stackrel{ˉ}{\sigma })\mathrm{=}{\mathrm{\int }}_{{\Omega }_{e}}{\delta \stackrel{ˉ}{\dot{\varepsilon }}}^{T}\text{.}\sigma d\Omega +\delta {\mathrm{\int }}_{{\Omega }_{e}}{\stackrel{ˉ}{\sigma }}^{T}\text{.}(\mathrm{\nabla }{}_{s}\text{}\upsilon \mathrm{-}\stackrel{ˉ}{\dot{\varepsilon }})d\Omega \mathrm{-}\delta {\dot{d}}^{T}\text{.}{f}^{\text{ext}}\mathrm{=}0\) eq 3.1-1
With
\(\stackrel{ˉ}{\dot{\varepsilon }}(x,t)\mathrm{=}B(x)\text{.}\dot{d}(t)\)
« Assumed strain » stabilization is based on the fact that the postulated stress is chosen orthogonal to the difference between the symmetric part of the speed gradient and the deformation rate. Therefore, it allows us to write:
\(\delta {\dot{d}}^{T}\text{.}{\int }_{{\Omega }_{e}}{B}^{T}\text{.}\sigma d\Omega -\delta {\dot{d}}^{T}\text{.}{f}^{\text{ext}}=0\)
And so that:
\({f}^{\text{int}}\mathrm{=}{\mathrm{\int }}_{{\Omega }_{e}}{B}^{T}\text{.}\sigma (\stackrel{ˉ}{\dot{\varepsilon }})d\Omega\)
3.2. Enrichment of the discretized gradient operator#
Enriching the discretized gradient operator \({B}_{c}\) consists in making it a new \(B\) operator by adding a third component \(\gamma\) which, like \({b}_{x}\) and \({b}_{y}\), is a vector of \({\mathrm{\Re }}^{4}\). However, since the initial operator \({B}_{c}\) correctly calculates the gradient of the linear displacement fields, the new \(\gamma\) component must be orthogonal to the linear displacement fields. The steps for calculating this rich operator are detailed in paragraph [§1.2] of the attached document entitled: « Bibliographical Report ».
Note: In this report, the new operator \(B\) relates the strain rate tensor \(\dot{\varepsilon }\) and the nodal velocity vector \({\dot{d}}_{i}\). In the rest of the report, this formulation allows us to reason in terms of displacement increments, and therefore to deal with incremental problems, carried out over several time steps. For elastic calculations (solutions obtained in a time step) we formulate the problem based on nodal displacements: the operator \(B\) then links deformation tensor \(\varepsilon\) and nodal displacements \({u}_{i}\).
The new \(B\) operator is based on the QUAD4 expression for the field of movement written by Belytschko and Bachrach [bib2]:
\({u}_{i}\mathrm{=}({\Delta }^{T}+x\text{.}{b}_{x}^{T}+y\text{.}{b}_{y}^{T}+{h}_{\gamma }^{T}){u}_{i}\) eq 3.2-1
And it is written:
\(B\mathrm{=}{B}_{c}+{B}_{n}\mathrm{=}\left[\begin{array}{cc}{b}_{{\text{x}}^{t}}+{h}_{\text{,x}}{\text{γ}}^{t}& 0\\ 0& {b}_{{\text{y}}^{t}}+{h}_{\text{,y}}{\text{γ}}^{t}\\ {b}_{{\text{y}}^{t}}+{h}_{\text{,y}}{\text{γ}}^{t}& {b}_{{\text{x}}^{t}}+{h}_{\text{,x}}{\text{γ}}^{t}\end{array}\right]\) eq 3.2-2
With:
\(\begin{array}{c}\Delta \mathrm{=}1\mathrm{/}4\left[t\mathrm{-}({t}^{T}\text{.}x){b}_{x}\mathrm{-}({t}^{T}\text{.}y){b}_{y}\right],\gamma \mathrm{=}1\mathrm{/}4\left[h\mathrm{-}({h}^{T}\text{.}x){b}_{x}\mathrm{-}({h}^{T}\text{.}y){b}_{y}\right]\\ {b}_{x}^{t}\mathrm{=}\frac{1}{\mathrm{2A}}\left[({y}_{2}\mathrm{-}{y}_{4}),({y}_{3}\mathrm{-}{y}_{1}),({y}_{4}\mathrm{-}{y}_{2}),({y}_{1}\mathrm{-}{y}_{3})\right]\\ {b}_{y}^{t}\mathrm{=}\frac{1}{\mathrm{2A}}\left[({x}_{4}\mathrm{-}{x}_{2}),({x}_{1}\mathrm{-}{x}_{3}),({x}_{2}\mathrm{-}{x}_{4}),({x}_{3}\mathrm{-}{x}_{1})\right]\\ h\mathrm{=}\xi \eta \end{array}\)
\(\xi ,\eta\) are the reference coordinates. \({u}_{i}\) are the nodal displacement vectors and \(h\) are the values taken by the \(h\) function at the four nodes.
\({d}^{T}\mathrm{=}\left[{u}_{x},{u}_{y}\right]\)
Note that \(\gamma\) can be expressed directly in terms of nodal coordinates:
\(\gamma \mathrm{=}\frac{1}{4}\left[\begin{array}{c}{x}_{2}({y}_{3}\mathrm{-}{y}_{4})+{x}_{3}({y}_{4}\mathrm{-}{y}_{2})+{x}_{4}({y}_{2}\mathrm{-}{y}_{3})\\ {x}_{3}({y}_{1}\mathrm{-}{y}_{4})+{x}_{4}({y}_{3}\mathrm{-}{y}_{1})+{x}_{1}({y}_{4}\mathrm{-}{y}_{3})\\ {x}_{4}({y}_{1}\mathrm{-}{y}_{2})+{x}_{1}({y}_{2}\mathrm{-}{y}_{4})+{x}_{2}({y}_{4}\mathrm{-}{y}_{1})\\ {x}_{1}({y}_{3}\mathrm{-}{y}_{2})+{x}_{2}({y}_{1}\mathrm{-}{y}_{3})+{x}_{3}({y}_{2}\mathrm{-}{y}_{1})\end{array}\right]\) eq 3.2-3
The [éq 3.2-2] form of the \(B\) operator is equivalent to that of the QUAD4 operator. However, this particular expression of the operator makes it possible to differentiate between the terms to be integrated into the center and the stabilization terms. It is only by intervening on the value of these stabilization terms that we will be able to improve the performance of the element.
3.3. Interpolation of the deformation field:#
The general form of an « assumed strain » field within a QUAD4 is given by Belytschko and Bindeman [bib3]. It has the following form:
\({\epsilon }_{\text{assumed}\text{strain}}=\left[\begin{array}{c}{\epsilon }_{x(\text{centre})}+{q}_{x}{e}_{1}{h}_{,x}+{q}_{y}{e}_{2}{h}_{,y}\\ {\epsilon }_{y(\text{centre})}+{q}_{x}{e}_{2}{h}_{,x}+{q}_{y}{e}_{1}{h}_{,y}\\ {\mathrm{2\epsilon }}_{\text{xy}(\text{centre})}+{q}_{x}{e}_{3}{h}_{,y}+{q}_{y}{e}_{3}{h}_{,x}\end{array}\right]=\left[\begin{array}{c}{\epsilon }_{x(\text{centre})}+{\epsilon }_{x(\text{stab})}\\ {\epsilon }_{y(\text{centre})}+{\epsilon }_{y(\text{stab})}\\ {\mathrm{2\epsilon }}_{\text{xy}(\text{centre})}+{\mathrm{2\epsilon }}_{\text{xy}(\text{stab})}\end{array}\right]\) eq 3.3-1
With:
\(\begin{array}{c}{q}_{x}\mathrm{=}\gamma \text{.}{u}_{x}\\ {q}_{y}\mathrm{=}\gamma \text{.}{u}_{y}\end{array}\)
and \({e}_{1},{e}_{2},{e}_{3}\) which vary according to the physical considerations of each author [Tableau 3.3-1].
Each triplet of values characterizes an element and gives rise to a particular interpolation of the deformation:
Element |
\(\mathit{e1}\) |
|
|
QUAD4 |
1 |
0 |
1 |
ASMD |
1/2 |
-1/2 |
1 |
ASBQI |
1 |
|
0 |
ASOI |
1 |
-1 |
0 |
ASOI (1/2) |
1/2 |
-1/2 |
0 |
Table 3.3-1
Thus we can deduce the expression for the discretized gradient operator arising from the assumed deformation field (\({\varepsilon }_{(\mathit{stab})}\)) of the element. We note this new operator \({B}_{n}\).
\(\begin{array}{}\text{QUAD4}:{B}_{n}=\left[\begin{array}{cc}{h}_{,x}{\gamma }^{t}& 0\\ 0& {h}_{,y}{\gamma }^{t}\\ {h}_{,y}{\gamma }^{t}& {h}_{,x}{\gamma }^{t}\end{array}\right]\\ \text{ASBQI}:{B}_{n}=\left[\begin{array}{cc}{h}_{\text{,x}}{\text{γ}}^{t}& -\stackrel{ˉ}{\nu }{\text{h}}_{\text{,y}}{\text{γ}}^{t}\\ -\stackrel{ˉ}{\nu }{\text{h}}_{\text{,x}}{\text{γ}}^{t}& {h}_{\text{,y}}{\text{γ}}^{t}\\ 0& 0\end{array}\right]\\ \text{ASOI}(\text{1/2}):{B}_{n}=\left[\begin{array}{cc}\frac{1}{2}{h}_{,x}{\gamma }^{t}& -\frac{1}{2}{h}_{,y}{\gamma }^{t}\\ -\frac{1}{2}{h}_{,x}{\gamma }^{t}& \frac{1}{2}{h}_{,y}{\gamma }^{t}\\ 0& 0\end{array}\right]\end{array}\) eq 3.3-2
This writing allows us to easily test the various elements.
The stiffness matrix is then written as follows:
\({K}_{e}={K}_{c}+{K}^{\text{stab}}\) eq 3.3-3
With:
\({K}_{c}\mathrm{=}\underset{{\Omega }_{e}}{\mathrm{\int }}{B}_{{c}^{T}}C{B}_{c}d\Omega \mathrm{=}A\text{.}{B}_{c}\text{.}C\text{.}{B}_{c}\) eq 3.3-4
And
\(\begin{array}{}{K}_{\text{stab}}=\underset{{O}_{e}}{\int }{B}_{{c}^{T}}C{B}_{n}\text{d}O+\underset{{O}_{e}}{\int }{B}_{{n}^{T}}{\text{C B}}_{c}\text{d}O+\underset{{O}_{e}}{\int }{B}_{{n}^{T}}{\text{C B}}_{n}\text{d}O\\ =\sum _{i=1}^{4}\text{JAC}(i)\text{.}({B}_{c}^{T}\text{.}C\text{.}{B}_{n}(i)+{B}_{n}^{T}(i)\text{.}C\text{.}{B}_{c}+{B}_{n}^{T}(i)\text{.}C\text{.}{B}_{n}(i))\end{array}\) eq 3.3-5
Finally we calculate the internal forces in the following way:
\(\begin{array}{}{F}_{\text{int}}=B\text{.}(C\text{.}({\epsilon }_{(\text{centre})}+{\epsilon }_{(\text{stab})}))\\ =\sum _{i=1}^{4}\text{JAC}(i)\text{.}(\text{}({B}_{c}+{B}_{n}(\text{i}))\text{}\text{.}\text{}(C\text{.}({\epsilon }_{c}+\text{}{\epsilon }_{\text{stab}}(\text{i})))\text{})\end{array}\) eq 3.3-6
Notes:
Although the calculation of \({K}_{\text{stab}}\) requires a sum on the four Gauss points, the integration of the law of behavior, which determines the value of the terms of \(C\), takes place at the center. Moreover, the equations [éq 3.3-4] and [éq 3.3-6] show us that the calculations remain relatively bulky. A solution to this problem (not yet implemented in Code_Aster) consists in performing the calculations by placing yourself in a rotating coordinate system with the element (cf. [§3.4] of the bibliographic report). This has the advantage of:
to remove the calculation of cross terms: \({B}_{c}\text{.}C\text{.}{B}_{n}\text{et}{B}_{n}\text{.}C\text{.}{B}_{c}\);
better treatment of transverse shear blockage;
a writing of the law of behavior better adapted to problems including geometric nonlinearities.