3. Construction of the recalculated stress field#

3.1. 1987 version#

Since solution \({u}_{h}\) resulting from equation [éq 2.1-1] is \({C}_{0}\) out of \(\Omega\) (due to the choice of functions of the form \({C}_{0}\)), it follows that \({\sigma }_{h}\) calculated by [éq 2.1-2] is discontinuous at the interfaces of the elements.

To obtain acceptable results on nodal stresses, a node average or a projection method is generally used. It is the latter method that is adopted here.

We assume that \({\sigma }^{\text{*}}\) is interpolated by the same form functions as \({u}_{h}\), i.e.:

\({\sigma }^{\text{*}}=N{\stackrel{ˉ}{\sigma }}^{\text{*}}\) eq 3.1-1

and we perform a global least-squares smoothing of \({\sigma }_{h}\), which amounts to minimizing the functional \(J(\tau )=\underset{\Omega }{\int }{}^{t}\text{}(\tau -{\sigma }_{h})(\tau -{\sigma }_{h})d\Omega\) in the space generated by \(N\).

By derivation, \({s}^{}\) should check \(\underset{\Omega }{\int }{}^{t}\text{}N({\sigma }^{\text{*}}-{\sigma }_{h})d\Omega =0\)

using the equation [éq 3.1-1], we get the linear system:

\(M\left\{{\stackrel{ˉ}{\sigma }}^{\text{*}}\right\}=\left\{b\right\}\)

with \(M=\underset{\Omega }{\int }{}^{t}\text{}\mathrm{NN}d\Omega\) and \(\left\{b\right\}=\underset{\Omega }{\int }{}^{t}\text{}N{\sigma }_{h}d\Omega\)

This global system must be solved on each of the components of the stress tensor. Matrix \(M\) is calculated and inverted only once.

3.2. Version 1992#

The constraint of field \({\sigma }^{\text{*}}\) differs from the 1987 version in the following way:

We assume \({\sigma }^{\text{*}}\) polynomial of the same degree as the displacements on all the elements having an internal vertex node \(S\) in common.

We note \({S}_{K}\mathrm{=}\underset{S\mathrm{\in }K}{\mathrm{\cup }}K\) this set called patch.

For each component of \({\sigma }^{\text{*}}\), we write:

\({\sigma }^{\text{*}}{\mid }_{{S}_{K}}={\text{Pa}}_{s}\) eq 3.2-1

where

\(P\) contains the appropriate polynomial terms

\({a}_{s}\) the unknown coefficients of the corresponding monomials

Example:
_images/Object_74.svg

\(\mathrm{P1}P=\left[\mathrm{1,}x,y\right]\) \({a}_{s}{}^{t}=\left[{a}_{1},{a}_{2},{a}_{3}\right]\)

\(\mathrm{Q1}P=\left[\mathrm{1,}x,y,\mathrm{xy}\right]\) \({a}_{s}={}^{t}\text{}\left[{a}_{1},{a}_{2},{a}_{3},{a}_{4}\right]\)

The determination of the coefficients of polynomial \({a}_{s}\) is done by minimizing the functional:

\(\begin{array}{cc}F(a)& \text{=}\sum _{i=1}^{N}{({\sigma }_{h}({x}_{i},{y}_{i})-{\sigma }^{\text{*}}{\mid }_{{S}_{K}}({x}_{i},{y}_{i}))}^{2}\\ & =\sum _{i=1}^{N}{({\sigma }_{h}({x}_{i},{y}_{i})-P({x}_{i},{y}_{i}){a}_{s})}^{2}\end{array}\)

(discrete local smoothing of \({\sigma }_{h}\) by least squares)

where

\(({x}_{i},{y}_{i})\) are the coordinates of the points GAUSS on \({S}_{K}\).

\(N\) is the total number of GAUSS points on all items in the patch

Solution \({a}_{s}\) verifies:

\(\sum _{i=1}^{N}{}^{t}\text{}P({x}_{i},{y}_{i})P({x}_{i},{y}_{i}){a}_{s}=\sum _{i=1}^{N}{}^{t}\text{}P({x}_{i},{y}_{i}){\sigma }_{h}({x}_{i},{y}_{i})\)

Hence \({a}_{s}={A}^{\text{-}1}b\) with \(A=\sum _{i=1}^{N}{}^{t}\text{}P({x}_{i},{y}_{i})P({x}_{i},{y}_{i})\)

\(A\) can be very poorly conditioned (especially on high-grade elements) and therefore impossible to reverse in this form. To remedy this problem, the authors [bib4] proposed normalizing the coordinates on each patch, which is the same as changing variables:

\(\begin{array}{c}\stackrel{ˉ}{x}\mathrm{=}\mathrm{-}1+2\frac{x\mathrm{-}{x}_{\text{min}}}{{x}_{\text{max}}\mathrm{-}{x}_{\text{min}}}\\ \stackrel{ˉ}{y}\mathrm{=}\mathrm{-}1+2\frac{y\mathrm{-}{y}_{\text{min}}}{{y}_{\text{max}}\mathrm{-}{y}_{\text{min}}}\end{array}\)

where \({x}_{\text{min}},{x}_{\text{max}},{y}_{\text{min}},{y}_{\text{max}}\) represent the minimum and maximum values of \(x\) and \(y\) on the patch.

This method significantly improves the packaging of \(A\) and completely eliminates the previous problem.

Once \({\mathrm{a}}_{\mathrm{s}}\) is determined, the nodal values are derived from equation [éq 3.2-1] only on the nodes internal to the patch, except in the case of patches with edge nodes.

Internal patches:

_images/100018F8000069D500003F75204B9781D82C043E.svg
_images/100000AE000000EE000000EE1CF4CBF523C49E33.svg

points in GAUSS where the \({\sigma }_{h}\) constraints are calculated according to the equation []

_images/100000A4000000D4000000D4640DCCB272317F12.svg

\({\sigma }^{\text{*}}\) \({\sigma }^{\text{*}}\) compute nodes

_images/1000010E0000017200000172D03BA2D27C22F56C.svg

internal vertex defining the patch

Edge patches:

_images/1000062C00000C9C00000CEB9CB48DC6D99E3A20.svg

The nodal values at the middle nodes belonging to 2 patches are averaged, the same for the internal nodes in the case of QUAD9.

Note:

In the case of finite elements of different types, the choice of \(P\) in equation [éq 3.2-1] is difficult (validity problems of \({a}_{s}\) if the space is too rich, loss of super‑convergence if it is not rich enough). A more thorough study seems indispensable.

This is why the estimator \(\mathit{ZZ2}\) is currently limited to meshes with only one type of element. This restriction does not exist for \(\mathit{ZZ1}\) .

The authors showed numerically [bib3] that with this choice of \({\sigma }^{\text{*}}\), their estimator was asymptotically accurate for elastic materials whose characteristics are independent of the domain and for all types of elements and that the convergence rates with \(h\) of \(\parallel {e}^{\text{*}}\parallel\) were improved compared to the previous version (especially for degree 2 elements): see test case SSLV110 Manual of Validation), resulting in a better estimation of the error.

An illustration of these convergence rates can be found in reference [bib 5].