u2.08.08 Utilisation de la Méthode des Solutions Manufacturées pour la vérification logicielle#

Résumé :

Ce document présente l’utilisation de la Méthode des Solutions Manufacturées pour la validation logicielle de Code_Aster . Le principe de cette méthode est détaillé sur un exemple simple de diffusion thermique. Cet exemple nous permet de préciser l’intérêt de cette approche et particulièrement de présenter la démarche de validation qui peut être mise en œuvre. Enfin, l’implantation informatique de la Méthode des Solutions Manufacturées, basée sur l’utilisation du calcul symbolique, est détaillée.

Principe général et utilisation pour la vérification#

La Méthode des Solutions Manufacturées en thermique#

Le principe de cette méthode est très simple; nous allons l’illustrer sur un exemple de thermique linéaire 2D. Avant tout, nous rappelons les équations d’un problème de thermique linéaire. Dans le domaine \(\Omega\) , on recherche le champ de température \(T(x,y)\) tel que:

(109)#\[\begin{split}\lbrace \begin{array}{c}-\lambda \Delta T=s\text{dans}\Omega \\ T={T}_{d}\text{sur}\partial {\Omega}_{d}\\ \lambda \overrightarrow{\nabla}T\cdot \overrightarrow{n}={q}_{n}\text{sur}\partial {\Omega}_{n}\end{array}\end{split}\]

où \({T}_{d}\) et \({q}_{n}\) sont respectivement la température et le flux imposés sur des parties de la frontière.

Les étapes de la Méthode des Solutions Manufacturées sont les suivantes:

  • On choisit un domaine d’étude 2D arbitraire \(\Omega\) comportant les frontières \(\partial {\Omega}_{d}\) avec conditions de Dirichlet et \(\partial {\Omega}_{n}\) avec conditions de Neumann. Par choisir, on entend «choisir la forme du domaine». Ainsi nous choisissons la forme de la figure .

  • On choisit une solution arbitraire, disons

    ../../../../_images/Object_1104.svg

  • En ré-injectant l’expression précédente dans , on déduit les expressions:

  • du terme source;

(110)#\[s(x,y)=-\lambda \Delta T=-\lambda (6x+6y)\]

  • du terme de Dirichlet;

(111)#\[{T}_{d}(x,y)={x}^{3}+{y}^{3}\]

  • du terme de Neumann:

(112)#\[{q}_{d}(x,y)=\lambda \overrightarrow{\nabla}T(x,y)\cdot \overrightarrow{n}\]

Si on veut imposer cette condition sur le bord \(\partial {\Omega}_{n}\) , alors on on connait l’expression de la normale \(\overrightarrow{n}\) à savoir \(\overrightarrow{n}={[1,0]}^{T}\) . On obtient donc:

(113)#\[{q}_{d}(x,y)=\lambda \overrightarrow{\nabla}T(x,y)\cdot \overrightarrow{n}=3\lambda {x}^{2}\]

Arrivé à ce stade, on dispose des informations nécessaires à la modélisation du problème étudié. Il est néanmoins essentiel de faire la remarque suivante: pour obtenir la solution recherchée, il faut bien affecter des conditions aux limites sur tout le bord \(\partial \Omega\) du domaine d’étude. Au choix, suite à une partition de la frontière, cela peut être un mélange de conditions de Dirichlet et de Neumann. L’obligation d’imposer des conditions sur toute la frontière est notamment dû au fait qu’une absence de condition aux limites est implicitement une condition de Neumann nul. Or l’équation montre que l’on n’a pas le choix de la valeur de cette condition. Dans l’exemple présent, nous faisons le choix d’imposer une condition de Neumann sur \(\partial {\Omega}_{n}\) et des conditions de Dirichlet sur \(\partial {\Omega}_{d}\) , \(\partial {\Omega}_{d}^{1}\) et \(\partial {\Omega}_{d}^{2}\) .

../../../../_images/Cadre116.gif

Dessin 1: Domaine d’étude

Une fois, les équations précédentes obtenues, on passe à la phase modélisation. On commence par choisir une discrétisation éléments finis \({V}_{h}\) (linéaire, quadratique). On modélise dans son logiciel d’éléments finis préféré le domaine \({\Omega}_{h}\) auquel on applique les chargements \({s}_{h}(x,y)\) , \({T}_{\mathrm{dh}}(x,y)\) et \({q}_{\mathrm{dh}}(x,y)\) . Après résolution, on obtient la solution approchée \({T}_{h}(x,y)\) .

Techniques de vérification#

Grâce à la MSM, différentes vérifications sont réalisables.

Choisir une solution de référence dans l’espace d’approximation \(T\in {V}_{h}\)#

C’est la vérification la plus simple car on doit alors obtenir \(T(x,y)={T}_{h}(x,y)\) . A noter que, du fait de résolution numérique, on n’aura pas égalité parfaite mais la différence doit être de l’ordre de la précision machine \(\text{~}1.E-15\) . De manière concrète, si on utilise des éléments linéaires, il faudra alors choisir une solution linéaire.

Choisir une solution de référence hors de l’espace d’approximation \(T\notin {V}_{h}\)#

Dans le cas de figure où \(T\notin {V}_{h}\) , on n’a plus l’égalité du paragraphe précédent. Par contre, on peut assez facilement mesurer la vitesse de convergence de \({T}_{h}\) vers \(T\) quand \(h\to 0\) .

On précise que la vitesse de convergence, avec la finesse \(h\) du maillage, de la solution calculée vers la solution analytique, dans une norme donnée, se définit comme le plus grand réel \(\alpha >0\) tel que \(\parallel T(x,y)-{T}_{h}(x,y)\parallel <C\ast {h}^{\alpha}\) où \(C\) est indépendant de \(h\) .

Pour mesurer la vitesse de convergence, on procède comme suit :

  • On réalise des maillages de plus en plus fins du domaine \(\Omega\)

  • Pour chacun des maillages, on détermine la solution \({T}_{h}\)

  • Pour chacun des maillages, on mesure dans une norme bien choisie l’erreur

    ../../../../_images/Object_1117.svg

  • On calcule l’ordre de convergence asymptotique de

    ../../../../_images/Object_2411.svg

vers \(0\) quand \(h\to 0\) par la formule \(-\log(\frac{{e}_{h/2}}{{e}_{h}})/\log(2)\)

  • On vérifie de l’ordre de convergence est en accord avec la théorie.

Par exemple, pour une solution \(T\) suffisamment régulière, en utilisant le norme \({L}_{2}\) , l’ordre de convergence d’un élément linéaire est au minimum de 2 tandis l’ordre de convergence d’un élément quadratique est au minimum de 3.

Aspects informatiques#

Calcul symbolique#

A la lecture des grands principe de la méthode, il est évident que le passage de l’expression analytique de la solution recherchée aux données d’entrées du problème éléments finis à résoudre nécessite de nombreuses opérations de dérivation et d’algèbre. Pour faciliter et automatiser cette étape, on s’appuie des fonctionnalités de calcul symbolique. Elles sont fournies par le module Sympy [bib4] du langage Python. Il s’agit d’un module intégralement écrit en Python d’une grande simplicité d’utilisation.

Pour faciliter la mise en œuvre de la MSM a été développée une classe Python TensorModule qui définit l’objet tenseur ainsi que ses méthodes (produit simplement et doublement contracté, trace, déterminant, …). Nous avons aussi développé les principaux opérateurs s’appliquant sur ces tenseurs (gradient, gradient symétrique, divergence, Laplacien). Sur la base de cette classe, nous avons défini un autre module HookTensor qui permet simplement d’implanter différents tenseurs de Hook (isotrope, anisotrope, orthotrope avec angles nautiques).

Ces classes sont placés dans le répertoire bibpyt/Utilitai des sources de Code_Aster.

Mise en œuvre dans Code_Aster#

Dans cette partie, nous décrivons comme la MSM est mise en œuvre dans le langage de commandes de Code_Aster. Pour ce faire, nous nous basons sur l’exemple de thermique linéaire précédent et détaillons les principales étapes nécessaires. La mise en oeuvre précise est réalisée dans le cas-test tplp107a.

import sympy

On importe le module de calcul symbolique Sympy, pour utiliser ses fonctionnalités.

X,Y=sympy.symbols(“XY”);

On définit les variables d’espace qui seront utilisées pour le calcul symbolique.

T=sympy.Function(“T”); S=sympy.Function(“S”);

On définit le nom des fonctions; \(T\) sera bien évidemment la température et \(f\) le terme source.

T=100*(X**6+Y**6);

On donne l’expression de la solution manufacturée que l’on cherche à retrouver.

S=-Lambda*(sympy.diff(sympy.diff(T,X),X)+sympy.diff(sympy.diff(T,Y),Y));

On déduit l’expression du terme source de l’équation .

SS=FORMULE(NOM_PARA=(“X”,”Y”),VALE=str(S));

Cette commande permet de transformer l’expression précédente en objet formule de Code_Aster qui pourra être évaluée dans le logiciel de la manière suivante:

CLIMIT=AFFE_CHAR_THER_F(MODELE=MO, SOURCE=_F(GROUP_MA=”SURFACE”, SOUR=SS,), );

Bibliographie#

  1. American Institute for Aeronautics and Astronautics, Guide for the Verification and Validation of Computational Fluid Dynamics Simulations , 1998.

  2. American Society of Mechanics Engineers, Guide for Verification and Validation in Computational Solid Mechanics , 2006.

  3. Eric Chamberland, André Fortin, Michel Fortin, Comparison of the performance of some finite element discretizations for large deformation elasticity problems , Computers & Structures, Vol. 88, no. 11-12, pp. 664-673, 2010

  4. Module Sympy de Python, http://code.google.com/p/sympy/ ` <http://code.google.com/p/sympy/>`_