Projet M3N

previous up next contents
Précédent : Logiciels Remonter : Projet M3N, Multi-Modèles et Méthodes Suivant : Actions industrielles

Résultats nouveaux

Modélisation, Optimisation et Contrôle des Fluides Visqueux Turbulents

Projet MICA

 

Participants : Frédéric Hecht , Marie-Hélène Lallemand , Patrick Le Tallec , Claire Martin


Mots-clés : équation de Navier-Stokes, logiciel numérique, maillage adaptatif, mécanique des fluides


Résumé : La mise au point d'un outil d'adaptation de maillage structuré (cartésien) et monobloc (l'adaptation ne concerne qu'une seule grille donnée) a été terminée fin janvier 97 et les logiciels correspondants fournis à CHAM, coordinateur du projet MICA qui devait insérer ce module dans le code Phoenics. Depuis, une stratégie d'interpolation plus précise a été mise au point, et une version multibloc a été développée qui permet de raffiner le maillage par blocs sur plusieurs niveaux de maillage emboités.


On rappelle que ce projet européen a pour but de mettre en place un réseau d'expertise en CFD ( Computational Fluid Dynamics) (mais qui peut être étendu également à d'autres secteurs) auprès des PME en Europe. Les partenaires industriels existants actuellement sont essentiellement concernés par le traitement numérique des problèmes de fours et ceux liés à l'environnement (interne ou externe), et pour cela ont besoin de stratégies de validation et d'adaptation de maillages structurés. La mise au point d'un outil d'adaptation de maillage structuré (cartésien) et monobloc (l'adaptation ne concerne qu'une seule grille donnée) a été terminée fin janvier 97 et les logiciels correspondants fournis à CHAM, coordinateur du projet MICA qui devait insérer ce module dans le code Phoenics.

Par rapport aux choix fixés lors de l'écriture du rapport d'activités de l'an passé, un certain nombre de choses ont été améliorées afin de finaliser l'outil d'adaption monobloc tel qu'il a été fourni début février 97. La première étape a consisté à améliorer les stratégies d'interpolation permettant de calculer les dérivées premières et secondes des champs physiques calculés lors d'une première exécution du code Phoenics. La nouvelle version calcule ces dérivées au centre des faces de chaque cellule par différences finies, en prenant en compte les différents cas particuliers qu'il est possible de rencontrer dans une grille de calcul structurée. Une fois les dérivées partielles secondes ainsi calculées, on a introduit une nouvelle définition d'un hessien pondéré adimensionné à chaque centre des faces des cellules, et pour chaque variable physique calculée. La pondération et l'adimensionnement utilisent une valeur de référence de l'enthalpie et permettent d'obtenir des définitions indépendantes du problème physique étudié, du choix de la variable physique effectué, ou de tout changement d'échelle physique. Une fois le hessien pondéré défini, on le symétrise et on le positivise. La métrique locale est finalement définie en intersectant toutes les métriques locales de chaque variable physique concernée (cf [63]). Cette métrique définit les longueurs des nouvelles cellules dans chaque direction pour le nouveau maillage adapté. Ceci nous permet de définir la taille du nouveau maillage i.e de donner le nombre de cellules à équirépartir dans chaque direction en s'assurant d'avoir un maillage cartésien (cf. [66] pour plus de détails). Par rapport à la stratégie adoptée l'année passée, on y a ajouté un contrôle par l'utilisateur, de la taille du nouveau maillage adapté, et une interpolation plus soignée des variables de sortie du calcul préalable sur le maillage adapté.

Cependant, la stratégie monobloc est d'emploi limité sur des grilles cartésiennes, car elle génére des zones raffinées superflues. C'est pourquoi une politique multibloc a été adoptée afin de déterminer des blocs de raffinement (où l'outil d'adaptation a jugé un raffinement nécessaire) sur éventuellement plusieurs niveaux, correspondant à des ``sous-grilles'' devant être elles aussi cartésiennes. L'adaptation multibloc a d'abord nécessité une généralisation des modules de l'adaptation monobloc en tenant compte des contraintes imposées par CHAM sur la structure multibloc. Ce travail a été initialisé par Claire Martin (ingénieur expert de janvier à mai 97) et Frédéric Hecht (qui a fourni un outil de programmation permettant l'allocation dynamique de mémoire).

Développement numérique du modèle de turbulence

 

Participants : Gorazd Médic , Bijan Mohammadi , Olivier Pironneau , Frédéric Valentin


Mots-clés : mécanique des fluides, turbulence


Résumé : Au cours de l'année 1997, un nouveau modèle de turbulence a été développé à partir des équations pour les tensions de Reynolds. Ce modèle a été introduit dans le solveur bidimensionnel compressible NSC2KE. Afin de préparer la validation de ce nouveau modèle pour les écoulements 3D, on a aussi développé un nouveau solveur des équations Navier-Stokes incompressibles 3D, en utilisant des techniques d'éléments finis stabilisés et un algorithme de projection. Enfin, deux nouveaux cadres mathématiques ont aussi été construits pour analyser des nouvelles lois de parois pour des domaines rugueux : le premier par une méthode de décomposition de domaine (M.D.D.) et le deuxième par une méthode de développement asymptotique à deux échelles.


Dans le cadre d'une collaboration avec le CEA-CESTA, nous avons poursuivi le développement de lois de paroi et modèles de tubulence pour la simulation d'écoulements turbulents instationnaires. Nous considérons toujours le modèle $k-\varepsilon$ et les lois de parois généralisées développés précédemment car cette technique est intéressante par sa simplicité et son faible coût.

Les tests ont été effectués cette année pour différentes configurations, stationnaires et instationnaires, en utilisant le modèle $k-\varepsilon$ et les lois de paroi, implémentées dans le solveur compressible NSC2KE. Une nouvelle implémentation des lois de paroi a été proposée et testée, et on a aussi développé des corrections pour les lois de paroi classique de Reichardt. Dans le cadre du contrat pour le CEA/CESTA, le calcul pour le cas de la marche descendante supersonique ($M_{\infty}=4$) a été réalisé, montrant enfin la nécessité d'introduire des corrections des lois de paroi pour le cas supersonique.

En général, les expériences numériques ont montré que, pour les écoulements 2D, le modèle $k-\varepsilon$ classique était capable de fournir des résultats satisfaisants, même pour les écoulements instationnaires et décollés. En 3D, nous avons constaté que sur un canal tournant les tourbillons secondaires contrarotatifs étaient correctement capturés par la même approche malgré des publications récentes affirmant l'impossibilité d'une telle éventualité. Nous pensons donc que trop souvent les mauvais résultats associés à $k-\varepsilon$ proviennent d'une défaillance de la méthode numérique utilisée ou bien un maillage non adapté. En outre, pour les écoulements 3D, il peut quand même être utile d'introduire des corrections d'ordre 2 aux tensions de Reynolds.

Au cours de l'année 1997, au niveau de la partie de modélisation, un nouveau modèle de turbulence a donc été développé à partir des équations pour les tensions de Reynolds. L'étude mathématique de ce modèle a commencé, et aussi l'implémentation numérique. Le nouveau modèle a été introduit dans le solveur bidimensionnel compressible NSC2KE, et les premiers tests montrent la supériorité de nouveau modèle par rapport au modèle $k-\varepsilon$ classique pour les écoulements 2D du type impinging jet. La validation du nouveau modèle devra être réalisée pour le cas de 3D-curved-duct, pour lequel on a déjà obtenu des premiers résultats avec le modèle $k-\varepsilon$.

Comme la validation de ce nouveau modèle devra être effectuée pour les écoulements 3D, le développement d'un solveur des équations Navier-Stokes incompressibles 3D, plus rapide, a été réalisé. La discrétisation est basée sur les éléments finis $P^1/P^1$ pour la vitesse et la pression. La méthode de projection du premier ordre (en temps) a été utilisée pour le découplage du calcul de la vitesse et de la pression. Le problème de Poisson pour la pression est résolu en utilisant la méthode de gradient conjugué.

Le calcul de vitesse est explicite, avec une stabilisation du terme de convection par le schéma PSI (positive streamwise invariant) du type fluctuation splitting. Les équations du modèle $k-\varepsilon$ sont aussi introduites, en utilisant le même schéma explicite, avec la stabilisation PSI pour la convection. Le solveur a été vérifié pour les cas tests standards : la cavité, la marche descendante, les canaux courbes, en 2D et en 3D. La vérification de l'implémentation du modèle $k-\varepsilon$ et des lois de paroi en 3D a été effectuée pour le cas de 3D-curved-duct (Figure [*]).


  Figure:   Cas test de 3D-curved-duct (calcul préliminaire avec le modèle $k-\varepsilon$ et lois de parois). Tourbillon secondaire à la sortie du coude.

\begin{figure}\centering\mbox{\subfigure[Domaine de calcul] {\epsfig {file=g... ... WL ] {\psfig {file=gorazd-nsike-vb.ps,width=0.25\textwidth}}} \end{figure}


Enfin, pendant l'année 1997, on a développé deux nouveaux cadres mathématiques pour construire et analyser des nouvelles lois de parois pour des domaines rugueux : le premier par une méthode de décomposition de domaine (M.D.D.) et le deuxième par une méthode de développement asymptotique à deux échelles (M.D.A.) ([75], [14], [58], [59], [42] et [41])

a)
L'approche par analyse asymptotique nous a permis de construire des lois de parois pour des problèmes nonlinéaires. En effet, on a pu développer des lois de parois d'ordre un et deux pour les équations de Navier-Stokes contenant des termes de gradient de pression et des termes nonlinéaires liés à la convection. Dans le cadre des équations linéaires, on a construit des lois de parois d'ordre un et deux pour le problème de Laplace, suivies d'une analyse d'erreur. Les lois de parois construites par cette nouvelle méthodologie sont indépendantes du choix de la hauteur, car les constantes d'homogénéisation varient en fonction de la hauteur du couplage. Les validations numériques ont été faites en résolvant plusieurs problèmes de Navier-Stokes dans différents domaines contenant différentes géométries de rugosités. On a remarqué l'importance de l'utilisation des lois de parois d'ordre élevé pour retrouver les bons profils de pression et de $C_f$ à l'interface plaque plane - plaque rugueuse et pour récupérer le décalage du profil de vitesse. On a mis en évidence la nécessité de l'utilisation des lois de parois d'ordre deux pour des problèmes contenant des forts gradients de pressions.
b)
Par la méthode M.D.D. on a construit des lois de parois pour les problèmes de Laplace et Stokes avec différents ordres de viscosité. On a mis en évidence l'équivalence entre la nouvelle approche M.D.D et celle par un développement asymptotique à deux échelles. L'équivalence est obtenue par un simple changement de variable, une fois la cellule tronquée à une certaine hauteur de la paroi. L'influence de la hauteur de couplage sur la qualité de l'approximation est explicitée par l'analyse d'erreur. Les validations numériques ont été faites en résolvant plusieurs problèmes de Stokes définis sur des domaines contenant différentes géométries de rugosités.

Les deux méthodologies sont très flexibles pour traiter différentes formes de rugosités et l'implémentation numérique des lois de parois, prises de façon faible, est immédiate. Le gain sur le maillage, et en conséquence sur le temps de calcul, est sensible. Les lois de parois nous permettent d'utiliser des maillages quatre à cinq fois plus grossiers que ceux utilisés dans un calcul direct, en gardant une excellente précision.

Eléments Finis mixtes pour les Fluides Incompressibles

 

Participants : Fadi El Dabaghi , Mohamed Amara


Mots-clés : élément fini, équation de Navier-Stokes, mécanique des fluides


Résumé : On a proposé et développé cette année une nouvelle formulation mixte $(\psi-\omega)$ en 2D/3D pour le Bilaplacien $\Delta^2$ avec une méthode d'approximation optimale de classe $C^0$ par éléments finis $P^k$ et une estimation d'erreur en $O(h^k)$ pour $k\ge 1$.


Dans un cadre variationnel naturel $\psi\in H^1(\Omega)$ et $\omega\inH^{-1}(\Delta,\Omega)$, sans condition de régularité supplémentaire, on a proposé et développé cette année une nouvelle formulation mixte $(\psi-\omega)$ en 2D/3D pour le Bilaplacien $\Delta^2$ avec une méthode d'approximation optimale de classe $C^0$ par éléments finis $P^k$ et une estimation d'erreur en $O(h^k)$ pour $k\ge 1$. L'originalité de ce résultat tient dans l'optimalité de l'approximation tandis que les nombreuses études sur ce sujet ne donnaient qu'une estimation en $O(h^{k-1})$ pour des éléments finis $P^k$ avec $k\ge 2$; cette nouvelle formulation par son coût peut être une alternative sérieuse à celle en $u-p$ car elle est de coût égal en 2D et inférieur en 3D. Cette étude s'est appliquée à l'analyse théorique et à la validation numérique bidimensionnelle des équations d'Euler incompressibles [60], à l'analyse théorique et à la validation numérique bidimensionnelle pour un problème de plaque avec fissure.

Optimisation et Contrôle

  La conception optimale de formes demande à utiliser les modèles numériques de comportement développés au préalable pour optimiser la forme de l'objet étudié. Enfin, plus généralement, les techniques de contrôle peuvent être utilisées comme outil de résolution ou de stabilisation de certaines équations.

Optimisation de formes en Aérodynamique Externe

 

Participants : Emmanuel Laporte , Patrick Le Tallec , Bijan Mohammadi


Mots-clés : algorithme numérique, mécanique des fluides, optimisation de forme


Résumé : Après avoir formalisé mathématiquement les différentes stratégies de calcul des gradients des fonctions coûts par rapport à la forme de l'objet à optimiser, une nouvelle méthode d'optimisation de formes a été développée permettant de pouvoir adapter le maillage de manière arbitraire au cours de la boucle d'optimisation.


Sous forme générale, le problème d'optimisation de formes s'écrit

\begin{displaymath}\mbox{Trouver une forme }\, \gamma^* \, \mbox{ telle que } ... ...^*\right)\le 0 \\ h\left(\gamma^*\right)=0 \end{array}\right.\end{displaymath}


où la fonction coût $j$ dépend de la forme $\gamma$ et d'une variable d'état $ W_\gamma$, c'est à dire, $j\left(\gamma\right)=J\left(\gamma,W_\gamma\right)$, $W_\gamma$étant la solution d'une Equation aux Dérivées Partielles (Euler, Navier-Stokes, ...) du type


\begin{displaymath}E\left(\gamma,W_\gamma\right)=0,\end{displaymath}


$g$ et $h$ étant des contraintes géométriques ou physiques.

Les travaux effectués cette année ont permis de formaliser mathématiquement le calcul du gradient de la fonction $j$ par rapport à $\gamma$ et de montrer que le gradient de la fonction approchée $j_h$ par rapport au maillage et obtenu par différentiation automatique etait une approximation, dans les cas réalistes, du même ordre que l'approximation de $W_\gamma$ obtenue par les solveurs classiques. On a également montré la différentiabilité de la fonctionnelle dissipation d'énergie pour les équations de Navier-Stokes incompressibles et stationnaires, et pour celles de Stokes instationnaires.

Au vu de ces résultats, on utilise donc maintenant une nouvelle méthode d'optimisation de formes consistant à minimiser la fonction $j$ en ne connaissant qu'une approximation de sa valeur et de celle de son gradient, la forme étant déterminée par un nombre fixe de paramètres. Cette méthode permet de travailler indépendamment du maillage au cours de l'optimisation, et donc de pouvoir adapter le maillage pour obtenir les meilleures approximations possibles. Pour l'adaptation, on utilise le logiciel BAMG, développé par le projet GAMMA. Le passage au couplage adaptation de maillage et optimisation de formes est donc immédiat et des simulations numériques sont en cours sur des problèmes stationnaires. Une des principales questions reste encore la stabilité des algorithmes de minimisation de type point intérieur  [21] lorsque le gradient utilisé n'est qu'une approximation.


  Figure:   Iso-mach après optimisation : la cambrure du profil et le choc intrados ont quasiment disparu.

\begin{figure} \begin{center} \mbox{\psfig {figure=laporte-fig1.eps,width=100mm}} \end{center} \end{figure}


Des problèmes instationnaires sont aussi à l'étude, tels que le couplage aéroélastique, ou le buffeting, ou des problèmes de stabilité lors de variations brusques de conditions de vol. L'étude des problèmes tridimensionels reste limitée par le coût mémoire des problèmes instationnaires et par les développements encore peu avancés de l'adaptation de maillage en 3D.

Equations de Maxwell 3D

 

Participants : Paul Ayoub , Marie-Odile Bristeau , Jacques Périaux


Mots-clés : algorithme numérique, élément fini, équations de Maxwell, logiciel numérique, modélisation, contrôle


Résumé : Un solveur explicite, entièrement écrit en C++, a été développé pour la simulation numérique de la diffraction d'ondes électromagnétiques par des obstacles tridimensionnels, parfaitement conducteurs et nonconvexes. En parallèle, diverses modifications ont été apportées pour améliorer l'efficacité d'une formulation "Contrôlabilité Exacte" qui permet d'obtenir la solution harmonique en temps du problème à partir du calcul explicite de diverses solutions instationnaires.


On s'intéresse à la résolution numérique des équations de Maxwell en trois dimensions d'espace, pour la simulation de phénomènes de diffraction d'ondes électromagnétiques par des obstacles que l'on suppose dans un premier temps parfaitement conducteurs et nonconvexes. Une première approche a consisté à développer un solveur explicite, utilisant une formulation de second ordre en temps et en espace. Le champ électrique et le multiplicateur de Lagrange sont liés par l'intermédiaire d'une équation des ondes vectorielle pour le premier et une équation de la chaleur pour le deuxième. Après avoir validé le code (écrit entièrement en C++) et le schéma sur un conducteur parfait sphérique, nous avons testé le comportement de l'algorithme sur différentes géométries réalistes. Les résultats obtenus, dans les cas d'une onde incidente harmonique ou d'une impulsion, sont satisfaisants. Cette méthode est maintenant en cours d'adaptation pour pouvoir traiter cette année la diffraction d'ondes électromagnétiques par un conducteur cylindrique de petit rayon (un centième de longueur d'onde), comme par exemple, les antennes.

On s'intéresse aussi au cas du régime harmonique. Dans des études précédentes, on a montré, pour les cas 2D et 3D acoustiques, l'efficacité d'une formulation "Contrôlabilité Exacte" pour accélérer la convergence de la solution instationnaire, calculée par un schéma explicite, vers la solution du régime harmonique. Cette méthodologie a été étendue aux cas vectoriels 2D et 3D. La formulation moindres carrés introduite conduit à la résolution du problème par un algorithme de gradient conjugué préconditionné. Nous avons amélioré l'efficacité du code principalement sur deux points :

- Le calcul du terme de la forme

\begin{displaymath}\int_\Omega [(\nabla\times u)(\nabla\times v) + (\nabla\cdotu)(\nabla\cdot v)] dx \end{displaymath}


qui intervient aussi bien dans les équations des ondes que dans l'opérateur de préconditionnement.

- La résolution du système linéaire associé au préconditionnement. Le code étant écrit pour une machine vectorielle (CRAY90), le calcul du terme ci-dessus par un produit matrice-vecteur où la matrice creuse serait stockée sous forme morse par ligne n'est pas du tout intéressant. On a donc introduit un stockage par diagonale ce qui conduit à de longs vecteurs et donc une très bonne vectorisation.

En ce qui concerne la résolution du système linéaire, on se trouve dans la situation où l'on doit résoudre successivement plusieurs systèmes avec matrice identique et seconds membres différents. Ce travail a été poursuivi en collaboration avec J. Erhel et F. Guyomarc'h du projet Aladin. La méthode est basée de nouveau sur un algorithme de gradient conjugué ; un sous-espace de Krylov est généré à la première résolution. Lors des résolutions suivantes, on applique un gradient conjugué modifié dans lequel la solution initiale et les directions de descente sont modifiées par des projections orthogonales sur le sous-espace de Krylov. Les travaux de J. Erhel et F. Guyomarc'h ont permis de définir un algorithme robuste sans pertes d'orthogonalité. Il a été testé sur des problèmes issus de calculs 3D acoustiques ou électromagnétiques.

Le code de résolution des équations de Maxwell a été validé et testé sur différentes géométries (ex : calcul dans et autour d'une entrée d'air). Par ailleurs, le code a été parallélisé sur CRAY T3E par T. Rossi (Université de Jyvaskyla et Dassault). Etant donné le caractère complètement explicite du code, la parallélisation est très efficace.

Couplage Fluide-Structure

 

Participants : Rodolfo Araya , Emmanuel Laporte , Patrick Le Tallec , Bijan Mohammadi , Mugurel Stanciu , Marina Vidrascu


Mots-clés : algorithme numérique, décomposition de domaines, équation de Navier-Stokes, calcul d'erreurs, calcul des structures, turbulence


Résumé : Le travail de cette année s'est organisée autour de trois actions. La première action a consisté à mettre au point un calcul d'erreur a-posteriori par résidus pour valider les calculs de déformations de matériaux fortement hétérogènes. La seconde action a porté sur la définition et les premières simulations numériques d'écoulements sanguins dans des anévrismes. La dernière action s'est intéressée à la simulation directe des phénomènes de buffeting.



  Figure:   Comparaison entre les solutions obtenues avec le maillage initial, le maillage adapté à l'aide de l'estimateur de résidu et la solution de référence sur un problème d'élasticité test.

\begin{figure}\begin{center}\mbox{\psfig {figure=araya-iso2_u1_3.ps,width=0.6\textwidth}}\end{center}\end{figure}


Le développement de calcul d'erreur pour la partie structure s'est effectué à deux niveaux. D'une part, en coopération avec Saloua Mani de l'Université de Tunis, nous avons formalisé un calcul d'erreur a priori garantissant la convergence au premier ordre d'une approximation des modèles de coques géométriquement exactes par des éléments finis de type DKT [23]. La seconde contribution a consisté à définir, analyser, implémenter et valider numériquement un nouvel estimateur d'erreur a posteriori du type résidu pour des problèmes elliptiques non homogènes (voir [61]). Cet estimateur est facile à implémenter, prend en compte le caractère hétérogène des coefficients, et surtout, il est correctement dimensionné par rapport aux données physiques. Nous avons appliqué notre estimateur aux problèmes d'élasticité avec coefficients élastiques fortement hétérogènes. L'estimation d'erreur ainsi obtenue peut être fournie à un mailleur adaptatif et permet alors d'obtenir à moindre coût et de manière automatique une solution fiable (Figure [*]) Les premières extensions prévues sont maintenant de passer au cas trimensionnel, incompressible et/ou anisotrope.


   Figure:  Champ de vitesse et de pression associés à des lignes de courant dans un modèle d'anévrisme sacculaire (rapport entre la hauteur de l'anévrisme et le diamètre du tronc de 3). L'anévrisme est situé à l'apex d'un embranchement asymétrique sur un coude (coude à 90 degrés, rapport de courbure de 1/10). Le col étroit couvre la paroi externe du coude et la paroi interne de la branche incurvée par la malformation (rapport entre la hauteur du col et le diamètre du tronc de 0.2 et rapport entre le diamètre du col et celui du tronc de 1.6). Ecoulement laminaire (Re=500).

\begin{figure}\begin{center}\mbox{\psfig {figure=THIRIET1_bis.ps,width=0.75\t... ...gure=THIRIET2_bis.ps,angle=-90.,width=0.75\textwidth}} \end{center}\end{figure}

Figure:  Champ de vitesse et de pression associés à des lignes de courant dans un modèle d'anévrisme artériel cérébral situé à l'apex d'une bifurcation symétrique (rapport des aires des sections droites de 0.72, rayon de courbure de la zone de transition de 1/4). Tronc rectiligne. Anévrisme sphérique (rapport entre le rayon de l'anévrisme et celui du tronc de 2) possédant un col large (rapport entre la longueur du col et son rayon de 1.8 et rapport entre la hauteur du col et le diamètre du tronc de 0.5). Ecoulement laminaire (Re=500).

Le problème de la simulation des anévrismes a été déjà décrit en section 3. Au cours du travail effectué cette année, deux configurations génériques ont été définies, et étudiées numériquement : (1) un anévrisme développé à l'apex d'une bifurcation symétrique type tronc basilaire, et (2) un anévrisme développé à l'apex d'un embranchement asymétrique entre un tronc courbe et une branche de calibre plus faible. Les principaux résultats obtenus sont résumés sur les figures [*], [*].

Enfin, le but du travail réalisé sur la simulation du buffeting était de montrer les capacités du code NSC2KE à simuler le phénomène d'oscillation instationnaire de l'onde de choc en écoulement transsonique, provoqué par une interaction instationnaire entre choc et couche limite turbulente. La première partie a consisté à prouver que ce phénomène pouvait s'observer même si le profil d'aile était parfaitement rigide et immobile. De bons résultats ont été obtenus en construisant des maillages adaptés au phénomène physique et en utilisant un modèle de turbulence de type $k$ -$\varepsilon$ couplé à une loi de paroi. Afin ensuite de comprendre l'influence possible de la déformation de la structure, on a également ajouté au code un modèle simplifé de couplage fluide-structure. L'optimisation de formes autour du buffeting avait été aussi envisagée, mais cette étude a également montré que la simulation numérique directe du phénomène restait encore trop exigeante pour pouvoir obtenir des résultats suffisament fiables et rapides susceptibles d'être utilisés à l'intérieur d'une boucle d'optimisation.

Semi-conducteurs

 

Participants : Americo Marrocco , Benoit Perthame , Philippe Montarnal , Frédéric Hecht


Mots-clés : algorithme numérique, élément fini, logiciel numérique, modélisation, semiconducteur


Résumé : Les recherches sur la simulation numérique des dispositifs semi-conducteurs se sont tout d'abord poursuivies selon les axes évoqués l'an passé : couplage multiéchelles d'éléments finis mixtes, adaptation de maillage, développement de nouveaux solveurs algébriques. Par ailleurs, nous avons étudié et développé différents modèles permettant de prendre en compte l'évolution énergétique des différents porteurs, en regardant en particulier les liens existants entre modèles de transport d'énergie et modèles hydrodynamiques. Ceci nous a permis d'aborder l'étude de nouveaux dispositifs : les H.E.M.T.


Nous nous intéressons toujours à la simulation numérique des dispositifs semi-conducteurs de type III-V gouvernés par les modèles de dérive-diffusion ou par les modèles de type hydrodynamiques (énergie-transport, hydrodynamique simplifié).

Les dispositifs qui nous intéressent sont d'une part les dispositifs bipolaires comme le T.B.H. (Transistor Bipolaire à Hétérojonction) mais aussi des dispositifs de conception plus récente comme les H.E.M.T. (High Electron Mobility Transistor); ces derniers sont essentiellement du type unipolaire mais les dimensions réduites de ces dispositifs rendent nécessaire l'utilisation des modèles de type hydrodynamique. Les recherches sur la simulation numérique des dispositifs semi-conducteurs se sont tout d'abord poursuivies selon les axes évoqués l'an passé.


  Figure:   Maillage initial ayant servi pour des calculs antérieurs (en haut)et maillage adapté (pour le problème d'équilibre) obtenu après 3 itérations.

\begin{figure}\begin{center} \mbox{\psfig {file=am-mail1.ps,width=0.6\textwid... ...ile=am-mail2.ps,width=0.6\textwidth,angle=-90.,clip=}}\end{center} \end{figure}


Au début de cette année a aussi été entreprise une étude concernant un nouveau type de dispositif: les H.E.M.T. La première structure (HEMT-pseudomorphique) nous a été proposée par l'I.E.F. (Institut d'Electronique Fondamentale d'Orsay -Laleh Rochebois et André de Lustrac) puis un autre dispositif du type HEMT nous a été proposé par le CNET. Pour simuler ce type de dispositif, il a été nécessaire d'introduire une modélisation des contacts de type SCHOTTKY; pour ce genre de dispositifs, les résultats provenant du modèle de dérive-diffusion sont très différents de ceux donnés par l'énergie-transport (ou l'hydrodynamique simplifié). Les courants traversant le dispositif sont sous évalués par le modèle de dérive-diffusion.

Enfin, dans l'optique d'effectuer des tests de validation, ou de pouvoir calculer des solutions de référence (sur des maillages suffisamment fins), nous avons essayé l'an passé de mettre en oe uvre des techniques de résolution moins gourmandes en place mémoire (gradient-conjugué préconditionné, GMRES préconditionné) mais les résultats obtenus n'ont pas été très satisfaisants. Une autre voie a été entamée, en collaboration avec le CERFACS (L.Giraud) en s'appuyant sur les techniques de décomposition de domaines et de résolution sur machines multiprocesseurs ou réseau de stations.

Le dernier volet des études de cette année concerne le développement et la mise au point de modèles prenant en compte les équations d'énergie des porteurs, modèles dont la nécessité se fait sentir quand la taille du dispositif à étudier diminue. Les résultats obtenus dans ce cadre cette année sont les suivants [13], [53]:

Ensuite, nous nous sommes intéressés aux liens existants entre les modèles de type transport d'énergie et ceux de type hydrodynamique. Pour cela, nous étudions la dérivation d'un modèle hydrodynamique à partir d'une équation de transport de Boltzmann contenant différentes échelles de collision. Nous effectuons un développement de type Chapman-Enskog en exploitant les deux échelles de temps des termes de collision élastique et électrons/électrons. Le modèle obtenu généralise le modèle hydrodynamique classique dans la mesure où nous faisons apparaître deux termes supplémentaires de "force de friction thermique" et de "flux de chaleur de friction". De plus, suivant le choix des paramètres du modèle, nous retrouvons à la fois le modèle de transport d'énergie et le modèle hydrodynamique d'Euler-Poisson. Enfin, grâce à la continuité du passage entre le modèle de Boltzmann et le modèle hydrodynamique proposé, nous déduisons de l'entropie $H$ de Boltzmann une loi d'entropie pour le modèle hydrodynamique. De cette loi d'entropie, nous obtenons, grâce aux résultats de Kawashima et Shizuta, un jeu de variables qui symétrise le système hydrodynamique. Ces variables, que nous appelons "variables entropiques", sont d'un grand intérêt sur le plan numérique [31].

Analyse Particulaire des Fluides Complexes

Mots-clés : équation de Navier-Stokes, équation de Boltzmann


Fluides multiespèces ou diphasiques

 

Participants : Jean-François Bourgat , Patrick Le Tallec , Bijan Mohammadi , Benoit Perthame


Résumé : Le premier thème de cette action porte sur la simulation de la séparation de gaz de masses différentes par centrifugation et circulation thermique. Le second aspect du travail a consisté à étudier le mouvement de bulles dans un liquide avec prise en compte de l'effet de masse ajoutée.


Le premier thème de cette action a été initié en 1997 dans le cadre d'un partenariat industriel et porte sur le problème de séparation de gaz de masses différentes par centrifugation et circulation thermique. D'un point de vue scientifique, il s'agit de résoudre les équations de Navier-Stokes compressibles avec une loi de viscosité fortement nonlinéaire variant sur dix ordres de magnitude. Pendant la première phase du travail, nous avons procédé à une réécriture des équations permettant à la fois une simplification et une généralisation du problème. Ensuite, nous avons développé un outil informatique pour résoudre ces équations dans une configuration particulière qui est celle utilisée dans l'industrie. Une conclusion de cette analyse est que l'approche classique incompressible reste valable pour ce type d'applications dans la plus grande partie du domaine de calcul. Mais un calcul compressible semble inévitable dans les régions où l'écoulement est perturbé par le dispositif d'extraction. Nous nous orientons maintenant vers l'analyse plus fine de l'écoulement dans ces régions.

Le second aspect du travail a consisté à étudier, en collaboration avec B. Lucquin (Paris 6), le mouvement de bulles dans un liquide avec prise en compte de l'effet de masse ajoutée qui est important dans ce cas. C'est un nouveau sujet de recherche qui débute.

Les bulles sont supposées sphériques, seules les forces de pression sont prises en compte et l'on suppose de plus qu'elles n'ont pas de collision. La modélisation de la dynamique des bulles repose sur une équation cinétique de type Vlasov portant sur la densité $f(t,x,p)$ de bulles d'impulsion $p$, à la position $x$, au temps $t$. L'ensemble de toutes les interactions entre toutes les bulles est pris en compte. Des modèles simplifiés sont obtenus grâce à un développement à l'ordre 1 ou 2 par rapport au paramètre naturel qu'est le taux volumique de bulles $\lambda$. L'approche numérique est faite par une méthode particulaire standard à partir du Hamiltonien discret et les équations obtenues au niveau particulaire sont résolues par un schéma de Runge-Kutta d'ordre 4. Le potentiel est calculé par une formule analytique sommant les potentiels de chaque bulle. Une matrice de dérivées secondes d'un noyau $K =1/4\pi\vert x\vert$ nécessite l'utilisation d'une régularisation pour les faibles valeurs de $\vert x\vert$ ($ <$ rayon des bulles).

Des résultats numériques ont été obtenus avec 900 bulles en 2D et 300 en 3D. Le nombre d'opérations proportionnel au carré du nombre de bulles et le suivi de trajectoires nécessitant des milliers de pas de temps ne permettent pas de faire mieux pour l'instant. Cependant 100 bulles suffisent déjà pour rendre les phénomènes d'agrégation visibles. Les bulles se meuvent lentement en s'agrégeant en couches qui à leur tour se regroupent en amas. Notre modèle cesse alors d'être valable, il faudrait introduire des collisions. Actuellement nous testons une formulation cinétique en position-vitesse au lieu de position-impulsion. Elle donne des résultats voisins et devrait permettre des développements du modèle vers les collisions de bulles et les réflexions sur des parois. La modélisation présentée ci-dessus s'applique aussi au mouvement de particules solides dans un liquide (sédimentation) bien que l'effet de masse ajoutée soit relativement moins important.

De Boltzmann à Navier-Stokes

 

Participants : Jean-François Bourgat , Patrick Le Tallec , François Mallinger , Benoit Perthame , Jean-Philippe Perlat


Résumé : Nous avons étudié numériquement différents outils pour identifier, à partir d'une solution de l'équation de Boltzmann, les régions où les modèles continus sont valides et celles où l'équation de Boltzmann doit être utilisée. Nous avons aussi analysé la cohérence mathématique de différentes approximations de l'équation de Boltzmann proposées par D. Levermore, et leurs liens avec les équations de Navier-Stokes. Grâce à cette étude, nous avons obtenu une nouvelle interprétation cinétique des équations de Navier-Stokes compressible en terme de limite asymptotique des équations de Levermore. Nous avons alors proposé et implémenté une stratégie de couplage explicite en temps, qui résout les équations de Levermore dans un domaine proche de l'obstacle, et les équations de Navier-Stokes dans un domaine fluide prédéfini. Les derniers résultats sur les Equations de Boltzmann concernent enfin le travail effectué en commun avec l'Institute of Theoretical and Applied Mechanics de Novosibirsk sur le développement et l'analyse de modèles de collision plus riches prenant en compte les différents niveaux d'énergie vibrationnelle.


Le but de cette étude est développer des modèles ou des techniques pour étudier la transition entre des écoulements raréfiés qui doivent être décrits par des modèles cinétiques appropriés, et des écoulements plus continus susceptibles d'être décrits par des modèles de la mécanique des milieux continus.

Nous avons étudié numériquement des outils pour identifier, à partir d'une solution de l'équation de Boltzmann, les régions où les modèles continus sont valides et celles où l'équation de Boltzmann doit être utilisée. Le plus précis est la comparaison de la distribution des vitesses avec la maxwellienne ou la distribution de Chapman-Enskog, mais ce calcul long ne peut être que ponctuel. Près des parois, le calcul de la vitesse tangentielle et du saut de température qui devraient être nuls pour un gaz en équilibre sont de bons indicateurs du déséquilibre. Dans l'ensemble de l'écoulement, les bonnes mesures sont fournies par le tenseur des contraintes visqueuses et par le flux de chaleur dont les normes croissent avec le déséquilibre. Enfin, pour les gaz diatomiques, l'écart entre la température de translation et la température de rotation apporte une mesure complémentaire du déséquilibre.

Après identification des différentes régions, il faut ensuite pouvoir utiliser un modèle numérique adapté dans chaque région identifiée. Nous avions décrit et validé numériquement l'année dernière le problème aux limites obtenu en associant le modèle cinétique asymptotique aux 14 moments développé par D.Levermore, et des conditions aux limites faibles écrites en terme de demi-flux. Cette année, nous avons d'abord démontré la cohérence mathématique des solutions ainsi obtenues. L'approche utilisée consiste à vérifier l'existence et l'unicité de la solution du problème aux limites linéarisé par une méthode de viscosité.


  Figure:   Isovaleur de la densité (en haut). Coupe du Mach (à gauche) et de la Temperature (à droite). Comparaison avec une méthode de Monte Carlo résolvant les équations de Boltzmann.

\begin{figure}\begin{center}\mbox{ \psfig {file=denscoupl.ps, width=\textwidt... ...\epsfig {figure=comptemp1d.eps,width=0.49\textwidth}}\end{center}\end{figure}


Nous avons également étudié un développement asymptotique en puissance du temps de relaxation des équations aux 14 moments non-conservatives. Grâce à cette étude, nous avons obtenu une nouvelle interprétation cinétique des équations de Navier-Stokes compressible en terme de limite asymptotique des équations de Levermore. Il en découle une écriture naturelle des conditions d'interface pour le couplage de ces deux systèmes d'équations. Nous avons alors proposé et implémenté une stratégie de couplage explicite en temps, qui résout les équations de Levermore dans un domaine proche de l'obstacle, et les équations de Navier-Stokes dans le domaine fluide prédéfini. Les configurations typiquement envisagées pour la validation du couplage sont des écoulements monoatomiques, froids, autour de plaques à divers incidences. Le cadre ainsi choisi correspond aux écoulements en soufflerie étudié par le CESTA et par l'ONERA, et pour lesquels on dispose aussi de résultats numériques obtenus à l'INRIA par une résolution précise des équations de Boltzmann. Nous avons alors comparé de façon systématique les résultats du problème couplé à ceux obtenus par le modèle cinétique des équations de Boltzmann (voir figure [*]).

Les derniers résultats sur les Equations de Boltzmann concernent le travail effectué en commun avec l'Institute of Theoretical and Applied Mechanics de Novosibirsk et dans le cadre du stage postdoctoral de François Mallinger sur le développement et l'analyse de modèles de collision plus riches. Pour calculer des écoulements à hautes températures nous avons implémenté un nouveau modèle de collision, introduit par le groupe du Pr. Ivanov de l'ITAM, qui tient compte des échanges d'énergie internes de type VT et VV entre particules [70]. Deux modélisations des processus de relaxation vibrationnelle dans un gaz ont été aussi développés: le modèle de Landau-Teller (LT) et les équations master (ME). Différents modèles de taux de relaxation VT et VV ont été étudiés et proposés. Un algorithme a été développé pour la résolution des équations de relaxation (LT ou ME) couplées aux équations de conservation de la quantité de mouvement et de l'énergie totale. Les deux modèles de relaxation ont été comparés pour le calcul d'un écoulement autour d'un cylindre infini, pour différents nombres de Mach [71].


previous up next contents Précédent : Logiciels Remonter : Projet M3N, Multi-Modèles et Méthodes Suivant : Actions industrielles