Figure 1.1-1: Géométrie du domaine et position de l’interface
v3.04.321 SSLV321 - Validation du pré-conditionnement XFEM pour une interface rasante#
Résumé:
Ce test valide le pré-conditionneur dédié aux éléments XFEM quadratiques. Le problème correspond à celui d’un bloc 3D coupé par une interface XFEM. On applique une compression multiaxiale sur toute la surface du bloc, y compris sur l’interface XFEM. La solution analytique attendue, est un déplacement \({P}_{1}\) à l’intérieur du domaine. Cette solution est captée naturellement par les éléments \({P}_{1}\) et \({P}_{2}\) de Code_Aster.
Lorsque l’interface passe près des nœuds du maillage, le conditionnement de la matrice de rigidité augmente [1,2]. Les stratégies usuelles, utilisées dans le Code_Aster (réajustement de l’interface aux nœuds et élimination de ddls Heaviside), s’appuyant sur des critères «estimatifs», ne sont pas toujours adaptées pour maîtriser le conditionnement et la précision des résultats [2]. Ces inconvénients sont d’autant plus pénalisant, avec des éléments quadratiques.
La nouvelle stratégie (construction d’un pré-conditionneur diagonal par blocs) permet de palier tous ces inconvénients [3].
Ce test comporte 5 modélisations:
\(A\) : une modélisation HEXA8 \(\mathrm{3D}\) avec une interface XFEM d’équation \(x+y+z+0.01=0\)
\(B\) : une modélisation HEXA20 \(\mathrm{3D}\) avec une interface XFEM d’équation \(x+y+z+0.1=0\)
\(C\) : une modélisation HEXA20 \(\mathrm{3D}\) avec une interface XFEM d’équation \(x+y+z+0.01=0\)
\(D\) : une modélisation HEXA8 \(\mathrm{3D}\) avec une interface XFEM d’équation \(x+y+z+0.011=0\)
\(E\) : une modélisation HEXA8 \(\mathrm{3D}\) avec une interface XFEM d’équation \(x+y+z+0.011=0\)
Solution de référence#
Méthode de calcul utilisée pour la solution de référence#
La solution analytique en déplacement, dans la base canonique \((\overrightarrow{{e}_{x}},\overrightarrow{{e}_{y}},\overrightarrow{{e}_{z}})\) , vaut:
\(U(x,y,z)=k\left\lbrace \begin{array}{c}x\\ y\\ z\end{array}\right\rbrace \pm 2\) où, \(k=-\frac{p(1-2\nu )}{E}={4.10}^{-3}\) .
Par la suite, on vérifie que cette expression du champ déplacement, est l’unique solution du problème d’équilibre de part et d’autre de l’interface XFEM, c’est-à-dire, que la solution proposée § 2.1 , respecte les conditions limites et les hypothèses du problème HPP.
Par construction, la solution en déplacement vérifie les conditions de Dirichlet du § 1.3 .
Le tenseur des déformations vaut: \(\underline{\underline{\epsilon}}=\frac{1}{2}(\nabla U+\nabla {U}^{T})=\left[\begin{array}{ccc}k& 0& 0\\ 0& k& 0\\ 0& 0& k\end{array}\right]=kI\) .
En prenant \(k=0,4\text{\%}\) , on respecte l’hypothèse des petites perturbations.
La loi de Hooke, pour un matériau isotrope, valant:
\(\underline{\underline{\sigma}}=\frac{E}{1+\nu }(\underline{\underline{\epsilon}}+\frac{\nu}{1-2\nu }\mathit{Tr}(\underline{\underline{\epsilon}})I)=k\frac{E}{1-2\nu }I\)
or, par construction, on a: \(k=-\frac{p(1-2\nu )}{E}\)
avec, \(\underline{\underline{\sigma}}=-pI\) .
Enfin, les conditions limites de Neumann à imposer sur les bords du domaine \(-\underline{\underline{\sigma}}.\overrightarrow{n}=p\) , correspondent bien au chargement du § 1.3 .
Comme, le tenseur des contraintes est uniforme, il respecte l’équation d’équilibre statique, \(div\underline{\underline{\sigma}}=0\) .
Résultats de référence#
On teste l’erreur maximale, sur les déplacements en quelques points situés sur l’interface \(x+y+z+\mathit{cte}=0\) . La valeur du déplacement en ces points correspond à l’interpolation du champ de déplacement calculé par Aster. Le déplacement calculé par Aster, est alors comparé à la valeur analytique du déplacement donnée au § 2.1 :
\(U(x,y,z)={4.10}^{-3}\left\lbrace \begin{array}{c}x\\ y\\ z\end{array}\right\rbrace +\mathit{sign}(x+y+z+\mathit{cte})\times 2\)
Modélisation A#
Caractéristiques de la modélisation#
C’est une modélisation \(\mathrm{3D}\) XFEM avec des éléments linéaires.
Cette modélisation démontre le bon fonctionnement des critères «estimatifs» avec des éléments linéaires. L’interface est positionnée telle que le seuil de réajustement de l’interface à 1% est atteint: l’interface est alors modifiée et réajustée aux nœuds les plus proches. Cette modification permet de maîtriser le conditionnement. En revanche, le réajustement de l’interface dégrade sensiblement la qualité des résultats pour les éléments linéaires (pour les éléments quadratiques, on supprime ce réajustement, pour le remplacer par un pré-conditionnement algébrique).
L’équation de l’interface XFEM est : \(x+y+z+0.01=0\) .
Caractéristiques du maillage#
Nombre de nœuds: 125
Maillage: 64 mailles de type HEXA8
Figure 3.2-1: maillage linéaire
Grandeurs testées et résultats#
On teste l’erreur maximale sur le déplacement, en valeur absolue, le long de l’interface XFEM. On extrait une liste de points situés sur l’interface, ainsi que le champ de déplacement interpolé en ces points.
Comme le champ de déplacement est discontinu, il existe en chaque point de l’interface 2 valeurs analytiques du champ de déplacement \({U}^{+}\left(x,y,z\right)=\lbrace {U}_{1}^{+},{U}_{2}^{+},{U}_{3}^{+}\rbrace\) et \({U}^{-}\left(x,y,z\right)=\lbrace {U}_{1}^{-},{U}_{2}^{-},{U}_{3}^{-}\rbrace\) .On compare ces valeurs analytiques, au champ de déplacement interpolé en chaque point.
En pratique, dans le Code_Aster pour tenir compte de la discontinuité lors de l’interpolation, chaque point sur l’interface est transformé en nœuds doublon (NP* et NM*) auxquels sont associées des valeurs de déplacements \({\mathit{DX}}_{i}(\mathit{NP})\) et \({\mathit{DX}}_{i}(\mathit{NM})\) .
Pour les nœuds «PLUS» (noté NP* par défaut dans le code_Aster) , on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NP})}_{i}=∣{U}_{i}^{+}({x}_{\mathit{NP}},{y}_{\mathit{NP}},{z}_{\mathit{NP}})-{\mathit{DX}}_{i}(\mathit{NP})∣\) ;
pour les nœuds «MOINS» (noté NM* par défaut dans le code_Aster), on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NM})}_{i}=∣{U}_{i}^{-}({x}_{\mathit{NM}},{y}_{\mathit{NM}},{z}_{\mathit{NM}})-{\mathit{DX}}_{i}(\mathit{NM})∣\) .
Identification |
Référence |
Type |
% tolérance |
\(\text{DIFF}{\text{(NP)}}_{X}\) (MAX) |
0.0 |
Analytique |
0.2 |
\(\text{DIFF}{\text{(NP)}}_{Y}\) (MAX) |
0.0 |
Analytique |
0.2 |
\(\text{DIFF}{\text{(NP)}}_{Z}\) (MAX) |
0.0 |
Analytique |
0.2 |
\(\text{DIFF}{\text{(NM)}}_{X}\) (MAX) |
0.0 |
Analytique |
0.2 |
\(\text{DIFF}{\text{(NM)}}_{Y}\) (MAX) |
0.0 |
Analytique |
0.2 |
\(\text{DIFF}{\text{(NM)}}_{Z}\) (MAX) |
0.0 |
Analytique |
0.2 |
Tableau 3.3-1 : résultats synthétiques
Remarque#
Les tolérances analytiques sont beaucoup plus élevées que pour les modélisations quadratiques \(B\) et \(C\) .
Modélisation B#
Caractéristiques de la modélisation#
C’est une modélisation \(\mathrm{3D}\) XFEM avec des éléments quadratiques.
Cette modélisation démontre la nécessité du pré-conditionneur pour le calcul avec des éléments quadratiques. L’interface est positionnée à 10% (de la longueur de l’arête) au-delà du seuil de réajustement de l’interface à 1%, mais le conditionnement de la matrice de rigidité est déjà élevé: sans le pré-conditionneur, le test ne peut pas tourner de manière robuste. Avec le pré-conditionneur, le solveur direct calcule une solution de manière robuste.
L’équation de l’interface XFEM est : \(x+y+z+0.1=0\) .
Caractéristiques du maillage#
Nombre de nœuds: 425
Maillage: 64 mailles de type HEXA20
Figure 4.2-1: maillage quadratique
Grandeurs testées et résultats#
On teste l’erreur maximale sur le déplacement, en valeur absolue, le long de l’interface XFEM. On extrait une liste de points situés sur l’interface, ainsi que le champ de déplacement interpolé en ces points.
Comme le champ de déplacement est discontinu, il existe en chaque point de l’interface 2 valeurs analytiques du champ de déplacement \({U}^{+}\left(x,y,z\right)=\lbrace {U}_{1}^{+},{U}_{2}^{+},{U}_{3}^{+}\rbrace\) et \({U}^{-}\left(x,y,z\right)=\lbrace {U}_{1}^{-},{U}_{2}^{-},{U}_{3}^{-}\rbrace\) .On compare ces valeurs analytiques, au champ de déplacement interpolé en chaque point par ASTER.
En pratique, dans le Code_Aster pour tenir compte de la discontinuité lors de l’interpolation, chaque point sur l’interface est transformé en nœuds doublon (NP* et NM*) auxquels sont associées des valeurs de déplacements \({\mathit{DX}}_{i}(\mathit{NP})\) et \({\mathit{DX}}_{i}(\mathit{NM})\) .
Pour les nœuds «PLUS» (noté NP* par défaut dans le code_Aster) , on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NP})}_{i}=∣{U}_{i}^{+}({x}_{\mathit{NP}},{y}_{\mathit{NP}},{z}_{\mathit{NP}})-{\mathit{DX}}_{i}(\mathit{NP})∣\) ;
pour les nœuds «MOINS» (noté NM* par défaut dans le code_Aster), on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NM})}_{i}=∣{U}_{i}^{-}({x}_{\mathit{NM}},{y}_{\mathit{NM}},{z}_{\mathit{NM}})-{\mathit{DX}}_{i}(\mathit{NM})∣\) .
Identification |
Référence |
Type |
% tolérance |
\(\text{DIFF}{\text{(NP)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NP)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NP)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NM)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-04 |
\(\text{DIFF}{\text{(NM)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-04 |
\(\text{DIFF}{\text{(NM)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-04 |
Tableau 4.3-1 : résultats synthétiques
Remarque#
Les tolérances analytiques sont beaucoup plus faibles que pour la modélisation \(A\) . Cette différence s’explique par le réajustement de l’interface aux nœuds pour les éléments linéaires de la modélisation \(A\) . Si l’interface est positionnée loin du seuil de réajustement à 1%, on constate une amélioration des résultats pour les éléments linéaires.
Modélisation C#
Caractéristiques de la modélisation#
C’est une modélisation \(\mathrm{3D}\) XFEM avec des éléments quadratiques.
L’interface est positionnée exactement comme dans la modélisation \(A\) (à 1% de l’arête). Comme, les critères estimatifs ne sont pas adaptés aux éléments quadratiques, seul le pré-conditionneur contrôle le conditionnement de la matrice de rigidité. Le test démontre que le pré-conditionneur rend d’une part, possible un positionnement de l’interface XFEM, totalement indépendant de la proximité des nœuds du maillage et d’autre part, augmente considérablement la précision des résultats par rapport à la modélisation \(A\) .
L’équation de l’interface XFEM est : \(x+y+z+0.01=0\) , dans la même configuration que la modélisation \(A\) .
Caractéristiques du maillage#
Même maillage que la modélisation \(B\) .
Grandeurs testées et résultats#
On teste l’erreur maximale sur le déplacement, en valeur absolue, le long de l’interface XFEM. On extrait une liste de points situés sur l’interface, ainsi que le champ de déplacement interpolé en ces points.
Comme le champ de déplacement est discontinu, il existe en chaque point de l’interface 2 valeurs analytiques du champ de déplacement \({U}^{+}\left(x,y,z\right)=\lbrace {U}_{1}^{+},{U}_{2}^{+},{U}_{3}^{+}\rbrace\) et \({U}^{-}\left(x,y,z\right)=\lbrace {U}_{1}^{-},{U}_{2}^{-},{U}_{3}^{-}\rbrace\) .On compare ces valeurs analytiques, au champ de déplacement interpolé en chaque point par ASTER.
En pratique, dans le Code_Aster pour tenir compte de la discontinuité lors de l’interpolation, chaque point sur l’interface est transformé en nœuds doublon (NP* et NM*) auxquels sont associées des valeurs de déplacements \({\mathit{DX}}_{i}(\mathit{NP})\) et \({\mathit{DX}}_{i}(\mathit{NM})\) .
Pour les nœuds «PLUS» (noté NP* par défaut dans le code_Aster) , on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NP})}_{i}=∣{U}_{i}^{+}({x}_{\mathit{NP}},{y}_{\mathit{NP}},{z}_{\mathit{NP}})-{\mathit{DX}}_{i}(\mathit{NP})∣\) ;
pour les nœuds «MOINS» (noté NM* par défaut dans le code_Aster), on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NM})}_{i}=∣{U}_{i}^{-}({x}_{\mathit{NM}},{y}_{\mathit{NM}},{z}_{\mathit{NM}})-{\mathit{DX}}_{i}(\mathit{NM})∣\) .
Identification |
Référence |
Type |
% tolérance |
\(\text{DIFF}{\text{(NP)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NP)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NP)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NM)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-03 |
\(\text{DIFF}{\text{(NM)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-03 |
\(\text{DIFF}{\text{(NM)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-03 |
Tableau 5.3-1 : résultats synthétiques
Modélisation D#
Caractéristiques de la modélisation#
C’est une modélisation \(\mathrm{3D}\) XFEM avec des éléments linéaires.
Cette modélisation permet d’effectuer une comparaison directe entre 2 méthodes de pré-conditionnement: l’élimination des ddls Heaviside avec le critère de rigidité (modélisation \(D\) ) et le pré-conditionneur XFEM (cf. modélisation \(E\) ). On se place tout juste avant le seuil de réajustement de l’interface, à 1.1% de la longueur de l’arête. Dans la commande MODI_MODELE_XFEM, on laisse l’option PRETRAITEMENTS=”AUTO” [U4.41.11], qui désactive le pré-conditionneur XFEM pour les éléments linéaires.
L’équation de l’interface XFEM est : \(x+y+z+0.011=0\) . Le pré-conditionneur n’est pas activé dans les PRETRAITEMENTS=”AUTO”[U4.41.11] pour cette modélisation. Le critère d’élimination de ddls Heaviside est donc activé.
Caractéristiques du maillage#
Même maillage que la modélisation \(\text{A}\) .
Grandeurs testées et résultats#
On teste l’erreur maximale sur le déplacement, en valeur absolue, le long de l’interface XFEM. On extrait une liste de points situés sur l’interface, ainsi que le champ de déplacement interpolé en ces points.
Comme le champ de déplacement est discontinu, il existe en chaque point de l’interface 2 valeurs analytiques du champ de déplacement \({U}^{+}\left(x,y,z\right)=\lbrace {U}_{1}^{+},{U}_{2}^{+},{U}_{3}^{+}\rbrace\) et \({U}^{-}\left(x,y,z\right)=\lbrace {U}_{1}^{-},{U}_{2}^{-},{U}_{3}^{-}\rbrace\) .On compare ces valeurs analytiques, au champ de déplacement interpolé en chaque point.
En pratique, dans le Code_Aster pour tenir compte de la discontinuité lors de l’interpolation, chaque point sur l’interface est transformé en nœuds doublon (NP* et NM*) auxquels sont associées des valeurs de déplacements \({\mathit{DX}}_{i}(\mathit{NP})\) et \({\mathit{DX}}_{i}(\mathit{NM})\) .
Pour les nœuds «PLUS» (noté NP* par défaut dans le code_Aster) , on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NP})}_{i}=∣{U}_{i}^{+}({x}_{\mathit{NP}},{y}_{\mathit{NP}},{z}_{\mathit{NP}})-{\mathit{DX}}_{i}(\mathit{NP})∣\) ;
pour les nœuds «MOINS» (noté NM* par défaut dans le code_Aster), on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NM})}_{i}=∣{U}_{i}^{-}({x}_{\mathit{NM}},{y}_{\mathit{NM}},{z}_{\mathit{NM}})-{\mathit{DX}}_{i}(\mathit{NM})∣\) .
Identification |
Référence |
Type |
% tolérance |
\(\text{DIFF}{\text{(NP)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-07 |
\(\text{DIFF}{\text{(NP)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-07 |
\(\text{DIFF}{\text{(NP)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-07 |
\(\text{DIFF}{\text{(NM)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-01 |
\(\text{DIFF}{\text{(NM)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-01 |
\(\text{DIFF}{\text{(NM)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-01 |
Tableau 6.3-1 : résultats synthétiques
Remarque#
On gagne plusieurs ordres de grandeurs sur la précision des résultats, par rapport à la modélisation \(\text{A}\) , sachant que, l’interface a été décalée de \(0.1\text{\%}\) .
Le réajustement de l’interface (avec un seuil à 1%) a donc une influence néfaste sur la qualité des résultats.
Néanmoins, pour contrôler le conditionnement dans la modélisation \(\text{D}\) , on constate l’activation de l’élimination de ddls Heavisides, grâce à un critère d’estimation de rigité [R7.02.12].
Modélisation E#
Caractéristiques de la modélisation#
C’est une modélisation \(\mathrm{3D}\) XFEM avec des éléments linéaires.
Cette modélisation permet d’effectuer une comparaison directe entre 2 méthodes de pré-conditionnement: l’élimination des ddls Heaviside avec le critère de rigidité (modélisation \(D\) ) et le pré-conditionneur XFEM (modélisation \(E\) ). Comme dans la modélisation \(D\) , on se place tout juste avant le seuil de réajustement de l’interface, à 1.1% de la longueur de l’arête. Dans la commande MODI_MODELE_XFEM, on passe l’option PRETRAITEMENTS=”FORCE”[U4.41.11],qui active le pré-conditionneur XFEM. L’utilisation du pré-conditionneur XFEM permet de gagner 3 ordres de grandeurs sur la précision des résultats.
L’équation de l’interface XFEM est : \(x+y+z+0.011=0\) . Le pré-conditionneur est activé par PRETRAITEMENTS=”FORCE”[U4.41.11] pour cette modélisation. Le critère d’élimination de ddls Heaviside n’est donc pas activé.
Caractéristiques du maillage#
Même maillage que la modélisation \(\text{A}\) .
Grandeurs testées et résultats#
On teste l’erreur maximale sur le déplacement, en valeur absolue, le long de l’interface XFEM. On extrait une liste de points situés sur l’interface, ainsi que le champ de déplacement interpolé en ces points.
Comme le champ de déplacement est discontinu, il existe en chaque point de l’interface 2 valeurs analytiques du champ de déplacement \({U}^{+}\left(x,y,z\right)=\lbrace {U}_{1}^{+},{U}_{2}^{+},{U}_{3}^{+}\rbrace\) et \({U}^{-}\left(x,y,z\right)=\lbrace {U}_{1}^{-},{U}_{2}^{-},{U}_{3}^{-}\rbrace\) .On compare ces valeurs analytiques, au champ de déplacement interpolé en chaque point.
En pratique, dans le Code_Aster pour tenir compte de la discontinuité lors de l’interpolation, chaque point sur l’interface est transformé en nœuds doublon (NP* et NM*) auxquels sont associées des valeurs de déplacements \({\mathit{DX}}_{i}(\mathit{NP})\) et \({\mathit{DX}}_{i}(\mathit{NM})\) .
Pour les nœuds «PLUS» (noté NP* par défaut dans le code_Aster) , on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NP})}_{i}=∣{U}_{i}^{+}({x}_{\mathit{NP}},{y}_{\mathit{NP}},{z}_{\mathit{NP}})-{\mathit{DX}}_{i}(\mathit{NP})∣\) ;
pour les nœuds «MOINS» (noté NM* par défaut dans le code_Aster), on calcule alors la table de différence suivante \(\mathit{DIFF}{(\mathit{NM})}_{i}=∣{U}_{i}^{-}({x}_{\mathit{NM}},{y}_{\mathit{NM}},{z}_{\mathit{NM}})-{\mathit{DX}}_{i}(\mathit{NM})∣\) .
Identification |
Référence |
Type |
% tolérance |
\(\text{DIFF}{\text{(NP)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NP)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NP)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-09 |
\(\text{DIFF}{\text{(NM)}}_{X}\) (MAX) |
0.0 |
Analytique |
1.E-03 |
\(\text{DIFF}{\text{(NM)}}_{Y}\) (MAX) |
0.0 |
Analytique |
1.E-03 |
\(\text{DIFF}{\text{(NM)}}_{Z}\) (MAX) |
0.0 |
Analytique |
1.E-03 |
Tableau 7.3-1 : résultats synthétiques
Remarque#
On gagne plusieurs ordres de grandeurs sur la précision des résultats, par rapport à la modélisation \(\text{D}\) , pour la même position de l’interface. Dans la modélisation \(\text{E}\) , le conditionnement est contrôlé par le pré-conditionneur et du même ordre de grandeur que pour la modélation \(\text{D}\) (autour de \({10}^{3}\) ).
Synthèse des résultats#
Ce test permet de valider le fonctionnement du pré-conditionnement des ddls Heavisides pour les éléments quadratiques XFEM.
On constate que cette stratégie de pré-conditionnement donne de meilleurs résultats, que le réajustement de la level-set avec des éléments linéaires, ainsi que l’élimination des ddls par le critère de rigidité.
Pour les éléments quadratiques, le pré-conditionneur est indispensable pour obtenir une solution robuste.