Modélisation et simulation numérique du changement de phase liquide–vapeur en cavité

  • Published on

  • View

  • Download


  • C. R. Mecanique 334 (2006) 2533

    Modlisation et simulation numrique du changement de phaseliquidevapeur en cavit

    Virginie Daru a,b,, Marie-Christine Duluc a,c, Olivier Le Matre a,d,Damir Juric a, Patrick Le Qur a

    a LIMSI-CNRS, UPR 3251, B.P. 133, 91403 Orsay cedex, Franceb Laboratoire SINUMEF, ENSAM, 151, boulevard de lhpital, 75013 Paris, France

    c CNAM, 292, rue Saint-Martin, 75141 Paris cedex 03, Franced LME, universit dEvry, 40, rue du Pelvoux, CE1455, 91020 Evry cedex, France

    Reu le 30 avril 2005 ; accept aprs rvision le 18 octobre 2005Disponible sur Internet le 6 dcembre 2005

    Prsent par Sbastien Candel


    Un modle numrique dcoulement avec changement de phase liquidevapeur en cavit ferme est prsent. Linterface liquidevapeur est dcrite par une mthode de front tracking. Le liquide est considr comme parfaitement incompressible, tandis que lavapeur est assimile un gaz parfait faiblement compressible. Ce modle tient compte de la courbe de saturation. Des simulationssont ralises partir dun cas-test 1D, assimilable un autocuiseur. Les rsultats sont compars avec une solution approcheutilisant un modle faible Mach. Pour citer cet article : V. Daru et al., C. R. Mecanique 334 (2006). 2005 Acadmie des sciences. Publi par Elsevier SAS. Tous droits rservs.Abstract

    Modeling and numerical simulation of liquidvapor phase change in an enclosed cavity. A model for the simulation ofboiling flow with phase change in a closed cavity is presented. A front-tracking method is used to deal with the liquidvaporinterface. The liquid phase is incompressible while the vapor phase is weakly compressible and obeys to the perfect gas law.This model can deal with large density ratio (l/v 1000) flows while accounting for the saturation curve. Computations areperformed on a 1D validation case, idealizing a pressure cooker. Results are compared with a low Mach number approximation.To cite this article: V. Daru et al., C. R. Mecanique 334 (2006). 2005 Acadmie des sciences. Publi par Elsevier SAS. Tous droits rservs.

    Mots-cls : Mcanique des fluides numrique ; coulement compressible ; Changement de phase ; Liquidevapeur ; Courbe de saturationKeywords: Computational fluid mechanics; Compressible flow; Phase-change; Liquidvapor; Saturation curve

    * Auteur correspondant.Adresses e-mail : (V. Daru), (M.-C. Duluc), (O. Le Matre), (D. Juric), (P. Le Qur).1631-0721/$ see front matter 2005 Acadmie des sciences. Publi par Elsevier SAS. Tous droits rservs.doi:10.1016/j.crme.2005.10.015

  • 26 V. Daru et al. / C. R. Mecanique 334 (2006) 2533

    Abridged English version

    Numerical simulation of boiling flows is a real challenge from both modeling and computational standpoints. It hasreceived much attention in the recent past and is still the subject of intensive ongoing research. However, simulationsof liquidvapor flows with realistic physical properties and large density ratio between the phases (1000) are stillunavailable in the literature. This is due to various modeling and numerical difficulties. For instance, it is necessaryto consider variations of the vapor density with pressure and/or temperature, while in a closed cavity, any vaporproduction will lead to a significant pressure increase, which in turn affects the saturation temperature. Therefore, themodel must account for the saturation curve which links pressure and temperature at the interface. Also, the numericalmethod for such problems needs to deal with incompressible (liquid) and weakly compressible (vapor) flows in thesame computational domain. This difficulty can be overcome using different computational sub-domains for the twophases, but this approach is limited to small deformations of the interface. In contrast, the front-tracking techniquesbased on a one-field (one domain) formulation have shown to be effective in much complex situations in the limit ofboth incompressible phases (e.g. [1,2]). Thus, we present an extension of the one-field formulation of [4] to includephase-change processes. We report here early numerical tests for a 1D configuration similar to a pressure cooker.

    The model developed in the present paper uses a one-field formulation for the flow governed by the compressibleNavierStokes equations (1). The motion of the interface is described using a front-tracking method [1]: the Heavisidefunction indicates the liquid (H = 0) or the vapor phase (H = 1). In this formulation, the fluid density is written =l + H(x, t)(v l), where the liquid density l is constant, and the vapor density v obeys to the perfect gaslaw. On the interface, classical jump conditions are used, Eqs. (2), (3), the last term in Eq. (3) being neglected in thecomputations. Other variables of the problem and fluid properties are similarly written. We consider water with a latentheat decreasing linearly with the temperature and a 2nd order polynomial approximation of the Dupr law Eq. (4). Thetime-stepping procedure is based on an explicit integration of the governing equations, followed by a projection step toenforce the mass conservation: at each time-step, the pressure field is solution of Eq. (7), which is solved using a directsolver. For validation purposes, an approximate solution for the test case is developed, considering a motionless liquidphase bounding a compressible vapor phase with phase change at the fixed interface (jump conditions are boundaryconditions). Due to the very low vapor velocity involved in the test case, a low-Mach number approximation is usedin the vapor domain where the equations to be solved are given by Eqs. (8). In this formulation the pressure in thevapor phase is split into two components: the thermodynamic pressure P(t), solely function of time and the dynamicpressure pd(x, t), function of both time and space to enforce the mass conservation. The latter is, again, solution ofan elliptic equation.

    The test case configuration is similar to a pressure cooker: a 1D cavity with length L contains a liquidvapormixture initially under equilibrium conditions (P0, T0 = Tsat(P0)). The liquid phase is set in the left part of the cavity(0 x xs0) and the vapor in the right one (xs0 x L). The initial location of the interface is denoted xs0. Thetemperature of the cavity wall on the liquid side (x = 0) is suddenly increased to Tw > T0, while the opposite wall(x = L) is thermally insulated. The system then undergoes a vaporization process to reach a new uniform equilibriumstate (Tf = Tw , Pf = Psat(Tf ), xsf < xs0). Numerical simulations are performed for a cavity with length L =100 m. The physical properties of the fluid (water) used in the calculations are indicated in Section 3. Two differentinitial locations of the interface were examined. Analytical expressions derived from classical thermodynamics allowthe comparison between the theoretical and numerical equilibrium states. The interface location in the final steady-state xsf is given by Eq. (11) and the amount of heat Q transfered to the fluid during the vaporization process is givenby Eq. (12). The comparison is reported in Table 1 and shows an excellent agreement. As expected for this situation,the interface location is essentially unchanged. The fixed interface approximation in the low-Mach number model istherefore relevant.

    Before starting the discussion on transient phenomena, it is worthwhile to consider the energy equation written inthe dimensionless form, Eq. (14). It shows that the problem is governed by a single dimensionless number, Dilv , whichcompares the two time-scales for the heat diffusion in the liquid and vapor domains. For xs0 = 0.9L, Dilv 1 andthe vapor dynamics is dominated by the heat diffusion. In contrast, when xs0 = 0.1L, Dilv 1 and advection, com-pressibility and phase change effects are expected to become significant too, yielding much more complex dynamics.Transient results for xs0 = 0.9L are presented in Figs. 14. Time evolutions of the temperature fields are plotted inFig. 1 for the two models. It is first shown that the two models (one-field formulation and low-Mach approximation)

    are in close agreement. The same agreement is observed for the thermodynamic pressure as a function of time re-

  • V. Daru et al. / C. R. Mecanique 334 (2006) 2533 27

    ported in Fig. 2. For the one-field model, plotted is the spatially averaged pressure, which however is nearly uniformin space (not shown). It is also reported that the pressure rise is not immediate but starts when the heat flux at theinterface becomes significant with a non negligible vaporization rate as a result. The interfacial velocity computed forthe one-field formulation is plotted in Fig. 3. The interfacial velocity presents a sharp acceleration at the early stage,followed by a longer relaxation period which ends when the equilibrium state is achieved. The maximal interfacialvelocity is roughly 107 m/s. The velocity fields at different times reported in Fig. 4 are exactly equal to zero in theliquid phase, and are almost linear function of x in the vapor phase with a maximum at the interface (and vanish onthe right-hand wall).

    For this test case, Fig. 1 shows that the vapor temperature evolves monotonically in time to its equilibrium value andis nearly uniform in space as expected from the dimensional analysis which indicates that the dynamics is dominatedby the diffusion of heat: interfacial and vapor temperatures are therefore identical at any time. A different behavioris reported for xs0 = 0.1L. For this case, the liquid layer separating the heated wall and the interface is thinnerleading to higher heat fluxes and a higher vaporization rate at the early stage of the process. As a consequence, themagnitude of the pressure increase is larger and heats the vapor through compressibility effects. Although, on theinterface the evolutions of the temperature and pressure are not independent but constrained by the saturation curve.Differentiating this dependence with time gives dTsat/dt = (dTsat/dP) (dP/dt). For water, the pressure-inducedtemperature increase is larger in the vapor (compressibility effect) than at the interface (saturation curve), giving riseto a temperature gradient. Then, if the heat diffusion is not fast enough to homogenize the vapor temperature thetemperature gradient is sustained and in turn enhances the vaporization rate, favoring the compressibility effects. Theconsequence of this dynamics is clearly visible in the curves given in Fig. 5, corresponding to xs0 = 0.1L, where themaxima of the vapor temperature are achieved at the insulated wall while the minima stand at the interface. It hasto be underlined that the maximal vapor temperature even temporarily exceeds Tw . Also, Fig. 5 exhibits small butnoticeable discrepancies between the one-field formulation and the low-Mach approximation. Since for this test casethe maximal Mach number is less than 103 and the interface motion negligible, the differences in the two modelspredictions are likely due to the front-tracking treatment. This claim is verified in Fig. 6 where an improvement inthe agreement between the two models is achieved by using a refined mesh (500 mesh points instead of 200) in theone-field model. This observation calls for an improvement of the front-tracking algorithm.

    To conclude, this Note has presented preliminary computations of 1D boiling flows. From the modeling point ofview the main focus was on the satisfaction of the saturation curve and realistic density radio. It is shown that thecompressible one-field formulation, coupled with a front-tracking algorithm, efficiently predicts equilibrium states.Comparisons with a fixed interface low-Mach approximation have shown that transient behaviors are also correctlyaccounted for by the model. However, the model requires very refined grids when compressibility effects are signifi-cant compared to heat diffusion in the vapor, i.e. when Dilv 1. For such situations, the description of the interfacehas to be improved. Current efforts include the extension of the model to 2D configurations and surface tension effects.

    1. Introduction

    La simulation numrique des coulements liquidevapeur dans des conditions ralistes est un vritable dfi eugard aux multiples difficults rencontres, relatives aussi bien la modlisation physique quaux mthodes num-riques. En effet, il sagit de reprsenter un coulement comportant une interface liquidevapeur de dynamique rapideavec changement de phase. Du point de vue numrique, le problme du suivi dune interface travers laquelle lesproprits du fluide sont discontinues a fait lobjet dun grand nombre de travaux. Parmi les mthodes sur grille fixe,la mthode de front-tracking que nous utilisons, fonde sur un suivi explicite de linterface reprsente par une suitede marqueurs, a fait la preuve dune remarquable prcision pour des simulations de diffrents coulements multipha-siques [1]. En particulier, elle a t rcemment utilise pour effectuer des simulations dcoulements multiphasiquesincompressibles tridimensionnels avec changement de phase [2]. Cependant, le modle incompressible nest pas ap-propri pour dcrire correctement, par exemple, une configuration o les phases liquide et vapeur coexistent dansune enceinte ferme ou encore une configuration dans laquelle la dtermination de la pression thermodynamique estncessaire (exemple : un autocuiseur avant ouverture de la soupape).

    On prsente ici un modle destin dcrire ces coulements avec changement de phase liquidevapeur lintrieurdune enceinte ferme. Le modle retenu est un modle un fluide dans lequel lcoulement est rgi par un jeu

    dquations unique pour les deux phases selon une approche comparable celle dveloppe par Juric et Tryggvason

  • 28 V. Daru et al. / C. R. Mecanique 334 (2006) 2533

    en incompressible [3]. La phase vapeur est compressible tandis que la phase liquide reste parfaitement incompressible,le suivi dinterface tant effectu par front tracking. La formulation retenue fait apparatre une seule pression, qui estla pression totale. Dans un travail prcdent, ce modle a t test et valid sur un coulement liquidegaz 2D sanschangement de phase lintrieur dune cavit chauffe [4].

    Les premiers rsultats obtenus dans un cas 1D pour un coulement liquidevapeur avec changement de phasesont compars pour les valeurs stationnaires aux relations classiques issues de la thermodynamique. Dans le but devalider les simulations transitoires, un deuxime modle a t mis en uvre, utilisant une rsolution spare dans leliquide et dans la vapeur. Dans la phase vapeur, ce modle utilise une approximation de type faible Mach . Unecomparaison est ralise entre les deux approches. Les premiers rsultats mettent notamment en vidence leffet de lacompressibilit sur la vaporisation linterface.

    2. Modles mathmatiques et rsolution numrique

    2.1. Modle pression unique

    Lcoulement obit aux quations de NavierStokes sous la forme compressible gnrale (les effets de dissipationvisqueuse, de gravit et de tension superficielle sont ngligs dans cette premire approche) :

    t+ div(v) = 0


    t+ div(v v) + p = div( )


    t+ cvv T = div(kT ) T p


    div(v) + mLv H(1)

    compltes par la loi dtat gnrale l = H(x, t)(v l). Dans ce travail, la vapeur est assimile un gaz parfait,v = p/(rT ) et le liquide, incompressible, est de masse volumique constante l . H(x, t) est une fonction de Heavisidereprant le liquide ou le gaz, telle que H(x, t) = 1 dans le gaz et H(x, t) = 0 dans le liquide. La mthode de suivilagrangien de linterface liquidegaz nous permet de construire une telle fonction sur le maillage eulrien discrtisantlensemble du domaine fluide [2]. Dans lquation de lnergie, le dernier terme modlise le changement de phase ;m est la densit de flux de masse interfacial et Lv la chaleur latente, donne par : Lv = Lv0 + (cpv cpl )(T T0),avec Lv0 , la chaleur latente la temprature de saturation T0.

    Ce modle est complt par les relations de saut linterface (suppose de masse nulle), qui permettent de calculersa vitesse. Si lon fait lhypothse que linterface est fine et sans masse, elles scrivent [5] :

    l(vl s) n = v(vv s) n = m (2)

    (qv q l ) n = mLv m3





    avec s la vitesse de linterface, n la normale linterface et q la densit de flux de chaleur. La relation de saut pourla quantit de mouvement nest pas explicite ici car non directement utilise dans la rsolution. Dans les simulationsnous ngligerons le dernier terme dans le membre de droite de (3).

    Enfin la fermeture du modle ncessite la donne de la temprature linterface, que nous imposons gale latemprature de saturation Tsat(p), rgie par la loi de Dupr. Pour le fluide (eau) et la plage de temprature tudis,cette loi peut tre lisse par un polynme du second degr :

    Tsat(p) = 5,92 1010p2 + 3,862 104p + 340,18 (4)La rsolution numrique de ces quations repose sur une mthode de type projection, classique pour les coulementsincompressibles. On ne prsente ici que la discrtisation temporelle, dans la mesure o la discrtisation spatiale utiliseclassiquement des formules aux diffrences finies centres. Pour ces premiers rsultats, nous avons employ une

    discrtisation spatiale centre dordre 2 et une discrtisation explicite dordre 1 en temps. Par ailleurs la discrtisation

  • V. Daru et al. / C. R. Mecanique 334 (2006) 2533 29

    de lquation de lnergie ne posant pas de problme particulier, on ne la prsentera pas. Lquation de continuit estdiscrtise de faon implicite, soit :




    rT n+1 l

    ) Hn


    rT n l

    ))+ div(v)n+1 = 0 (5)

    Lquation de quantit de mouvement est discrtise comme :

    (v)n+1 (v)nt

    = pn+1 + div( v v)n = pn+1 + An (6)

    En prenant la divergence de (6) et en combinant avec (5), on obtient lquation pour la pression :

    2pn+1 1t2



    rT n+1 l

    ) Hn


    rT n l

    ))= 1

    tdiv(v)n + div(A)n. (7)

    Cette quation dHelmholtz pour pn+1 est rsolue laide dun solveur direct. Lalgorithme gnral de rsolutionpour une itration en temps est ainsi le suivant :

    calcul de Tsat(p) par (4) ; calcul de la densit de flux de chaleur de part et dautre de linterface et du terme de flux de masse m H ; calcul de Hn+1 laide de la mthode de suivi dinterface ; calcul de la temprature ; calcul de la pression par (7) et de la masse volumique par la loi dtat ; calcul de la vitesse par (6).

    Bien que lacoustique soit prise en compte dans le gaz, la rsolution implicite de la pression permet de saffranchirde la condition de stabilit explicite lie la vitesse du son, qui serait extrmement restrictive sur le pas de tempsdintgration. De plus, la diffusion numrique induite par limplicitation a pour effet bnfique de filtrer trs rapide-ment les ondes acoustiques sans intrt pour cette tude. Enfin, le liquide peut tre trait comme un fluide strictementincompressible grce cette implicitation.

    2.2. Modle faible Mach

    Ce modle dissocie le traitement des deux phases dont linterface est fixe, ce qui est lgitime du fait des faiblesvariations de volume pour les cas traits dans la suite. Il y a donc deux domaines avec des conditions de raccord portantsur la temprature et les flux. On limine ainsi toute source derreur lie ltalement artificiel de linterface induit parle front-tracking. Dans le domaine liquide une quation de diffusion sur T est employe avec comme condition surlinterface Tl = Tsat(P ). Dans le domaine vapeur, le nombre de Mach est petit et une approximation de type faibleMach (voir par exemple [6,7]) est employe. Dans cette approximation la pression est spare en deux composantes,la pression thermodynamique P(t) uniforme en espace et la pression dynamique pd(x, t) ; la seconde tant dordreM2 relativement la premire. Les quations dans la phase vapeur scrivent ainsi :


    t+ div(vvv) = 0


    t+ div(vvv vv) + pd = div( v)


    t+ vcpvvv Tv = div(kvTv) +



    compltes de la loi dtat P = rvTv . En drivant en temps lquation dtat et en combinant avec les quationsprcdentes, il est possible dobtenir une nouvelle expression de la conservation de la masse :[ ]v

    t= 1

    Tv 1

    rdiv(kvTv) + vvv Tv + 1




  • 30 V. Daru et al. / C. R. Mecanique 334 (2006) 2533

    Tableau 1Comparaison des tats dquilibre thoriques et calculs (modle pression unique)Table 1Comparison of the theoretical and computed (one-field formulation) equilibrium statesxs0 (m) xsf (m) Q (J/m2)

    Thorie Sim. Num. Thorie Sim. Num.10 9,953708 9,953699 895,7 894,750 49,97430 49,97426 4071,5 4070,890 89,99486 89,99482 7247,3 7246,9

    Par intgration de cette dernire quation sur le volume de vapeur v , on obtient une expression reliant le tauxinstantan de variation de P celui de la masse de vapeur m :

    m =v


    [ 1

    rdiv(kvTv) + vvv Tv

    ]dv + 1




    dv (10)

    Les conditions aux limites sur linterface sont Tv = Tsat(P ) et les relations approches de saut : m = vvv n etmLv(T ) = (q l qv) n. Connaissant la solution au temps tn, le flux m est calcul par la relation de saut. m estalors utilis pour calculer dP/dt par Eq. (10). On procde ensuite une intgration explicite jusque tn+1 = tn + t(schma dAdamsBashforth dordre 2 en temps) de la pression P et de v selon Eq. (9). Tv est actualise dans ven utilisant la relation dtat, sauf sur linterface o la condition Tv = Tsat(P ) est employe. Lquation de quantitde mouvement est alors intgre, par une technique de prdiction-projection : pd est dtermine par une quationelliptique pour satisfaire div(vvv) = t, la seconde relation approche de saut vvv n donnant la condition auxlimites. Pour finir, la temprature dans le liquide est calcule, avec linterface Tl donne par la nouvelle pression Pau temps tn+1.

    3. Rsultats et discussion

    On prsente les rsultats obtenus dans le cas 1D cartsien suivant : une cavit linique de longueur L contientune phase liquide (0 x xs ) en quilibre avec sa vapeur (xs x L) ; la position de linterface est repre parla coordonne xs . Le systme est initialement la pression P0, la temprature de saturation T0 = Tsat(P0) et laposition de linterface est xs0. La paroi en x = L est adiabatique tandis la paroi x = 0 est brusquement porte latemprature Tw . Ltat stationnaire obtenu correspond une temprature et une pression uniformes dans la bote,respectivement gales Tf = Tw et la pression de vapeur saturante Pf = Psat(Tw). Les calculs sont conduits dansune cavit de longueur L = 100 m. Le fluide est de leau avec les proprits thermophysiques suivantes : P0 =101 325 Pa, Lv(T0) = 2 251 200 J kg1, l = 958,8 kg/m3, kl = 0,68 W m1 K1, cpl = cvl = cl = 4216 J kg1 K1,l = 2,79104 kg m1 s1, r = 461,89 J kg1 K1, kv = 0,0248 W m1 K1, cpv = 2034 J kg1 K1, v = 1,21105 kg m1 s1. La temprature impose en paroi vaut Tw = 393,15 K. Certains rsultats peuvent tre comparsavec des valeurs thoriques issues de la thermodynamique classique. La position de linterface dans ltat final estobtenue partir de la relation suivante = f vf + (1 f )l , o f est le taux de vide dans ltat final dfini parf = (L xsf )/L. La masse volumique reste constante tout au long de la transformation, ce qui permet dcrire :

    xsf = L l v0l vf (L xs0) (11)

    La quantit de chaleur consomme au cours de la transformation, Q, est calcule en utilisant le premier principeappliqu une volution isochore : Q = Uf U0, o U dsigne lnergie interne du fluide contenu dans la cavit. Ladtermination de Q, indpendante du chemin suivi, est ralise en considrant dans un premier temps la condensationcomplte de la vapeur sur le palier de liqufaction P = P0, suivie dun chauffage isochore du liquide saturant jusqula temprature de saturation Tf et enfin de la vaporisation sur le palier de liqufaction P = Pf dune masse de vapeurgale m . La quantit m reprsente la masse de vapeur prsente ltat final, gale m = (L x ). Onvf vf vf vf sfobtient finalement pour Q :

  • V. Daru et al. / C. R. Mecanique 334 (2006) 2533 31

    Q = (L xs0)[v0F(T0) vf l v0

    l vf F (Tf )]

    + [(lxs0 + v(L xs0))cl(Tf T0)] (12)o

    F(Tk) = pk(



    ) Lv(Tk), k = 0, f (13)

    Les rsultats obtenus pour trois valeurs de la position initiale de linterface, xs0 = 0,1L, 0,5 L et 0,9L, sont prsen-ts dans le Tableau 1. Laccord entre les simulations et les valeurs thoriques prdites par la thermodynamique esttrs satisfaisant. Dans les deux cas, le dplacement de linterface est infime, justifiant ainsi lhypothse dinterfaceimmobile retenue pour le modle faible Mach. Avant dexaminer lvolution instationnaire du systme, commenonspar mettre en vidence les paramtres sans dimension caractrisant son comportement dans la phase vapeur. Un adi-mensionnement convenable est obtenu en choisissant pour temps de rfrence le temps caractristique de diffusionthermique dans le liquide, soit tref = L2liq/l , avec l = kl/(lcl) et Lliq la longueur de liquide. En effet, cest bienla diffusion dans le liquide qui est lorigine de la vaporisation sur linterface et par suite de laugmentation de lapression. La vitesse dans la vapeur est une vitesse de transport convectif. Nous poserons donc vref = (Lvap)/tref avecLvap la longueur de vapeur, soit encore vref = l(Lvap)/L2liq. Les autres grandeurs de rfrence sont fixes aux valeursinitiales (P0, T0) exception faite de la chaleur latente dont la valeur de rfrence est prise gale cvvT0. Lquation delnergie sans dimension scrit ainsi dans la phase vapeur :


    + v T)

    = Dilv div(T ) ( 1)p div(v) + mLv H (14)avec

    Dilv =(







    v0, v0 = P0


    Dilv , seul paramtre sans dimension caractrisant le problme, reprsente le rapport des temps caractristiques dediffusion thermique respectivement dans la lame liquide et dans la lame de vapeur. Ainsi, quelle que soit la positioninitiale de linterface, les termes dadvection, de compressibilit et de changement de phase sont du mme ordre degrandeur. Le facteur Dilv devant le terme div(T ) est gal 1,29 104 pour xs0 = 0,9L, et 1,97 pour xs0 = 0,1L.La diffusion thermique est donc leffet dominant dans la vapeur lorsque xs0 = 0,9L, tandis quelle est du mme ordrede grandeur que les autres termes lorsque xs0 = 0,1L.

    On prsente Figs. 1 4, lvolution instationnaire pour le cas xs0 = 0,9L. Les rsultats obtenus sur le champ detemprature avec les deux modles sont superposs sur la Fig. 1 et montrent un excellent accord. Il en est de mme pourla pression thermodynamique (Fig. 2). Dans le cas du modle pression unique, il sagit dune pression moyenne enespace (toutefois la pression est quasiment uniforme). La pression commence crotre lorsque la vaporisation dmarre,soit lorsque le flux de chaleur sur linterface devient significatif. La vitesse de linterface, de lordre de 107 m/s,montre une phase dacclration suivie dune phase de dclration plus lente (Fig. 3). La vitesse, rigoureusementnulle dans le liquide, est chaque instant pratiquement linaire dans la vapeur (Fig. 4). Dans ce cas (xs0 = 0,9L),le profil de temprature prsente une volution monotone (Fig. 1), la temprature tant quasiment uniforme dansla vapeur. Ceci est une consquence du fait que, comme nous lavons montr prcdemment, la convection et lacompressibilit ont peu deffet sur la dynamique dans la vapeur, essentiellement domine par la diffusion thermique.Comme par ailleurs sur linterface la temprature est fixe gale Tsat(p), le profil de temprature dans la vapeur estun profil de conduction thermique pure avec une condition aux limites de Dirichlet gauche et de Neumann droite,donc constant en espace et gal Tsat(p).

    Dans le deuxime cas pour lequel xs0 = 0,1L, prsent dans la Fig. 5 lpaisseur du domaine liquide plus faibleconduit des phnomnes diffrents du cas prcdent. Le flux de chaleur sur linterface, plus lev dans ce cas,engendre un flux de vaporisation bien suprieur. Laugmentation de pression qui en rsulte conduit un chauffageplus rapide dans la vapeur. La temprature de linterface, contrainte par la courbe de saturation, crot en revanchemoins vite. Il stablit ainsi un gradient de temprature dans la vapeur, le maximum de temprature, obtenu en x = L ;la temprature de vapeur pouvant mme temporairement dpasser la temprature de chauffage Tw . En consquence,les flux de chaleur linterface provenant du liquide et de la vapeur renforcent le taux de vaporisation du liquide.

    Cet exemple illustre le rle favorable des effets de compressibilit sur la vaporisation. La comparaison des rsultats

  • 32 V. Daru et al. / C. R. Mecanique 334 (2006) 2533

    Fig. 1. Champs de temprature diffrents instants : comparaisondes deux modles et xs0 = 0,9L (intervalle de temps 10 ms entrechaque courbe).Fig. 1. Temperature fields at different times for the two models andxs0 = 0.9L (time interval between curve is 10 ms).

    Fig. 2. Evolution de la pression thermodynamique en fonction dutemps pour les deux modles et xs0 = 0,9L.Fig. 2. Time evolutions of the thermodynamic pressure for the twomodels and xs0 = 0.9L.

    Fig. 3. Evolution temporelle de la vitesse de linterface pourxs0 = 0,9L (modle pression unique).Fig. 3. Time evolution of the interfacial velocity for xs0 = 0.9L(one-field formulation).

    Fig. 4. Champs de vitesse diffrents instants pour xs0 = 0,9L(modle pression unique).Fig. 4. Velocity fields at different times for xs0 = 0.9L (one-fieldformulation).

    obtenus par les deux modles fait apparatre de lgres diffrences sur les champs de temprature. Un raffinement dumaillage pour le modle pression unique (Fig. 6) permet dattnuer ces diffrences, notamment prs de linterface.Lcart rsiduel dcoule de la reprsentation de linterface en front tracking sur une zone de largeur contrle dequelques mailles. Cette zone est ainsi une zone de raccord des modles correspondant chacune des phases (i.e.pour H = 1 vapeur ou H = 0 liquide). Ce raccord, actuellement effectu par simple interpolation, devra fairelobjet damliorations. Enfin, puisque la configuration xs0 = 0,1L est celle qui engendre les plus grandes valeurs devitesse dans la vapeur, le nombre de Mach le plus lev peut tre estim partir des rsultats de cette simulation. Savaleur nexcde jamais 103. Il est ainsi pertinent de considrer la vapeur comme faiblement compressible et a fortioridutiliser un modle faible Mach pour la dcrire.

    4. Conclusions

    Le modle prsent ici, destin la simulation numrique dcoulements avec changement de phase liquidevapeura pour double objectif de rendre compte dcoulements rels grand rapport de densit en cavit ferme et de satisfairela courbe de saturation. Les rsultats obtenus sur un cas-test 1D cartsien assimilable un autocuiseur, attestent de la

    capacit du modle atteindre ce double objectif, du moins dans des conditions o la dynamique de linterface reste

  • V. Daru et al. / C. R. Mecanique 334 (2006) 2533 33Fig. 5. Champs de temprature pour xs0 = 0,1L : modle pressionunique (maillage 200 points) et approximation faible Mach.Fig. 5. Temperature field for xs0 = 0.1L: one-field formulation(pression unique) using 200 mesh points and low-Mach approxi-mation (faible Mach).

    Fig. 6. Champs de temprature pour xs0 = 0,1L : modle pressionunique (maillage 500 points) et approximation faible Mach.Fig. 6. Temperature field for xs0 = 0.1L: one-field formulation(pression unique) using 500 mesh points and low-Mach approxi-mation (faible Mach).

    lente. La comparaison avec une solution approche, obtenue partir dun code utilisant une approche faible Machmontre un excellent accord et rvle que le traitement de linterface doit tre amlior. Les rsultats obtenus dansle cas xs0 = 0,1L montrent que la description spatiale de la phase vapeur est indispensable pour rendre compte delintgralit des phnomnes. Le travail en cours porte sur limplantation de la tension superficielle (1D sphrique)ainsi que sur lextension un code 2D. Des simulations avec une dynamique dinterface plus rapide vont galementtre dveloppes.


    [1] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the compu-tation of multiphase flows, J. Comput. Phys. 169 (2001) 708759.

    [2] S. Shin, D. Juric, Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connec-tivity, J. Comput. Phys. 180 (2002) 427470.

    [3] D. Juric, G. Tryggvason, Computations of boiling flows, Int. J. Multiphase Flow 24 (1998) 387410.[4] D. Juric, V. Daru, M.C. Duluc, Simulation dcoulements liquidegaz lintrieur dune cavit chauffe, in : Congrs Franais de Thermique,

    SFT, 2004, pp. 95100.[5] J.M. Delhaye, Jump conditions and entropy sources in two-phase systems. Local instant formulation, Int. J. Multiphase Flow 1 (1974) 395409.[6] D.R. Chenoweth, S. Paolucci, Natural convection in an enclosed vertical layer with large horizontal temperature differences, J. Fluid Mech. 169

    (1986) 173210.[7] A. Majda, J. Sethian, The derivation and numerical solution of the equations for zero Mach number combustion, Combust. Sci. Technol. 42

    (1985) 185205.


View more >