3. The various refinement/deraffination criteria#
We saw in paragraph 2.3 that we can identify the areas to be refined:
or with an error indicator;
or with a field (or its derivative) that is considered to be a measure of the error (temperature, stress, deformation,…). This field can be expressed on knots or on meshes (at Gauss points, at nodes per element or constant per element). If the criterion is expressed at Gauss points or at the nodes per element, then HOMARD chooses the maximum of the value in the mesh in question;
or by boxes.
For error indicators, the user has several choices. This is the purpose of the following paragraphs.
3.1. Error indicators#
3.1.1. What is an error indicator?#
It is important here to redefine what is meant by error indicator (or error estimator) so that the user is well aware of what can and cannot be made to say to this indicator.
A function e is a posteriori error indicator if it verifies the inequality:
\(\parallel u-{u}^{h}\parallel \le e(h,{u}^{h},d)\)
where u is the displacement field of the exact solution, \({u}^{h}\) the displacement field of the finite element solution, \(\parallel \cdot \parallel\) a norm on the displacement fields, h the size of the elements, and d the problem data. Also if \(e(h,d,{u}^{h})\) can be located in the form:
\(e(h,{u}^{h},d)={(\sum _{E}{e}_{E}({u}^{h},d{)}^{2})}^{\frac{1}{2}}\)
then quantities \({e}_{E}({u}^{h},d)\) are called local error indicators.
For the two types of error indicators (therefore two different e functions) present in Code_Aster, the \(\parallel \cdot \parallel\) standard is an energy standard.
3.1.2. What to remember#
These error indicators lead to the estimation of an overall error. The choice of global precision based on the energy standard is often difficult to interpret because there is no direct quantitative link with an error in local mechanical quantities (constraints at a point or on a line, maximum constraints, etc.) and the error in movement. In other words, X% of the error indicator in a mesh does not mean that there is an X% error in stress or deformation.
Used alone, i.e. without the objective of carrying out mesh adaptation, the indicator does not make it possible to quantify the quality of the initial mesh, except possibly to identify problem areas, such as geometric singularities or strong gradients.
3.1.3. Scope of use#
The error indicators used in*Code_Aster* are ex-post indicators and are therefore calculated by post-processing Code_Aster by the CALC_ERREUR command. These are constant fields per element (so one value per finite element). Their field of use is delimited by the following constraints (refer to the reference documents given in paragraph for more details):
the errors taken into account are spatial discretization errors (therefore the size of the elements used); in particular, the temporal discretization errors (or pseudo-temporal in the case of non-linear materials) are outside this perimeter;
physical phenomena are limited to mechanics (linear or non-linear) and thermal (idem. );
in mechanics as in thermics, the behavior can be linear or non-linear (except for the Zhu-Zienkiewicz error estimator in mechanics, which only treats linear behavior), knowing that the theoretical results of the error indicators are obtained in the linear domain (their use in the non-linear domain is therefore not based on theoretical results but on an empirical observation of their interest);
the elements used can be any for the use of error indicators (except for the Zhu-Zienkiewicz error estimator in mechanics, which only treats 2D elements; the estimator \(\mathrm{ZZ2}\) only accepts meshes composed of either triangles or quadrangles); These elements can be linear or quadratic.
3.1.4. Mechanical error estimator#
There are two types of estimator for mechanics:
The Zhu-Zienkiewicz error estimator [bib5] which is calculated directly in the calc_e RREUR operator with the options” ERZ1_ELEM “for the estimator \(\mathrm{ZZ1}\) and” ERZ2_ELEM “for the estimator \(\mathrm{ZZ2}\). This field has three components per element \(E\):
“ERREST”: the total absolute error \({e}_{E}\);
“NUEST”: the total relative error \({e}_{{E}_{\text{REL}}}\mathrm{=}\text{100}\mathrm{\times }\frac{{e}_{E}}{\sqrt{{e}_{E}^{2}+{\mathrm{\parallel }{\sigma }_{h}\mathrm{\parallel }}_{E}^{2}}}\);
“SIGCAL”: the energy norm of the calculated solution \({\parallel {\mathrm{\sigma }}_{h}\parallel }_{E}\).
The residual error estimator [bib6] which is also calculated in the calc_e RREUR operator with the options” ERME_ELEM “for the calculation in the element and” ERME_ELNO “for the calculation at the nodes by elements. This field has eight components per element \(E\):
“ERREST”: the total absolute error \({e}_{E}\);
“TERMRE”: the absolute error of the balance residue;
“TERMSA”: the absolute error of the stress jumps between two meshes;
“TERMNO”: the absolute error of the jumps between the normal stress and the Neumann force;
“NUEST”: the total relative error;
“TERMR2”: the relative error of the balance residue;
“TERMS2”: the relative error of the stress jumps between two meshes;
“TERMN2”: the relative error of the jumps between the normal stress and the Neumann force;
“SIGCAL”: the energy norm of the calculated solution;
For these two types of error estimators, the user also has access to impressions in the result file containing the same information at the global level.
A few tips for making good use of the option **** ERME_ELEM :**
To properly calculate the residual indicator (within the theoretical limits of the formula developed in the elliptical framework with regular border…), it must be calculated across the entire model (CALC_ERREUR (TOUT = “OUI” (value by default))). Note that the model is not necessarily defined over all geometry.
Only loading types are taken into account: PESANTEUR, ROTATION,, FORCE_INTERNE,,, PRES_REP, FORCE_FACE, FORCE_ARETE. Only the last three can be variable. For other types of loads, execution stops in a fatal error.
It is advisable to use finite elements of order 2 in the case of volume forces, otherwise this term is very poorly calculated since DIV (SIGMA) is almost zero!
To take into account the error relating to a null CL of the Neumann type, it must be imposed as a function via a AFFE_CHAR_MECA_F. Via a constant, it will not be taken into account.
In 2D, it takes into account errors on (and between) isoparametric elements SEG2 /3, TRIA3 /6/7, QUAD4 /8/9. In 3D, same with FACE3 /4/6/8/9, TETRA4 /10, PENTA6 /13/15 and HEXA8 /20/27. Only the error on the structural elements (shell, plate, beam…) and the PYRAM are not taken into account.
On the other hand, care must be taken not to intersperse segments between two quadrangles or two triangles (resp. quadrangle or triangle between two hexahedra), otherwise one cannot calculate the jump term relating to this neighborhood.
3.1.5. Thermal error estimator#
There is only one thermal error indicator [bib7] which is of the residual type. It is calculated in operator CALC_ERREUR, by the options” ERTH_ELEM “for element-based calculation and” ERTH_ELNO “for element-based calculation at the nodes. For the “ERTH_ELEM” option, it is mandatory to calculate the flow per element at the “FLUX_ELNO” nodes.
The estimator provides the following components:
ERTABS: total absolute error
ERTREL: total relative error
TERMNO: normalization term for the total relative error
TERMVO: absolute error of the equilibrium residue (volume term)
TERMV2: relative error of the equilibrium residue (volume term)
TERMV1: normalization term for volume error
TERMSA: absolute error of the jumps in the heat flow between two meshes
TERMS2: relative error of jumps in heat flow between two meshes
TERMS1: normalization term for the relative jump error
TERMFL: absolute error on imposed flow type boundary conditions
TERMF2: relative error on boundary conditions such as imposed flow
TERMF1: normalization term on the relative error on the flow
TERMEC: absolute error on exchange-type boundary conditions
TERME2: relative error on exchange-type boundary conditions
TERME1: normalization term on the relative error on exchanges
A few tips for proper use:
specific overload rules concerning loads (alarm generation <A>in case of non-compliance);
calculation on the entire mesh associated with the model (error generation <F>in case of non-compliance) between two contiguous or non-contiguous time steps (generation of alarm <A>in case of non-compliance);
<A>* all 2D-plan/AXI and 3D elements are processed (except PYRAM: alarm generation);
all limit conditions except ECHANGE_PAROI, FLUX_NL and RAYO are taken into account (alarm generation <A>when using ECHANGE_PAROI, FLUX_NL or RAYO);
<A><F>* the mesh tolerates « sketches » but needs to be « cleaned up » a bit (no SEG/FACE interspersed in surfaces/volumes, symmetrization problem, double points: alarm or error generation).
3.2. Interest quantity error estimators#
3.2.1. Warning#
Unlike the classical error indicators above, so-called quantity of interest error indicators allow access to a local error on physical quantities (constraint, displacement), and therefore in the unit of this quantity.
Their use is more delicate and probably requires a quick reading of the documentation [bib9] and a more careful reading of the recommendations given for the use of CALC_ERREUR.
Moreover, as their development is very recent, there is little feedback on their use. The following is therefore entirely based on the reference documentation and a report [bib10].
3.2.2. What do these indicators calculate?#
These error indicators make it possible to estimate a physical quantity in a sub-region of the mesh, for example:
An average of displacement (component \(x\) for example) in this subregion, mathematically expressed by equation \(Q(u)=\frac{1}{\mid \Omega \mid }\underset{\Omega }{\int }{u}_{x}d\omega\).
An average of constraints (for example component \(\mathit{xx}\)) \(Q(u)=\frac{1}{\mid \Omega \mid }\underset{\Omega }{\int }{\sigma }_{\text{xx}}d\omega\).
Once you have chosen the quantity of interest, the indicator estimates the quantity \({\mathrm{\varepsilon }}^{Q}=Q(u)-Q({u}^{h})\).
This is where the difficulty arises for the user because this quantity of interest must be expressed in terms of loading known in AFFE_CHER_MECA.
For example, for the average quantity of interest of the move, it is equivalent to a loading volume:
\(Q(u)\mathrm{=}\frac{1}{\mathrm{\mid }\Omega \mathrm{\mid }}\underset{\Omega }{\mathrm{\int }}{u}_{x}d\Omega \mathrm{=}\frac{1}{\mathrm{\mid }\Omega \mathrm{\mid }}\underset{\Omega }{\mathrm{\int }}fud\Omega\) with \(\overrightarrow{f}=(\mathrm{1,0}\mathrm{,0})\)
By choosing vector \(\overrightarrow{f}\) correctly, we can access each component of the displacement and possibly linear combinations of the components. This load can be accessed by the keyword “FORCE_INTERNE” in “AFFE_CHAR_MECA”.
In the same way, for the average quantity of interest of the stress, it is equivalent to an initial deformation loading:
\(Q(u)\mathrm{=}\frac{1}{\mathrm{\mid }\Omega \mathrm{\mid }}\underset{\Omega }{\mathrm{\int }}{\sigma }_{\text{xx}}d\Omega \mathrm{=}\frac{1}{\mathrm{\mid }\Omega \mathrm{\mid }}\underset{\Omega }{\mathrm{\int }}A\mathrm{:}\sigma (u)d\Omega \mathrm{=}\frac{1}{\mathrm{\mid }\Omega \mathrm{\mid }}\underset{\Omega }{\mathrm{\int }}A\mathrm{:}K\mathrm{:}\varepsilon (u)d\Omega\)
The choice of the component of the stress tensor is controlled by the choice of the components of the \(A\) tensor. This load can be accessed in by the keyword “PRE_EPSI” in “AFFE_CHAR_MECA”.
These two quantities are quite simple to implement because they are represented by a loading existing in Code_Aster. For more complicated quantities (stress tensor functions for example) and possibly non-linear quantities, it is very difficult to express them according to known loads. This is the main difficulty to overcome in order to program new quantities of interest.
3.2.3. Scope of use#
These quantities of interest are calculated using the classical error indicators presented in § 3.1 and only for mechanics. Their field of use is therefore identical to that of the mechanical indicator chosen to calculate them.
We recall, again, that the quantity of interest chosen depends on the degree of development of the loads in AFFE_CHAR_MECA.
3.2.4. The different options in CALC_ERREUR#
The various options available in CALC_ERREUR are:
“QIZ1_ELEM” (respectively “QIZ2_ELEM”): error estimator in quantity of interest based on the Zhu-Zienkiewicz method. You must first calculate the option “ERZ1_ELEM” or “ERZ2_ELEM”.
“QIRE_ELEM” (per element) or “QIRE_ELNO” (per element at the nodes): error estimator in quantity of interest based on residues in mechanics. The option “ERME_ELEM” must be calculated beforehand.
Advice for use:
The field of use of the options” QIZ1_ELEM “and” QIZ2_ELEM “is the same as for the options” ERZ1_ELEM “and” ERZ2_ELEM “and that of the option” QIRE_ELEM “is the same as that of the option” ERME_ELEM “in mechanics.
To date, only two quantities of interest are available: average of a displacement component and average of a component of the stress tensor.
The calculation of the quantity of interest involves the resolution of an elastic calculation (command MECA_STATIQUE) whose loading is on the one hand, the same as that of the main problem on Dirichlet conditions (therefore in imposed displacement), and on the other hand, the one representing the quantity of interest chosen (“FORCE_INTERNE” or “PRE_EPSI”). Below is an example of a command file.
RESOLUTION FROM PROBLEME
CHARG = AFFE_CHAR_MECA (
MODELE =ME,
FACE_IMPO =( _F (GROUP_MA =” GRMA1 “, DY=0.0,),
_F (GROUP_MA =” GRMA2 “, DY=0.91333,),),
FORCE_CONTOUR =( _F (GROUP_MA =” GRMA3 “, FX=1.0,),
_F (GROUP_MA =” GRMA4 “, FX=-2.0,),),);
RESU = MECA_STATIQUE (
MODELE =ME,
CHAM_MATER =CM,
EXCIT =_F (CHARGE = CHARG,),);
RESU = CALC_ERREUR (reuse = RESU,
RESULTAT = RESU,
TOUT_ORDRE =” OUI “,
OPTION =( “SIGM_ELNO”, “ERME_ELEM”,),);
# CALCUL OF THE QUANTITE D INTERET IN MOYENNE OF OF CONTRAINTE SUR THE
GROUPE OF MAILLES GRMAQI
CHARQ = AFFE_CHAR_MECA (
MODELE =ME,
FACE_IMPO =( _F (GROUP_MA =” GRMA1 “, DY=0.0,),
_F (GROUP_MA =” GRMA2 “, DY=0.91333,),),
FORCE_INTERNE =_F (GROUP_MA =” GRMAQI “, FX=1.0,),);
RESUQ = MECA_STATIQUE (MODELE =MO,
CHAM_MATER =CM,
EXCIT =_F (CHARGE = CHARQ,),);
RESUQ = CALC_ERREUR (reuse = RESUQ,
RESULTAT = RESUQ,
TOUT_ORDRE =” OUI “,
OPTION =( “SIGM_ELNO”, “ERME_ELEM”,),);
CALCUL OF ESTIMATEUR OF ERREUR IN QUANTITE OF INTERET
RESU = CALC_ERREUR (reuse = RESU,
RESULTAT = RESU,
OPTION =( “QIRE_ELEM”),
RESU_DUAL = RESUQ,);
3.3. Singularity detection#
Another possible field (only in mechanics) that can be used as a criterion to refine the mesh is the use of the option “SING_ELEM” (constant field per element) or “SING_ELNO” (fields by elements at the nodes) in CALC_ERREUR. This field is calculated from error indicators and it aims to improve the treatment of singularities in mesh adaptation strategies. The reader will find in [bib8] a theoretical description of the methods used.
In practice, the error indicators are high in the singular areas so that quickly only the singular areas are refined. This therefore masks the other sensitive areas (areas with a high gradient) that one would like to refine. This constant field per element has two components:
“DEGRE” which indicates the level of singularity of a mesh. In practice, this component is equal to the degree of interpolation of the finite elements chosen if the finite element is not connected to any singularity and is equal to the order of the singularity if the finite element is connected to a node considered by the method to be singular (for example, for an element close to the tip of a crack, this value is equal to 0.5). This value cannot be greater than the degree of interpolation and therefore 1 or 2.
“RAPPORT” which corresponds to the map for modifying the size of the finite elements in case of remeshing for a given global error. This component is equal to the ratio between the current size and the new size of the finite element.
“TAILLE” which corresponds to the finite element size map in case of remeshing for a given global error.
Advice for use:
The calculation of this option requires, first, the calculation of an error indicator (the absolute component is used and is hard-coded in Code_Aster) and the total deformation energy. If one of these options is not calculated, an alarm message is sent and the option “SING_ELEM” is not calculated.
For the total deformation energy, “ETOT_ELEM” with STAT_NON_LINE or “EPOT_ELEM” with MECA_STATIQUE is used.
The user must also enter the keyword “PREC_ERR” (a fatal message is sent in case of absence) which makes it possible to calculate the desired precision on the global error to determine the size change map. The value of “PREC_ERR” is strictly between 0 and 1 and a fatal message is sent if this condition is not met.
The scope of use is the same (but smaller) as that of the error indicator chosen, namely:
For the residual indicator: finite elements of continuous media in 2D (triangles and quadrangles) or 3D (only tetrahedra) for elastoplastic behavior,
For the Zhu-Zienkiewicz indicator: finite elements of 2D continuous media (triangles and quadrangles) for elastic behavior.
Rigorously, the calculation of the order of the singularity is obtained from the theoretical energy at the crack point, an equation valid only in elasticity. The use of this option in elastoplasticity should therefore be handled with caution.