Télécharger difasyk2Dax.dgibi

Retour à la liste

* fichier :  difasyk2Dax.dgibi
opti trac ps ;
opti dime 2 elem qua8 mode axis ;
*'OPTION' 'echo' 0 ;
*
*
**************************************************************************
*
*        CAS TEST DIPHASIQUE LIQUIDE-LIQUIDE 2D PLAN
*
***************************************************************************
*
*
*        Présentation: Ce cas test permet la simulation de la sédimention
*                      de gouttes sphériques de liquide immergé dans
*                      un milieu continu (autre liquide) de densité moindre.
*                      Bien entendu les densités peuvent être modifiées.
*                      les deux liquides sont dans un extracteur
*                      centrifuges et sont soumis à une
*                      force centrifuge horizontale
*
*
*        Le modèle : Le modèle simulé consiste à résoudre 2 équations de la
*                    quantité de mouvement sur les flux volumiques des 2
*                    phases JC et JD. Ensuite on résout une équation de la
*                    conservation de la masse qui porte sur la fraction
*                    volumique de la phase dispersée : AD (pour alpha_D).
*                    La conservation de la fraction voumique de la phase
*                    continue AC s'en déduite par AD + AC = 1.
*                    L'équation d'incompressibilité du fluide est
*                    Div(JC+JD) = 0.
*
*                    Des équations de mouvements sont rajoutées pour
*                    les composantes orthoradiales pour la prise en
*                    compte des forces de Coriolis.
*
*                    le modèle de coalescence est simplifié, ce qui
*                    revient à jouer sur la viscosité efficace de l'émulsion
*                    MUEFF qui tend vers la valeur de la phase pure
*                    reconstituée lorsque la fraction volumique dépasse
*                    un certain seuil. 
*                    Le front de sédimentation est rès raide
*                    notre problème revient regarder évoluer le front
*                    de sédimentation de gouttes en sédimentation dans
*                    un extracteur perpétuellement alimenté en mélange
*                    et pourvu de deux sorties expulsant des phases pures
*                    séparées.
*
*
*         Auteur     : Gilles Bernard-Michel
*
*         Révision   : version 1
*
*         date       : 21/12/99
*
*         REMARQUES :
*
*  Nous somme en diphasique et en mode axis
*
*  différence entre mueff et hinder pas à jour pour exposant.
*
* AD : ALPD au Sommet 
* AC : ALPC au Sommet 
*
* implicitation en alpha pour les frottements.
*
* attention meme mueff pour les dex phases ce qui est faux.
* on doit tendre vers les bons mu homogenes.
*
*
************************ Options générales *********************************
*
*
** Sortie graphique.
*
 
GRAPH = 'O';
 
*
*test court uniquement.
*
 
COMPLET = FAUX;
 
*
*période des sorties graphiques : tous les NITMA/PERIO pas de temps
*
 
PERIO = 1 ;
 
*
* nombre d'itérations en temps
*
 
NITMA = 1;
 
*
* nombre de points fixes max.
*
 
NITER = 2;
 
*
* Précision du point fixe
*
 
PREC = 1.D-4;
 
* 
* Pas de temps
*
 
DELTATPS = 0.2 ;
 
*
********************* On traite les aspects géométriques ********************
*
** viscosité dynamique mu/rho (utilisé uniquement pour le maillage).
*
 
nu =1.E-6 ;
 
*
** La vitesse de rotation initiale. (rad/s)
*
 
omeg = 100 ;
 
*
** Les paramètres de la géométrie.
*
 
rint     = 0.1 ;
rext     = 0.2;
hh       = 0.044 ;
inj1     = 0.145 ;
inj2     = 0.155;
exi1     = 0.102;
exi2     = 0.198;
 
*
** Coefficients affectés aux couches d'Ekman et Stewartson.
*
 
 coefekma = 8. ;
 coefste1 = 4. ;
 coefste2 = 2. ;
 
*
** Densité du maillage.
*
 
npts = 12. ;
d1 = 0.1 / npts;
d2x1 = d1 / 10. ;
d2x2 = d1 / 1.4 ;
d2y = d1 /  2. ;
d2y2 = d2y / 10. ;
 
*
** Calcul des épaisseurs des couches d'Ekman et de Stewartson.
** Ces épaisseurs sont des ordres de grandeurs sans dimension.
** on épaissit un peu les mailles en prenant nu = 10 *nu
 
nu = 10.*nu;
ekma = nu / (omeg * (rext**2));
ddek = ekma**0.5;
stew = ekma**0.25;
 
*
** Modifications de ces épaisseurs en tenant compte des coefficents
** introduits.
*
 
dbek = coefekma * ddek * hh ;
fiek = hh * (1. - (coefekma * ddek)) ;
dste = rint * coefste1 * stew ;
fste = rext * coefste2 * stew ;
 
*
** Les sommets du bol
*
 
p1 = rint 0. ;
p2 = rext 0. ;
p3 = rext hh ;
p4 = rint hh ;
 
*
** Les milieux des cotés.
*
 
mil1 = ((rint + rext)/2.0);
mil2 = hh/2.0 ;
p6 = rext mil2 ;
p8 = rint mil2 ;
p68 = inj1 mil2 ;
*
** Les points délimitant les couches d'Ekman.
*
 
p11 = rint dbek ;
p22 = rext dbek ;
p33 = rext fiek ;
p44 = rint fiek ;
 
*
** Les points d'entrée du fluide et leurs symétriques.
** Sont rajoutées les couches de Stewartson qui les 
** entourent.
*
 
ent1 = inj1 0. ;
ent2 = inj2 0. ;
sym1 = inj1 hh ;
sym2 = inj2 hh ;
en11 = (inj1 - fste) 0. ;
en22 = (inj2 + fste) 0. ;
sy11 = (inj1 - fste) hh ;
sy22 = (inj2 + fste) hh ;
 
*
ents = inj1 dbek ;
syms = inj1 fiek ;
*
 
*
** Les points de sortie du fluide et leurs symétriques.
** Sont rajoutées les couches de Stewartson qui les 
** entourent.
*
 
out1 = exi1 hh ;
out2 = exi2 hh ;
syo1 = exi1 0. ;
syo2 = exi2 0. ;
ou11 = (exi1 '+' fste) hh;
ou22 = (exi2 '-' fste) hh;
so11 = (exi1 '+' fste) 0.D0;
so22 = (exi2 '-' fste) 0.D0;
 
 
*
** Assemblage des bords du domaine.
*
 
bas1 = p1 'DROIT' dini d2x1 dfin d2x1 syo1
          'DROIT' dini d2x1 dfin d1 so11
          'DROIT' dini d1 dfin d1 en11
          'DROIT' dini d1 dfin d2x2 ent1;
entr = ent1 'DROIT' dini d2x2 dfin d2x2 ent2 ;
bas2 = ent2 'DROIT' dini d2x2 dfin d1 en22
            'DROIT' dini d1 dfin d1 so22
            'DROIT' dini d1 dfin d2x1 syo2
            'DROIT' dini d2x1 dfin d2x1 p2; 
cdro  = p2 'DROIT' dini d2y dfin d2y p22
           'DROIT' dini d2y dfin d1 p6 
           'DROIT' dini d1 dfin d2y p33
           'DROIT' dini d2y dfin d2y p3;
sor2  = p3 'DROIT' dini d2x1 dfin d2x1 out2;
hau   = out2 'DROIT' dini d2x1 dfin d1 ou22
          'DROIT' dini d1 dfin d1 sy22
          'DROIT' dini d1 dfin d2x2 sym2 
          'DROIT' dini d2x2 dfin d2x2 sym1
          'DROIT' dini d2x2 dfin d1 sy11 
          'DROIT' dini d1 dfin d1 ou11
          'DROIT' dini d1 dfin d2x1 out1;
sor1  = out1 'DROIT' dini d2x1 dfin d2x1 p4; 
cgau = p4 'DROIT' dini d2y dfin d2y p44
           'DROIT' dini d2y dfin d1 p8
           'DROIT' dini d1 dfin d2y p11
           'DROIT' dini d2y dfin d2y p1; 
cot1 = bas1 'ET' entr 'ET' bas2;
cot2 = cdro;
cot3 = sor2 'ET' hau 'ET' sor1;
cot4 = cgau;
cnt  = cot1 'ET' cot2 'ET' cot3 'ET' cot4 ;
*
*
*
*borbg =  p1 'DROIT' dini d2x1 dfin d2x1 syo1
*            'DROIT' dini d2x1 dfin d1 ent1;
*borbd = ent1 'DROIT' dini d1 dfin d2x1 syo2
*            'DROIT' dini d2x1 dfin d2x1 p2;
*           
*bord  =  p2 'DROIT' dini d2y2 dfin d2y2 p22
*           'DROIT' dini d1 dfin d1 p6 
*           'DROIT' dini d1 dfin d1 p33
*           'DROIT' dini d2y2 dfin d2y2 p3;
 
*borhd =  p3 'DROIT' dini d2x1 dfin d2x1 out2
*            'DROIT' dini d2x1 dfin d1 sym1;
*borhg = sym1 'DROIT' dini d1 dfin d2x1 out1
*            'DROIT' dini d2x1 dfin d2x1 p4;
*           
*borg =  p4 'DROIT' dini d2y2 dfin d2y2 p44
*           'DROIT' dini d1 dfin d1 p8
*           'DROIT' dini d1 dfin d1 p11
*           'DROIT' dini d2y2 dfin d2y2 p1; 
 
*milieud = ent1  'DROIT' dini d2y2 dfin d2y2 ents
*           'DROIT' dini d1 dfin d1 p68 
*           'DROIT' dini d1 dfin d1 syms
*           'DROIT' dini d2y2 dfin d2y2 sym1;
*milieug = 'INVERSE' milieud ;
*
** On maille avec des rectangles.
*
*mt1 = borbg  milieud borhg  borg  'DALLER' ;
*mt2 = borbd  bord borhd milieug   'DALLER' ;
*mt = mt1 'ET' mt2 ; 
mt= cot1 cot2 cot3 cot4 'DALLER' ;
'TASSER' mt ;
*'TITRE' 'Maillage de la centrifugeuse';
*'TRACER' mt ;
 
*
** On réoriente les éléments.
*
 
mt = 'ORIENTER' mt ;
 
 
* définition du modèle.
 
&mt    = CHANGER mt   QUAF ;
*&bord  = CHANGER bord QUAF ;
*&borg  = CHANGER borg QUAF ;
*&borhd  = CHANGER borhd QUAF ;
*&borhg  = CHANGER borhg QUAF ;
*&borbd  = CHANGER borbd QUAF ;
*&borbg  = CHANGER borbg QUAF ;
*&milieud = CHANGER milieud QUAF ;
 
&sor1  = CHANGER sor1 QUAF ;
&sor2  = CHANGER sor2 QUAF ;
&bas1  = CHANGER bas1  QUAF ;
&bas2  = CHANGER bas2  QUAF ;
&cdro  = 'CHANGER' cdro QUAF ;
&cgau  = 'CHANGER' cgau QUAF ;
&entr  = 'CHANGER' entr QUAF ;
&hau  = CHANGER hau QUAF ;
&cot1 = CHANGER cot1 QUAF;
&cot2 = CHANGER cot2 QUAF;
&cot3 = CHANGER cot3 QUAF;
&cot4 = CHANGER cot4 QUAF;
&cnt  = CHANGER cnt  QUAF ;
 
'ELIMINATION' 1.e-6 (&mt  'ET' &bas1 'ET' &bas2 'ET' &cdro 'ET' &cgau
'ET' &entr 'ET' &hau 'ET' &cot1 'ET' &cot2
'ET' &cot3 'ET' &cot4 'ET' &cnt 'ET' &sor1 'ET' &sor2);
 
*'ELIMINATION' 1.e-6 (&mt 
*'ET' &borbd 'ET' &borbg 'ET' &bord 'ET' &borhd 'ET' &borhg 'ET' &borg
*'ET' &milieud ) ;
 
$mt    = MODE &mt    'NAVIER_STOKES' QUAF ;
$sor1  = MODE &sor1    'NAVIER_STOKES' QUAF ;
$sor2  = MODE &sor2    'NAVIER_STOKES' QUAF ;
 
*
** On crée le domaine associé au maillage.
*
 
cnt = 'CONTOUR' mt;
 
mt   = doma $mt    maillage ;
 
lmt = 'CHANGE' mt LIGNE ;
coup = lmt 'POINT' 'DROITE' p8 p6 1.D-5 ;
alpdco = lmt 'ELEM' 'APPUYE' 'STRICTEMENT' coup ;
 
 
 
*
** on crée une frontière moins un noeud pour le cas
** de cavités carrés. Il ne faut pas retirer un coin
** car on peut avoir un champ sortant.
*
 
*$cnt =  'MODELISER' &cnt 'NAVIER_STOKES' QUAF;
*cntcc = 'DOMA' $cnt maillage ;
*cntcc = 'CHANGER' cntcc POI1;
 
*retir = p8 'MANUEL' POI1;
*cntcc = retir 'DIFF' cntcc;
 
sor2c = 'DOMA' $sor2 maillage ;
sor2c = 'CHANGER' sor2c POI1;
 
retir = p3 'MANUEL' POI1;
sor2c = retir 'DIFF' sor2c;
 
*
 
 
*
* *********************************
* PROCEDURE CALCULANT LE MODLO 
* *********************************
*
DEBPROC MODULO i*ENTIER j*ENTIER ;
k*ENTIER = i / j ;
mod = i - ( k * j ) ;
FINPROC mod ;
 
*
* ***************************************
* PROCEDURE DE CALCUL DES COEFFICIENTS
* DES EQUATIONS DE QUANTITE DE MOUVEMENT       
* ***************************************
*
DEBPROC ACTU  RV*'TABLE' $mt*'MMODEL' ;
 
* constantes du probleme
roc     = RV.INCO.'ROC' ;
rod     = RV.INCO.'ROD' ;
muc     = RV.INCO.'MUC' ;
mud     = RV.INCO.'MUD' ;
expo    = RV.INCO.'EXPOS' ;
farcc   = RV.INCO.'FARCC' ;
farcd   = RV.INCO.'FARCD' ;
 
* LECTURE DES INCONNUES *
 
JC = RV.INCO.'JC' ;
JD = RV.INCO.'JD' ;
 
* implicitation en alpha pour les frottements.
ALPM  = RV.INCO.'AM' ;
*ALPD = 'COPIER' RV.INCO.'AD';
ALPD = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
RV.INCO.'AC' = kops 1.D0 '-' RV.INCO.'AD';
 
MASK1D = ALPD MASQUE 'INFERIEUR' (1.D0 '-' 1.D-15)  ;
MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D  ;
ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
          '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
MASK1D = ALPD MASQUE 'SUPERIEUR' 1D-15  ;
MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
      '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
 
ALPC =  'KOPS' 1.D0 '-' ALPD;
 
*** COUPLAGE DYNAMIQUE ***
 
* VISCOSITE DYNAMIQUE EFFICACE
MUEFF = MUC * RV.INCO.'MUCOEF'  ;
RV.INCO.'MUEFF' = 'KCHT' $mt SCAL 'SOMMET' MUEFF ;
RV.INCO.'SNUFFC' = 'KOPS' MUEFF '/' ROC ;
RV.INCO.'SNUFFC' = 'KOPS' ALPC * RV.INCO.'SNUFFC';
RV.INCO.'NUEFFC' = 'KCHT' $mt scal centre
             ('NOEL' $mt RV.INCO.'SNUFFC');
RV.INCO.'SNUFFD' = 'KOPS' MUEFF '/' ROD ; 
RV.INCO.'SNUFFD' = 'KOPS' ALPD * RV.INCO.'SNUFFD';
RV.INCO.'NUEFFD' = 'KCHT' $mt scal centre
             ('NOEL' $mt RV.INCO.'SNUFFD');
 
 
** UC ET UD : VITESSES INSTANTANEES **
 
SEUIL = RV.INCO.'AD' 'MASQUE' 'SUPERIEUR' 0.05D0;
SEUIL = 'KCHT' $mt SCAL SOMMET SEUIL            ;
 
ALPCX     = NOMC 'UX'    ALPC                  ;
ALPCY     = NOMC 'UY'    ALPC                  ;
ALPCV     = ALPCX 'ET' ALPCY                    ;
UC = 'KOPS' RV.INCO.'JC' '/' ALPCV                       ;
UC = 'KOPS' SEUIL '*' UC                       ;
RV.INCO.'UC'  = 'KCHT' $mt VECT 'SOMMET' UC     ;
 
ALPDX     = NOMC 'UX'     ALPD                 ;
ALPDY     = NOMC 'UY'     ALPD                 ;
ALPDV     = ALPDX 'ET' ALPDY                    ;
UD = 'KOPS' RV.INCO.'JD' '/' ALPDV                       ;
UD = 'KOPS' SEUIL '*' UD                       ;
RV.INCO.'UD'  = 'KCHT' $mt VECT 'SOMMET' UD     ;
 
 
* COEFF DE FROTTEMENT
*on calcule une fraction volumique ALPD prenant
*en compte les parois verticales = max (alpd alpdm).
maskp  = ALPD 'MASQUE' 'SUPERIEUR' ALPM ;
maskp  = 'KCHT' $mt scal sommet maskp;
maskpi = 'KOPS' 1.D0 '-' maskp ;
ALPDF  = 'KOPS' maskp * ALPD ;
ALPFF  = 'KOPS' maskpi * ALPM ;
ALPDF  = 'KOPS' ALPDF '+' ALPFF ;
 
*HINDER =  ALPDF ;
HINDER  =  0.5*((KOPS 0.001D0 '+' ALPC) 
           + ( abs (KOPS 0.001D0 '-' ALPC)));
*HINDER = ALPC;
HINDER = 'KCHT' $mt 'SCAL' 'SOMMET' HINDER ;
HINDER = 'KOPS' HINDER '**' expo ;
HINDER = 'KOPS' MUC '*' HINDER ;
kft = 'KCHT' $mt SCAL 'SOMMET' (9.D0/2.D0) ;
kft = 'KOPS' ('KOPS' RV.INCO.'AI' '**' 2) '*' kft;
 
RV.INCO.'FRCD' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' -1.D0 
    '*'  ('KOPS' kft '/' roc) ) ;    
RV.INCO.'FRCD1' = 'KOPS' RV.INCO.'FRCD' '*' HINDER ;
RV.INCO.'FRCD2' = -1.0D0 * ('KOPS' RV.INCO.'FRCD1' '*' ALPD) ;
RV.INCO.'FRCD' = 'KOPS' RV.INCO.'FRCD1' '+' RV.INCO.'FRCD2' ;
RV.INCO.'MFRCD' = -0.0D0 * ('KOPS' RV.INCO.'FRCD1' '*' RV.INCO.'JD') ;
 
RV.INCO.'FRDC' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' -1.D0
    '*' ('KOPS' kft '/' rod)) ;
RV.INCO.'FRDC' = 'KOPS' RV.INCO.'FRDC' '*' HINDER ; 
RV.INCO.'MFRDC' = 0.0D0 * ('KOPS' RV.INCO.'FRDC' '*' RV.INCO.'JC') ;
RV.INCO.'FRDC' = 1.0D0 * ('KOPS' RV.INCO.'FRDC' '*' ALPD) ;
 
RV.INCO.'FRCC' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' 1.D0 
    '*'  ('KOPS' kft '/' roc) ) ;    
RV.INCO.'FRCC' = 'KOPS' RV.INCO.'FRCC' '*' HINDER ;
RV.INCO.'MFRCC' = 0.0D0 * ('KOPS' RV.INCO.'FRCC' '*' RV.INCO.'JC') ;
RV.INCO.'FRCC' = 1.0D0 * ('KOPS' RV.INCO.'FRCC' '*' ALPD) ;
 
RV.INCO.'FRDD' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' 1.D0
    '*' ('KOPS' kft '/' rod)) ;
RV.INCO.'FRDD1' = 'KOPS' RV.INCO.'FRDD' '*' HINDER ; 
RV.INCO.'FRDD2' = -1.0D0 * ('KOPS' RV.INCO.'FRDD1' '*' ALPD) ;
RV.INCO.'FRDD' = 'KOPS' RV.INCO.'FRDD1' '+' RV.INCO.'FRDD2' ;
RV.INCO.'MFRDD' = -0.0D0 * ('KOPS' RV.INCO.'FRDD1' '*' RV.INCO.'JD') ;
 
RV.INCO.'KFT' = 'KCHT' $mt SCAL SOMMET kft;
 
 
* termes de pression 
 
PUC = 'KOPS' -1.D0 '*' RV.INCO.'AC'  ;
PUC = 'KOPS' PUC '/' ROC ;
RV.INCO.'PUC' = kcht $mt SCAL SOMMET PUC ;
 
PUD = 'KOPS' -1.D0 '*' RV.INCO.'AD'  ;
PUD = 'KOPS' PUD '/' ROD ;
RV.INCO.'PUD' = kcht $mt SCAL SOMMET PUD ;
RV.INCO.'PUD' = 1.0 * RV.INCO.'PUD';
 
PNN = ELNO $mt RV.INCO.PRES CENTREP1;
PNN = KCHT $mt SCAL SOMMET PNN;
GPNN =  KOPS GRAD PNN $MT;
GPNN = KCHT $mt VECT CENTRE GPNN;
GPNN = elno $mt GPNN;
RV.INCO.'GPUD' = kcht $mt VECT SOMMET (0.5D0 * (1/ROC) * GPNN);
 
* terme source pesanteur signe + si TOIMP
* signe - si MDIA ou dans DFDT. on met de l'ordre
* 2, au lieu de alpc alpd(n+1), on a
* (alpc-alpd)alpd(n+1) + alpd*alpd
 
farchc = 'KOPS' ('KOPS' RV.INCO.'AC' '-' RV.INCO.'AD') '*' farcc;
farchci = -1.D0 * farchc ;
*partie expl
farchce = 'KOPS' ('KOPS' RV.INCO.'AD' * RV.INCO.'AD') '*' farcc;
farchce = -1.D0 * farchce;
RV.INCO.'FARCHCI' =  'KCHT' $mt VECT 'SOMMET' farchci ;
RV.INCO.'FARCHCE' =  'KCHT' $mt VECT 'SOMMET' farchce ;
 
farchd = 'KOPS' ('KOPS' RV.INCO.'AC' '-' RV.INCO.'AD') '*' farcd;
farchdi = -1.D0 * farchd ;
*partie expl
farchde = 'KOPS' ('KOPS' RV.INCO.'AD' * RV.INCO.'AD') '*' farcd;
farchde = -1.D0 * farchde;
RV.INCO.'FARCHDI' =  'KCHT' $mt VECT 'SOMMET' farchdi ;
RV.INCO.'FARCHDE' =  'KCHT' $mt VECT 'SOMMET' farchde ;
 
* terme correction laplacien sur ortho radiale
 
RV.INCO.'TCORIC' = 'KOPS' RV.INCO.'SNUFFC' '*' RV.INCO.'INVCAR';
RV.INCO.'TCORID' = 'KOPS' RV.INCO.'SNUFFD' '*' RV.INCO.'INVCAR';
 
FINPROC RV;
 
*************************
*** DONNEES PHYSIQUES ***
*************************
* vitesse angulaire en rad/s
omeg = 100.0D0 ;
* rotation, coef de la force de Coriolis à gauche de l'égalité pour r,z
domeg = ((-2.D0*omeg) 0.D0) ;
* rotation, coef de la force de Coriolis à gauche de l'égalité pour theta
domegt = ((2.D0*omeg) 0.D0);
* champs des 1/r2
invrcar = 'KCHT' $mt SCAL SOMMET ('COORDONNEE' 1 &mt);
invrcar = 'KOPS' invrcar * invrcar;
uni = 'KCHT' $mt SCAL SOMMET 1.D0;
invrcar = 'KOPS' uni '/' invrcar ;
* accélération centrifuge
accelr = 'KCHT' $mt SCAL SOMMET ('COORDONNEE' 1 &mt);
omeg2  = omeg * omeg;
omeg2 = (omeg2 0.D0);
accelr = 'KOPS' omeg2 '*' accelr;
* formule Ishii trainée
expo = -4.1D0;
* formule viscosité
expomu = -0.0D0;
*coeff de viscosité (pour lisser avec mauvais maillage)
MUCOEF = 10.D0 ;
*Densité volumique
ROC = 1000.D0   ;
ROD = 1100.D0   ;
epsilc = (ROC - ROD) / ROC ;
epsild = (ROD - ROC) / ROD ;
* Viscosite cinématique et dynamique de la phase liquide
MUC  = 0.001D0 ;
NUC  = MUC / ROC ;
MUD  = 0.001D0 ;
NUD  = MUD / ROD ;
* Rayon initial des gouttes
ray0 = 1.D-5 ;
* Concentration initiale de la phase dispersée
ALPD0 = 0.3D0 ;
* Concentration initiale de la phase continue
ALPC0 = 1 - ALPD0 ;
*C0 = 1.D-2 ;
C0 = 1.D-2;
* Coeff. de Pénalisation *
EPSS=1.D-9;
* Valeur minimale de EPSN *
MIEPS = 1.e-5 ;
* Valeur maximale de EPSN *
MAEPS = 0.99999 ;
* Vicosité maximale et truc maximal *
ALPMAX = 0.9D0 ;
MUMAX = 1.D-1 ;
RAY2  = 1.D-3;
* fraction vol max
alphadm = 0.9900000001D0 ;
alphacm = MIEPS ;
* profil de fraction volumique
corx = 'COORDONNEE' 1 &mt;
alpd1 = 'TANH' (20000.D0 * (corx '-' (7.95*rext/10.D0))) ;
alpd1 = (0.998 '-' ALPD0) *  (1.D0 '+' alpd1) '/' 2.D0;
alpd2 = 'TANH' (20000.D0 * (corx '-' (7.05D0 * rext/10.D0))) ;
alpd2 =  ALPD0 *  (1.D0 '+' alpd2) '/' 2.D0;
alpd = 0.001D0 '+' alpd1 '+' alpd2;
alpd = 'KCHT' $mt scal sommet alpd;
alpini = alpd;
* concentration aire interf
ALPI0 = 3.D0 * alpd / ray0 ;
*'TRACER' ALPD &mt;
* poussée d'Archimède multipliée par +1 car TOIMP à droite de l'égaZZlité
gravos = 'KCHT' $mt VECT SOMMET (0.D0 -0.D0) ;
archi  = 'KOPS' gravos '+' accelr ;
farcc = epsilc* archi ;
farcd = epsild* archi ;
 
* PAS DE TEMPS *
 
DT = DELTATPS ;
 
************************************
*** SYSTEMES D'EQUATIONS ET C.L. ***
************************************
*
*
*************************************
***    SYSTEME PORTANT SUR :      ***
***          UN, VN, P            ***
*************************************
 
* EQUATION DE CONSERVATION DE LA QUANTIE DE MVT PHASE CONT
*
* composantes r,z
RV = 'EQEX' $mt ITMA 1  
OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'CENTREE'
ZONE $mt OPER 'DFDT' 1.D0 'JCP' 'JCPP' 'DT' 'UC' 1.D-15 INCO 'JC'
ZONE $mt OPER 'MDIA'  'FARCHCI'        INCO 'AD' 'JC'
;
RV = 'EQEX' RV 
OPTI 'EFM1' 'IMPL' 'CENTREE'
ZONE $mt OPER 'MDIA' 'FRCD' INCO 'JD' 'JC'
ZONE $mt OPER 'MDIA' 'FRCC' INCO 'JC' 'JC'
ZONE $mt OPER 'MDIA' domeg  INCO 'CT' 'JC'
OPTI 'EF' 'IMPL' 'CENTREE'
ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JC'
ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JD'
*ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JC'
;
 
* composante theta, otho radiale
RV = 'EQEX' RV  
OPTI 'EFM1' 'IMPL' 'BDF2' 'CENTREE'
ZONE $mt OPER 'DFDT' 1.D0 'CTP' 'CTPP' 'DT' 'CT' 1.D-15  INCO 'CT'
ZONE $mt OPER 'MDIA' 'FRCD' INCO 'JDT' 'CT'
ZONE $mt OPER 'MDIA' 'FRCC' INCO 'CT' 'CT'
'OPTI' 'EFM1' 'IMPL' 'CENTREE'
ZONE $mt OPER 'MDIA' 'domegt' INCO 'JC' 'CT'
ZONE $mt OPER 'MDIA' 'TCORIC' INCO 'CT' 'CT'
ZONE $mt OPER 'MDIA' 'TCORIC' INCO 'JDT' 'CT'
;
 
RV = EQEX  RV
OPTI 'EF' 'IMPL' 'CENTREE'
ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'CT'
ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JDT'
*ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'CT'
OPTI 'EF'  'IMPL' 'CENTREP1' 'CENTREE'
ZONE $mt OPER 'KMBT'           'PUC'        INCO 'PRES' 'JC'
;
 
*
* EQUATION DE CONSERVATION DE LA QUANTITE DE MVT PHASE DISP
*
* composantes r, z
RV = EQEX RV  
OPTI 'EFM1' 'IMPL' 'BDF2' 'CONS' 'CENTREE'
ZONE $mt OPER 'DFDT' 1.D0 'JDP' 'JDPP' 'DT' 'UD' 1.D-15 INCO 'JD'
ZONE $mt OPER 'MDIA'  'FARCHDI'        INCO 'AD' 'JD'
ZONE $mt OPER 'MDIA' 'FRDC' INCO 'JC' 'JD'
ZONE $mt OPER 'MDIA' 'FRDD' INCO 'JD' 'JD'
ZONE $mt OPER 'MDIA' domeg  INCO 'JDT' 'JD'
;
 
* ATTENTION LES LAPLACIENS PAS BONS, LAPN JC + JD -> JC
* ET RIEN POUR JD (A VERIF)
* composante ortho-radiale theta
RV = EQEX RV  
OPTI 'EFM1' 'IMPL' 'BDF2' 'CENTREE'
ZONE $mt OPER 'DFDT' 1.D0 'DTP' 'DTPP' 'DT' 'CT' 1.D-15  INCO 'JDT'
ZONE $mt OPER 'MDIA' 'FRDC' INCO 'CT' 'JDT'
ZONE $mt OPER 'MDIA' 'FRDD' INCO 'JDT' 'JDT'
ZONE $mt OPER 'MDIA' 'domegt' INCO 'JD' 'JDT'
ZONE $mt OPER 'MDIA' 'TCORID' INCO 'CT' 'JDT'
ZONE $mt OPER 'MDIA' 'TCORID' INCO 'JDT' 'JDT'
;
 
 
RV = EQEX RV
OPTI 'EF' 'IMPL' 'CENTREE'
ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'CT'
ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JDT'
*ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JDT'
OPTI 'EF' 'IMPL'  'CENTREE'
ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JD'
ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JC'
*ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JD'
;
 
 
RV = EQEX  RV
OPTI 'EF'  'IMPL' 'CENTREP1' 'CENTREE'
ZONE $mt OPER 'KMBT'          'PUD'         INCO 'PRES' 'JD'
;
 
*
* EQUATION DE CONSERVATION DU VOLUME (INCOMPRESSIBLE)
*
 
RV = EQEX  RV
OPTI 'EFM1'  'IMPL' 'CENTREP1' 'CENTREE' 
ZONE $mt OPER 'KMAB'            1.D0  0.        INCO 'JD' 'PRES'
ZONE $mt OPER 'KMAB'            1.D0  0.        INCO 'JC' 'PRES'
;
 
 
* 
* EQUATION DE CONSERVATION DE LA FRACTION VOLUMIQUE
*
 
RV = EQEX RV  
OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'SUPG' 'CMD' 0.2 
ZONE $mt OPER 'DFDT' 1.D0 'ADP' 'ADPP' 'DT' 'UD' 1.D-15 INCO 'AD'
OPTI 'EF'  'IMPL' 'CENTREP1' 'CONS' 'SUPG' 'CMD' 0.2
ZONE $mt OPER 'KMAB' 1.D0 0.D0 'SOMM' 1.D-15 INCO 'JD' 'AD'
;
 
 
*** FIN SECOND SYSTEME ***
**************************
 
*
* CONDITION AUX LIMITES
*
 
RV = EQEX  RV
  'CLIM'
         'JC' 'VIMP' (&bas1 'ET' &bas2 'ET' &hau) 0.  
         'JC' 'UIMP' (&bas1 ET &bas2 ET &hau ET &cot2 ET &cot4) 0.D0
         'JD' 'VIMP' (&bas1 'ET' &bas2 'ET' &hau) 0.  
         'JD' 'UIMP' (&bas1 ET &bas2 ET &hau ET &cot2 ET &cot4) 0.D0
         'JD' 'VIMP' (&entr) 1.65D-4
         'JC' 'VIMP' (&entr) 1.65D-4
                           ;
RV = EQEX  RV
  'CLIM'
          'CT' 'TIMP' &cnt 0.D0
          'JDT' 'TIMP' &cnt 0.D0
;
 
RV = EQEX  RV
   'CLIM'
         'AD' 'TIMP' (&cot2) 0.99  
         'AD' 'TIMP' (&cot4) 0.
         'AD' 'TIMP' &sor1   0.D0
         'AD' 'TIMP' sor2c   0.99D0  
         'AD' 'TIMP' &entr  ALPD0
   ;
 
 
*
*******************************************
* INITIALISATION DES TABLES DES INCONNUES
*******************************************
*
 
*
RV.INCO        = TABLE INCO                       ;
 
* constantes du probleme
RV.INCO.'ROC'   = ROC   ;
RV.INCO.'ROD'   = ROD   ;
RV.INCO.'MUC'   = MUC   ;
RV.INCO.'MUD'   = MUD   ;
RV.INCO.'EXPOS' = expo  ;
RV.INCO.'FARCC' = farcc ;
RV.INCO.'FARCD' = farcd ;
RV.INCO.'ALPHADMAX' = alphadm ;
 
*
RV.'DT'       =  DELTATPS ;
RV.INCO.'DT'  = RV.'DT' ;
RV.INCO.'GDT' = -1.D0 * epsild * 9.81D0 * RV.'DT' ;
 
RV.INCO.'JC'   = 'KCHT' $mt VECT 'SOMMET'  (1.e-25 1.e-25) ;
RV.INCO.'JD'   = 'KCHT' $mt VECT 'SOMMET'  (1.e-25 1.e-25) ;
RV.INCO.'JCP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC'    ;
RV.INCO.'JCPP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC'    ;
RV.INCO.'JDP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD'    ;
RV.INCO.'JDPP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD'    ;
RV.INCO.'CT'   = 'KCHT' $mt SCAL 'SOMMET'  1.e-25  ;
RV.INCO.'CTP'   = 'KCHT' $mt SCAL 'SOMMET'  1.e-25  ;
RV.INCO.'CTPP'   = 'KCHT' $mt SCAL 'SOMMET'  1.e-25  ;
RV.INCO.'JDT'   = 'KCHT' $mt SCAL 'SOMMET'  1.e-25  ;
RV.INCO.'DTP'   = 'KCHT' $mt SCAL 'SOMMET'  1.e-25  ;
RV.INCO.'DTPP'   = 'KCHT' $mt SCAL 'SOMMET'  1.e-25  ;
 
RV.INCO.'PRES' = 'KCHT' $mt SCAL 'CENTREP1'   1.e-5         ;
 
 
RV.INCO.'AC' = 'KCHT' $mt SCAL 'SOMMET'   (1.D0 '-' ALPD)  ;
RV.INCO.'AD' = 'KCHT' $mt SCAL 'SOMMET'   ALPD         ;
RV.INCO.'RALPD' = 'KCHT' $mt SCAL 'SOMMET'   ALPD        ;
RV.INCO.'RALPC' = 'KCHT' $mt SCAL 'SOMMET'   (1.D0 '-' ALPD) ;
RV.INCO.'ADP'=  'KCHT' $mt SCAL 'SOMMET' RV.INCO.'AD'    ;
RV.INCO.'ADPP'= 'KCHT' $mt SCAL 'SOMMET' RV.INCO.'AD'    ;
 
RV.INCO.'RAY'   = 'KCHT' $mt SCAL 'SOMMET'  ray0            ; 
 
* AI est l'inverse du rayon ;
RV.INCO.'AI'  = 'KCHT' $mt SCAL 'SOMMET' (1.D0 '/' ray0)    ;
 
RV.INCO.'FARCHCI' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
RV.INCO.'FARCHDI' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
RV.INCO.'FARCHCE' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
RV.INCO.'FARCHDE' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
 
* INITIALISATION DES COEFFICIENTS
 
RV.INCO.'DR2' = kcht $mt VECT 'CENTRE' ((d1*d1) (d1*d1));
RV.INCO.'MUCOEF' = kcht $mt SCAL 'SOMMET' MUCOEF;
*
RV.INCO.'CNUL' = 'KCHT' $mt SCAL 'SOMMET' 0.D0            ;
RV.INCO.'JNUL' = 'KCHT' $mt VECT 'SOMMET' (1.D-25 1.D-25) ;
RV.INCO.'JUNI' = 'KCHT' $mt VECT 'SOMMET' (1.D0 1.D0)     ; 
 
RV.INCO.'CNUL1' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
RV.INCO.'CNUL2' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
RV.INCO.'CNUL3' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
RV.INCO.'CNUL4' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
RV.INCO.'SOUC' = 'KCHT' $mt VECT 'CENTRE' (1.e-25 1.e-25) ;
RV.INCO.'SOUD' = 'KCHT' $mt VECT 'CENTRE' (1.e-25 1.e-25) ;
 
RV.INCO.'PUC'  = 'KCHT' $mt SCAL 'SOMMET'  1.e-5          ;
RV.INCO.'PUD'  = 'KCHT' $mt SCAL 'SOMMET'  1.e-5          ;
RV.INCO.'GPUD'  = 'KCHT' $mt VECT 'SOMMET'  (0.D0 0.D0)   ;
 
RV.INCO.'UC'   = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
RV.INCO.'UD'   = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
 
RV.INCO.'NUEFFC' = 'KCHT' $mt SCAL 'CENTRE' 0.0D0         ;
RV.INCO.'NUEFFD' = 'KCHT' $mt SCAL 'CENTRE' 0.0D0         ;
RV.INCO.'SNUFFC' = 'KCHT' $mt SCAL 'SOMMET' 0.0D0         ;
RV.INCO.'SNUFFD' = 'KCHT' $mt SCAL 'SOMMET' 0.0D0         ;
 
RV.INCO.'MUDIFC' = 'KCHT' $mt SCAL 'CENTRE' 0.D0          ;
RV.INCO.'MUDIFD' = 'KCHT' $mt SCAL 'CENTRE' 0.D0          ;
 
RV.INCO.'FRCC'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRCD1'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRCD2'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRCD'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRDC'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRDD1'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRDD2'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'FRDD'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
RV.INCO.'MFRCC'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
RV.INCO.'MFRCD'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
RV.INCO.'MFRDC'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
RV.INCO.'MFRDD'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
 
RV.INCO.'TCORIC' = 'KCHT' $mt VECT SOMMET (0.D0 0.D0)      ;
RV.INCO.'TCORID' = 'KCHT' $mt VECT SOMMET (0.D0 0.D0)      ;
 
RV.INCO.'INVCAR' = 'KCHT' $mt SCAL SOMMET invrcar           ;
RV.INCO.'domegt' = domegt;
 
*COEFFICIENT DE RELAXATION
ALFI = 1.0D0 ;
ALFID = 1.0D0 ;
ALFII = 1.D0 ;
 
* Compteur des pas de temps
ntemp = 0 ;
 
*** Calcul de ALPM frac vol pour la paroi ***********************
*** distance à la paroi la plus proche **************************
 
corx = 'COOR' &mt 1;
corx = 'KCHT' $mt scal sommet corx;
corxi = 'KOPS' rext '-' corx;
corx = 'KOPS'  corx '-' rint; 
masko  = corx 'MASQUE' 'INFERIEUR' corxi ;
masko = 'KCHT' $mt scal sommet masko;
maskoi = 'KOPS' 1.D0 '-' masko;
corx   = 'KOPS' masko * corx ;
corxi  = 'KOPS' maskoi * corxi ;
corx  = 'KOPS' corx '+' corxi ;
 
ratio = 'KOPS' corx '/' RV.INCO.'RAY' ;
ratio = 'KOPS' ratio '+' 3.D0 ;
drel  = 2.D0 * (alphadm**0.34);
ratio = 'KOPS' ratio '/' drel;
* on fabrique une fraction volumique telle que les écarts
* entre particules soient identiques à la distance à la paroi
* la plus proche.
RV.INCO.'AM' = 'KOPS' ratio ** (-3);
 
*'TRACER' RV.INCO.'AM' &mt ('CONTOUR' &mt);
 
 
*
**** BOUCLES DES ITERATIONS INTERNES ET EN TEMPS ****************
*
REPETER BLOCKT NITMA ;
*
ntemp = ntemp + 1 ;
'MESSAGE' 'Itération externe N°' ntemp;
 
ntempi = 0 ;
 
* On stocke les variables au pas de temps ntemp
VJCP =  'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JCP' ;
VJCPP =  'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JCPP' ;
VJDP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JDP' ;
VJDPP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JDPP' ;
VJCTP = 'COPIER' RV.INCO.'CTP' ;
VJCTPP = 'COPIER' RV.INCO.'CTPP' ;
VJDTP = 'COPIER' RV.INCO.'DTP' ;
VJDTPP = 'COPIER' RV.INCO.'DTPP' ;
VALPDP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADP' ;
VALPDPP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADPP';
 
* On effectue les itérations de point fixe
*
'REPETER' BLOCKI NITER ;
*
ntempi = ntempi + 1 ;
*
*
*************************************
***    SYSTEME PORTANT SUR :      ***
***          UN, VN, P            ***
*************************************
*
* On charge les résultats du vrai pas de temps précédent
RV.INCO.'JCP' = 'KCHT' $mt 'VECT' SOMMET VJCP;
RV.INCO.'JCPP' = 'KCHT' $mt 'VECT' SOMMET VJCPP;
RV.INCO.'CTP' = 'KCHT' $mt SCAL SOMMET VJCTP;
RV.INCO.'CTPP' = 'KCHT' $mt SCAL SOMMET VJCTPP;
RV.INCO.'JDP' = 'KCHT' $mt 'VECT' SOMMET VJDP;
RV.INCO.'JDPP' = 'KCHT' $mt 'VECT' SOMMET VJDPP;
RV.INCO.'DTP' = 'KCHT' $mt SCAL SOMMET VJDTP;
RV.INCO.'DTPP' = 'KCHT' $mt SCAL SOMMET VJDTPP;
RV.INCO.'ADP' = 'KCHT' $mt SCAL SOMMET VALPDP ;
RV.INCO.'ADPP' = 'KCHT' $mt SCAL SOMMET VALPDPP ;
* resultats de l'iter interne préc.
JCIP = 'KCHT' $mt VECT SOMMET RV.INCO.'JC' ;
JDIP = 'KCHT' $mt 'VECT' SOMMET RV.INCO.'JD' ;
ALPDIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AD' ;
ALPCIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AC' ;
 
*RAYC = 'KOPS' ALPCIP '/' ('KOPS' 1.D-15 '+' ALPIIP) ;
*RAYC = 'KOPS' ('KOPS' ray0 '+' RAYC) '+' ('ABS' ('KOPS' ray0 '-' RAYC)); 
*RAYC = 'KCHT' $mt SCAL 'SOMMET' RAYC ;
*RAYC = 'KOPS' RAYC '*' 3 ;
*RV.INCO.'RAYC' = RAYC ;
 
 
*
* On réactualise les données.
*
 
RV = ACTU RV $mt;
 
* On affecte le terme source expl de QDM par le biais de DFDT
* coef 2 car ordre 2 en temps dans dfdt. JCPP sinon la poussee
* archi est rechargéee * fois lors des reinitialisations
RV.INCO.'JCPP'  = 'KOPS' RV.INCO.'JCP' '+'
                 ('KOPS' (2.D0*RV.INCO.'DT') * RV.INCO.'FARCHCE') ;
 
RV.INCO.'JDPP'  = 'KOPS' RV.INCO.'JDP' '+'
                 ('KOPS' (2.D0*RV.INCO.'DT') * RV.INCO.'FARCHDE') ;
 
 rv . 'METHINV' . 'TYPINV'  = 2 ;       
EXEC RV ;
'MESSAGE' 'Boucle interne N°' ntempi; 
 
ALPD = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
MASK1D = ALPD MASQUE 'INFERIEUR' (1.D0 '-' 1.D-15)  ;
MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D  ;
ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
MASK1D = ALPD MASQUE 'SUPERIEUR' 1D-15  ;
MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
 
* On masque AD
RV.INCO.'AD' = ALPD;
 
 
 
* CALCUL DE L'ERREUR RELATIVE MAX. *
ERR = 0.D0 ;
ERRC = 'KOPS' ('KOPS' RV.INCO.'JC' '-' JCIP)
          '/' ('MAXIMUM' ('ABS' JCIP) '+' 1.D-10) ;
ERRD = 'KOPS' ('KOPS' RV.INCO.'JD' '-' JDIP)
          '/' ('MAXIMUM' ('ABS' JDIP) '+' 1.D-10) ;
ERRC = 'MAXIMUM' ('ABS' ERRC) ;
ERRD = 'MAXIMUM' ('ABS' ERRD) ;
'SI' (ERRC > ERR) ;
ERR = ERRC ;
VAR = 'MOTS' JC ;
'FINSI' ;
 
'SI' (ERRD > ERR) ;
ERR = ERRD ;
VAR = 'MOTS' JD ;
'FINSI' ;
ERRA = 'KOPS' ('KOPS' RV.INCO.'AD' '-' ALPDIP)
          '/' ('MAXIMUM' ('ABS' ALPDIP) '+' 1.D-10) ;
ERRA = 'MAXIMUM' ('ABS' ERRA) ;
'SI' (ERRA > ERR) ;
ERR = ERRA ;
VAR = 'MOTS' ALPHAD ;
'FINSI' ;
 
* CALCUL DE L'ERREUR RELATIVE MAX. *
 
 
*'LIST' VAR;
'MESSAGE' ERR; 
'MESSAGE' ERRC ERRD ERRA;
'SI' (ERR <EG PREC) ;
QUITTER BLOCKI ;
'FINSI'       ;
 
*
'FIN' BLOCKI ;
*
 
'MESSAGE' 'NBR ITER INTERNES' ntempi;
'MESSAGE' 'ERREUR RELATIVE MAX' ERR ;
'MESSAGE' 'Max JC' ('MAXIMUM' (RV.INCO.'JC'))   ;
'MESSAGE' 'Min JC' ('MINIMUM' (RV.INCO.'JC'))   ;
'MESSAGE' 'Max JD' ('MAXIMUM' (RV.INCO.'JD'))   ;
'MESSAGE' 'Min JD' ('MINIMUM' (RV.INCO.'JD'))   ;
'MESSAGE' 'Max AlpD' ('MAXIMUM' (RV.INCO.'AD')) ;
'MESSAGE' 'Min AlpD' ('MINIMUM' (RV.INCO.'AD')) ;
 
 
SI (EGA GRAPH O);
   'SI' ('EGA' (ntemp MODULO (NITMA/PERIO)) 0 ) ;
 
   ung = vect rv.INCO.'JD' 25. ux uy jaune;
   ungu = 'EXCO' rv.inco.'JD' UX ;
   ungv = 'EXCO' rv.inco.'JC' UX ;
   ungd = 'EXCO' rv.inco.'UD' UX ;
 
   Ealpd  = 'EVOL' (BLEU) 'CHPO' rv.inco.'AD' alpdco ;
   ealpini = 'EVOL' (BLEU) 'CHPO' alpini alpdco;
   Ejd    = 'EVOL' (JAUNE) 'CHPO' ungu alpdco ;
   Ejc    = 'EVOL' (JAUNE) 'CHPO' ungv alpdco ;
   Eeff = 'EVOL' (VERT) 'CHPO' rv.inco.'MUEFF' alpdco ;
   Etransp = 'EVOL'(TURQ) 'CHPO' ungd alpdco ;
 
   'DESSIN' Ejd 'TITRE' 'Jd';
   'DESSIN' Ejc 'TITRE' 'Jc';
   'DESSIN' Etransp 'TITRE' 'Champ transportant Ud' ;
 
   toto = RV.INCO.'AD';
   Ealpd = 'EVOL' 'CHPO' toto alpdco ;
   'DESSIN' (Ealpd  'ET' ealpini) 'TITR'
    'Fraction volumique de la phase Dispersee' ;
 
   uref = 2.0D0;
   ung = vect RV.INCO.'JD' uref ux uy jaune ;
   trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JD' 'NCLK' ;
 
   uref = 2.0D0;
   ung = vect RV.INCO.'JC' uref ux uy jaune ;
   trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JC' 'NCLK' ;
 
   pnn = elno $mt RV.INCO.PRES CENTREP1;
   Epp = 'EVOL'(TURQ) 'CHPO' pnn alpdco ;
   'DESSIN' Epp 'TITRE' 'pression';
 
   ALPD   = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
   MASK1D = 'MASQUE' ALPD 'INFERIEUR' (1.D0 '-' 1.D-15)  ;
   MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D  ;
   ALPD   = 'KOPS' ('KOPS' MASK1D '*' ALPD)
            '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
   MASK1D = 'MASQUE' ALPD 'SUPERIEUR' 1D-15  ;
   MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
   ALPD   = 'KOPS' ('KOPS' MASK1D '*' ALPD)
            '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
   MONTAGNE ALPD mt 0.12 ('PROG' -2. -2. 1.) 'TITRE' 'tot' 'SUPER' ;
 
  'FINSI' ;
 
'FINSI';
 
*
'FIN' BLOCKT ;
 
fin;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales