La modélisation tridimensionnelle et la simulation numérique avec REM3D

La modélisation tridimensionnelle et la simulation numérique avec REM3D

Le logiciel REM3D

Le logiciel REM3Dr développé au CEMEF par [PICHELIN, 1998], [DABOUSSY, 2000], [BIGOT, 2001], [BATKAM, 2002], [ROCHA DA SILVA, 2004] pour la simu75 elisation tridimensionnelle et la simulation num´ Chapitre 5. La modélisation tridimensionnelle et la simulation numérique avec REM3Dr lation de l’injection de thermoplastiques est l’outil que nous allons utiliser afin de simuler notre procédé. Ce logiciel est constitué de trois principaux solveurs qui forment sa base numérique : un solveur mécanique de Stokes, un solveur de transport et un solveur thermique. ♦ Le solveur de Stokes Destiné à la résolution du probl`eme mécanique de Stokes en vitesse et pression, il est basé sur une méthode éléments finis mixtes utilisant l’élément stable P1 + /P1 ([COUPEZ, 1996]). La résolution se fait par une méthode itérative de type résidu conjugué préconditionné. Dans le cadre de la version multidomaine du logiciel, différents solveurs ont vu le jour, chacun étant associé à un domaine pris en compte dans une simulation. Ainsi, le code poss`ede à ce jour : X un solveur fluide viscoplastique dans lequel sont résolues les équations de Stokes pour les domaines fluides ; X un solveur gaz : le gaz est représenté par un fluide compressible de faible viscosité ; X un solveur outil permettant de traiter certains aspects du couplage fluide / structure. Ces outils peuvent ˆetre rigides, faiblement compressibles ou élastiques. Le fluide est considéré comme un fluide newtonien incompressible. Les inconnues du probl`eme sont la vitesse et la pression notées respectivement v et p. Les tenseurs des vitesses de déformation et des contraintes sont définis par :  = 1 2 (∇v + ∇v t ) σ = 2η(|(v)|)(v) − p11 (5.1) La viscosité dépend du taux de cisaillement γ¯˙ = s 2 X i,j ˙i,j suivant une loi qui peut ˆetre de la forme : η(γ¯˙) = η0((c + λ 2 γ¯˙ 2 ) (m−1)/2 (5.2) Un comportement newtonien correspond à m = 1, une loi puissance est obtenue en prenant c = 0 et c = 1 donne une loi de Carreau. Finalement, pour des fluides incompressibles, les équations du probl`eme sont, dans le sous-domaine Ωf : ( 2∇ · (η(|(v)|)(v) − ∇p = 0 ∇ · v = 0 (5.3) 76 5.1. Le logiciel REM3Dr ♦ Le solveur équation de transport Il permet de déterminer l’avancée du front de mati`ere. A chaque sous-domaine est ` associée une fonction caractéristique, SΩi , définie par : SΩi (x, t) = ( 1 si x ∈ Ωi(t) 0 si x 6∈ Ωi(t) ∀ t (5.4) Elle est obtenue à chaque incrément de temps par résolution de l’équation de transport : dSΩi dt = ∂SΩi ∂t (x, t) + v · ∇SΩi (x, t) ∀x ∈ Ω, ∀t ∈ R + (5.5) Dans le cas de l’utilisation d’une description ALE (Arbitrary Eulerian-Lagrangian), une vitesse arbitraire u est associée à chaque point géométrique x. L’équation de transport doit prendre en compte ce mouvement à travers l’équation : dSΩi dt = ∂SΩi ∂t (x, t) + (v − u) · ∇SΩi (x, t) ∀x ∈ Ω, ∀t ∈ R + (5.6) Il apparaît dans l’équation (5.6) que si la vitesse arbitraire u co¨ıncide avec la vitesse du fluide, la variation spatiale de la fonction caractéristique sera nulle : c’est une formulation lagrangienne. Au contraire, si u est nulle, nous serons dans le cas d’une description eulerienne. En pratique, dans le cas qui nous intéresse, nous choisirons une vitesse de maillage égale à la vitesse du fluide dans la direction de déformation, et nulle dans les deux autres directions.

Intérˆet de l’utilisation de REM3D

Bien que ce logiciel ait été initialement con¸cu pour la simulation en injection des polym`eres, sa flexibilité le rend adaptable à un bon nombre d’applications : co-extrusion, écoulements avec interaction fluide-structure [BATKAM, 2002], fonderie [SAEZ, 2003], expansion de mousses . Dans le cas qui nous préoccupe, c’est la capacité de REM3Dr à gérer les domaines multiples, les interfaces et les surfaces libres qui va nous intéresser. A cela, s’ajoute la ` possibilité d’adapter le maillage ce qui permet de mieux capturer les interfaces. 

Résolution du probl`eme de Stokes pour des fluides incompressibles 

Formulation variationnelle du probl`eme de Stokes 

Comme nous l’avons vu en introduction, le code REM3Dr est constitué d’un solveur de Stokes. Il permet de résoudre les équations mécaniques écrite sous la forme de Stokes généralisée dans le cas incompressible et compressible [ROCHA DA SILVA, 2001]. Pour un domaine fluide incompressible, le syst`eme de Stokes a été introduit précédemment (5.3). Si la fonction caractéristique du fluide est connue à l’instant t, il est possible de déterminer le sous domaine Ωf contenant le fluide. La fronti`ere Γ de Ω est divisée en deux parties : Γm correspondant à la fronti`ere avec le moule et Γi qui est le plan d’entrée (dans le cas de l’injection). Nous avons alors Γ = Γm ∪ Γi et nous prenons un champ de vitesse virtuel w ∈ H 1 Γ (Ω)3 défini par, dans le cas newtonien : H 1 Γ (Ω)3 =  w ∈ H 1 Γ (Ω)3 : w = 0 sur Γm, w − (w · n) n = 0 sur Γi (5.9) Ainsi, le probl`eme à résoudre est : trouver (v, p) tels que :    Z Ωf 2η(|(v)|)(v) : (w)dΩ − Z Ωf p ∇ · w dΩ = Z Ωi piw · n dΓ − Z Ωf q∇ · v dΩ = 0 ∀(w, q) ∈ H 1 Γ (Ωf ) 3×L 2 (Ωf ) (5.10) O`u L 2 (Ωf ) =  v : Ω → R; Z Ω v 2 dx < +∞  Ce syst`eme 5.10 peut se généraliser sur tout le domaine en multipliant les différents termes par la fonction caractéristique des fluides. 5.1.3.2 Formulation dans le domaine vide Le domaine vide est le domaine complémentaire du (ou des) domaine(s) fluide(s) dans la cavité. Les champs de vitesse et de pression étant définis uniquement dans le domaine 79 Chapitre 5. La modélisation tridimensionnelle et la simulation numérique avec REM3Dr fluide, ils sont prolongés dans la partie vide de mani`ere à ne pas perturber le fluide. La pression est nulle dans le vide mais non nulle à l’interface fluide/vide. La vitesse est continue dans le domaine de calcul, elle est donc prolongée dans le vide en posant : div(2ηv(v)) = 0 ; ηv étant la viscosité du vide, grandeur purement numérique. Ainsi, avec le mod`ele compressible introduit par [ROCHA DA SILVA, 2004] la formulation variationnelle dans le domaine vide s’écrit : Z Ω 2ηv (v) : (w) dΩ − Z Ω p ∇ · w dΩ = 0 ∀ w ∈ H 1 (Ωv) 3 Z Ω q ∇ · vdΩ + χ Z Ω pq dΩ = 0 ∀q ∈ L 2 (Ωv) (5.11) 

L’adaptation de maillage 

Le module d’adaptation de maillage a initialement été développé dans REM3Dr par E. Bigot. Dans le cas statique, le probl`eme est le suivant et est illustré sur la figure (5.2) : nous voulons déformer le maillage initial aux fronti`eres d’un sous-domaine (carré sur la figure 5.2(a)) afin de bien visualiser l’interface entre les deux domaines figure 5.2(b). (a) Maillage initial (b) Maillage adapté Fig. 5.2 – Principe de l’adaptation de maillage

L’adaptation de maillage

La méthode

Pour définir le sous-domaine Ωf , une fonction caractéristique lui est associé et est définie sur Ω par : 11Ωf (x, t) = ( 1 si x ∈ Ωf 0 si x 6∈ Ωf ∀ t (5.12) En notant 11h Ωf = π h (M)11Ωf la projection de la fonction caractéristique sur le maillage, le probl`eme peut se mettre sous la forme [BIGOT, 2001] : min M kπ h (M)11Ωf − 11Ωf k (5.13) La méthode d’adaptation de maillage est une méthode de type r-adaptation ce qui signifie que les éléments sont déformés mais que la topologie de maillage reste inchangée au cours de la simulation. L’erreur exacte d’interpolation commise peut s’écrire [BIGOT, 2001] :  = k11h Ωf − 11Ωf k0,K =   X K∈Th(Ω) (1 − 11K)11K|K|   1 2 (5.14) o`u 11K est la valeur de l’interpolé sur l’élément K. En observant l’expression 5.14, il apparaît que seuls les éléments dont la fonction caractéristique est ni 0 ni 1 contribuent à l’erreur. Ces éléments sont ceux qui sont traversés par l’interface entre les deux domaines. Leur contribution à l’erreur est proportionnelle à leur volume |K| et la valeur de la fonction caractéristique. Ainsi, pour réduire l’erreur d’interpolation, il y a trois solutions : 1. réduire le nombre d’éléments traversés par l’interface ; 2. modifier la valeur de la fonction caractéristique sur ces éléments ; 3. réduire le volume de ces éléments. Intuitivement, l’idée la plus simple à mettre en œuvre est de réduire le volume des éléments. Comme nous n’avons à priori aucune information sur la fa¸con de réduire ce volume, [BIGOT, 2001] a choisi de le réduire intelligemment selon une homothétie centrée par rapport au barycentre de l’élément.