* fichier : difasyk2D.dgibi
opti dime 2 elem qua8 mode plan trac ps ;
**************************************************************************
*
* 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.
*
*
*
* 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.
*
* 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 présent cas-test, pour des raisons
* d'obtention de résultats analytiques, laisse la
* viscosité efficace constante quelque soit la valeur de
* AD. Le front de sédimentation est alors très raide
* et notre problème revient regarder évoluer le front
* de sédimentation de billes tombant toujours à la
* même vitesse puis s'arrêtant brutalement lorsqu'elles
* rencontrent la zone sédimentée.
*
*
* Le maillage fait une maille d'épaisseur pour accélérer
* la simulation dans le cadre des tests du répertoire
* castem/dgibi.
*
*
* Auteur : Gilles Bernard-Michel
*
* Révision : version 1
*
* date : 21/12/99
*
************************ Options générales *********************************
*
*
** Sortie graphique.
*
GRAPH = 'O';
*
*test long ou court.
*
COMPLET = FAUX;
*
*période des sorties graphiques : tous les NITMA/PERIO pas de temps
*
PERIO = 1 ;
*
* nombre de mailles, pas de temps , intervalle de temps.
*
SI (NON COMPLET) ;
nptx = 1. ;
npty = 10. ;
DELTATPS = 2.;
* nombre pas de temps
NITMA = 10 ;
* nombre points fixes
NITER = 3 ;
* convergence point fixe.
PREC = 1.D-4 ;
SINON;
nptx = 1;
npty = 30;
DELTATPS = 0.5D0;
* nombre pas de temps
NITMA = 40 ;
* nombre points fixes
NITER = 5 ;
* convergence point fixe.
PREC = 1.D-4 ;
FINSI;
*
************************ On définit les propriétés physiques ***************
*
*
** viscosité dynamique mu/rho.
*
nu =1.E-6 ;
*
********************* On traite les aspects géométriques ********************
*
*
** Les paramètres de la géométrie.
*
rint = 0. ;
rext = 0.1;
hh = 0.5 ;
*
** Densité du maillage.
*
dx = (rext - rint) / nptx;
dy = hh / npty ;
*
** Les sommets du bol
*
p1 = rint 0. ;
p2 = rext 0. ;
p3 = rext hh ;
p4 = rint hh ;
p6 = 1.0*dx 0.0;
p8 = 1.D0*dx hh;
*
** Assemblage des bords du domaine.
*
bas = p1 'DROIT' dini dx dfin dx p2 ;
cotd = p2 'DROIT' dini dy dfin dy p3 ;
haut = p3 'DROIT' dini dx dfin dx p4 ;
cotg = p4 'DROIT' dini dy dfin dy p1 ;
cnt = bas 'ET' cotd 'ET' haut 'ET' cotg ;
*
** On maille avec des rectangles.
*
mt= bas cotd haut cotg 'DALLER' ;
'TASSER' mt ;
*
** On réoriente les éléments.
*
mt = 'ORIENTER' mt ;
&mt = CHANGER mt QUAF ;
&bas = CHANGER bas QUAF ;
&cotd = CHANGER cotd QUAF ;
&haut = CHANGER haut QUAF ;
&cotg = CHANGER cotg QUAF;
&cnt = CHANGER cnt QUAF ;
'ELIMINATION' 1.e-6 (&mt 'ET' &cnt 'ET' &cotg 'ET' &haut
'ET' &cotg 'ET' &cnt);
$mt = MODE &mt 'NAVIER_STOKES' QUAF ;
*
** On crée le domaine associé au maillage.
*
mt = doma $mt maillage ;
* On crée une droite pour tracer des profils de fraction volumique.
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;
* *********************************
* PROCEDURE CALCULANT LE MODULO
* *********************************
*
DEBPROC MODULO i*ENTIER j*ENTIER ;
k*ENTIER = i / j ;
mod = i - ( k * j ) ;
FINPROC mod ;
*------------------------------------------------------------------------
*-----------------------------------------------------------------------
* - LIREXY3 - -
* PROCEDURE POUR OBTENIR UNE LISTE CONTENANT LES VALEURS DU CHPO AUX -
* POINTS DE DONNEES : premiere étape pour le tracé d'histogramme -
* PROCEDURE EN TOUT POINT COMPARABLE A LIREXY SAUF QUE L'ON A UNE -
* TABLE EN SORTIE : POUR TRACE D'HISTOGRAMMES. -
* (VOIR REMARQUE SUR POINT PROCHE DANS LIREXY). -
*-----------------------------------------------------------------------
* En entrée ... -
* - CHPOD : champ point dont on extrait les valeurs, -
* - MAILD : maillage des points centres où prendre les valeurs -
* - K : le nom (MOT) de la composante du champ que l'on extrait. -
*-----------------------------------------------------------------------
* En sortie ... -
* la liste des valeurs f(xi,yi) -
*-----------------------------------------------------------------------
'DEBPROC' LIREXY3 CHPOD*'CHPOINT' MAILD*'MAILLAGE' K*'MOT';
NPT = NBEL MAILD;
REPETER BACQ NPT;
I = &bacq;
Pt = MAILD POINT I;
Pty = 'COOR' 2 Pt ;
KI = 'EXTRAIRE' CHPOD Pt K;
SI (I EGA 1);
LIPK = 'PROG' KI;
LIPP = 'PROG' Pty;
SINON;
LIPK = LIPK 'ET' ('PROG' KI);
LIPP = LIPP 'ET' ('PROG' Pty);
FINSI;
FIN BACQ;
evo = 'EVOL' 'MANU' LIPP LIPK;
dessin evo;
'FINPROC';
* Idem pour points sommet.
'DEBPROC' LIREXY2 CHPOD*'CHPOINT' MAILD*'MAILLAGE' K*'MOT';
NPT = NBEL MAILD;
NPT = NPT + 1 ;
REPETER BACQ NPT;
I = &bacq;
Pt = MAILD POINT I;
Pty = 'COOR' 2 Pt ;
KI = 'EXTRAIRE' CHPOD Pt K;
SI (I EGA 1);
LIPK = 'PROG' KI;
LIPP = 'PROG' Pty;
SINON;
LIPK = LIPK 'ET' ('PROG' KI);
LIPP = LIPP 'ET' ('PROG' Pty);
FINSI;
FIN BACQ;
evo = 'EVOL' 'MANU' LIPP LIPK;
dessin evo;
'FINPROC';
*
* ***************************************
* 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 *
RJC = RV.INCO.'RJC' ;
RJD = RV.INCO.'RJD' ;
JC = RV.INCO.'JC' ;
JD = RV.INCO.'JD' ;
* implicitation en alpha pour les frottements.
ALPM = RV.INCO.'AM' ;
*ALPD = 'COPIER' RV.INCO.'ALPD';
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 '-' RV.INCO.'AD';
RAY = RV.INCO.'RAY' ;
RAYC = RV.INCO.'RAYC' ;
ALPI = RV.INCO.'RALPI' ;
COAL = RV.INCO.'COAL' ;
*** COUPLAGE DYNAMIQUE ***
* VISCOSITE DYNAMIQUE EFFICACE
MUEFF = ALPD ;
MUEFF = 'KCHT' $mt 'SCAL' 'SOMMET' MUEFF ;
MUEFF = 'KOPS' 1.0D0 '-' MUEFF ;
MUEFF = 'KOPS' MUEFF '**' expo ;
MUEFF = 'KOPS' MUC '*' MUEFF ;
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 relaxées **
SEUIL = (('TANH' (1.D20 *
(ALPC '-' 1.D-3))) '+' 1)/2.D0;
ALPCX = NOMC 'UX' ALPC ;
ALPCY = NOMC 'UY' ALPC ;
ALPCV = ALPCX 'ET' ALPCY ;
UC = 'KOPS' RJC '/' ALPCV ;
UC = 'KOPS' SEUIL '*' UC ;
RV.INCO.'UC' = 'KCHT' $mt VECT 'SOMMET' UC ;
SEUIL = (('TANH' (1.D20 *
(ALPD '-' 1.D-3))) '+' 1)/2.D0;
ALPDX = NOMC 'UX' ALPD ;
ALPDY = NOMC 'UY' ALPD ;
ALPDV = ALPDX 'ET' ALPDY ;
UD = 'KOPS' RJD '/' ALPDV ;
UD = 'KOPS' SEUIL '*' UD ;
RV.INCO.'UD' = 'KCHT' $mt VECT 'SOMMET' UD ;
** VUD : VITESSES INSTANTANEES lissées**
VUD = 'NOEL' $mt UD ;
VUD = 'KCHT' $mt VECT 'CENTRE' VUD ;
VUD = 'ELNO' $mt VUD ;
RV.INCO.'VUD' = 'KCHT' $mt VECT 'SOMMET' VUD ;
* 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.01D0 '+' ALPC)
+ ( abs (KOPS 0.01D0 '-' 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' kft '*' ('KOPS' RV.INCO.'RAY' '**' -2) ;
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' '*'
RV.INCO.'AD') ;
RV.INCO.'FRCD' = 'KOPS' RV.INCO.'FRCD1' '+' RV.INCO.'FRCD2' ;
RV.INCO.'MFRCD' = -0.5D0 * ('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.5D0 * ('KOPS' RV.INCO.'FRDC' '*' RV.INCO.'JC') ;
RV.INCO.'FRDC' = 1.0D0 * ('KOPS' RV.INCO.'FRDC' '*'
RV.INCO.'AD') ;
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.5D0 * ('KOPS' RV.INCO.'FRCC' '*' RV.INCO.'JC') ;
RV.INCO.'FRCC' = 1.0D0 * ('KOPS' RV.INCO.'FRCC' '*'
RV.INCO.'AD') ;
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' '*'
RV.INCO.'AD') ;
RV.INCO.'FRDD' = 'KOPS' RV.INCO.'FRDD1' '+' RV.INCO.'FRDD2' ;
RV.INCO.'MFRDD' = -0.5D0 * ('KOPS' RV.INCO.'FRDD1' '*' RV.INCO.'JD') ;
RV.INCO.'KFT' = 'KCHT' $mt SCAL SOMMET kft;
* COEFFICIENT DE COALESCENCE
DIST = 'KOPS' ALPD '**' (-1.D0/3.D0);
DIST = 'KCHT' $mt SCAL 'SOMMET' DIST ;
DIST = 'KOPS' DIST '*' 2.03D0 ;
DIST = 'KOPS' ((3.D0 ** 0.5D0)/2.D0) '*' DIST ;
DIST = 'KOPS' DIST '-' 2.D0 ;
MDIST = DIST MASQUE 'SUPERIEURE' 1.D-4 ;
MDIST = 'KCHT' $mt SCAL 'SOMMET' MDIST ;
DIST = 'KOPS' DIST '*' MDIST ;
DIST = 'KOPS' DIST '+' ('KOPS'
('KOPS' 1.D0 '-' MDIST) '*' 1.D-4) ;
COAL = 'KOPS' ALPD '/' ALPC ;
COAL = 'KOPS' COAL '*' ECSTA2 ;
COALC = 'KOPS' ALPC '/' ALPD ;
COALC = 'KCHT' $mt SCAL 'SOMMET' COALC ;
COALC = 'KOPS' COALC '*' 1.D-2 ;
COALC = 'KOPS' ('KOPS' 1.D0 - ECSTA2) '*' COALC ;
COAL = 'KOPS' COALC '+' COAL ;
COAL = 'KOPS' COAL '*' C0 ;
RV.INCO.'COAL' = 'KCHT' $mt SCAL 'SOMMET' COAL ;
RV.INCO.'C0' = 'KCHT' $mt SCAL 'SOMMET' C0 ;
RV.INCO.'DIST' = 'KCHT' $mt SCAL 'SOMMET' DIST ;
* 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 ;
* equation du rayon
RV.INCO.'RAYALPI' = (1.D0/3.D0)* RV.INCO.'ALPI';
* terme de stabilisation QDM
RV.INCO.'GRALPC' = kops GRAD RV.INCO.'AD' $mt;
RV.INCO.'GRALPC' = kcht $mt VECT CENTRE RV.INCO.'GRALPC';
RV.INCO.'GRALPC' = KOPS RV.INCO.'DR2' '*' RV.INCO.'GRALPC';
RV.INCO.'GRALPC' = 'ABS' RV.INCO.'GRALPC';
RV.INCO.'GRALPC' = .1D0 * RV.INCO.'GRALPC';
RV.INCO.'GRALPS' = kcht $mt VECT SOMMET (ELNO $mt RV.INCO.'GRALPC');
RV.INCO.'GRALJD' = 'KCHT' $mt scal sommet ('EXCO' uy RV.INCO.'JD');
RV.INCO.'GRALJD' = kops GRAD RV.INCO.'GRALJD' $mt;
RV.INCO.'GRALJD' = kcht $mt VECT CENTRE RV.INCO.'GRALJD';
RV.INCO.'GRALJD' = KOPS RV.INCO.'DR2' '*' RV.INCO.'GRALJD';
RV.INCO.'GRALJDM' = 0.1D0 * RV.INCO.'GRALJD';
RV.INCO.'GRALJD' = -10.D0 * RV.INCO.'GRALJD';
RV.INCO.'GRALJS' = kcht $mt VECT SOMMET (ELNO $mt RV.INCO.'GRALJD');
FINPROC RV;
*************************
*** DONNEES PHYSIQUES ***
*************************
* formule Ishii
*expo= -4.1;
* ON NE FAIT PAS VARIER MUEFF VOIR DESCRIPTION DU MODELE
expo = -0.0D0;
*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-4 ;
* Concentration initiale de la phase dispersée
ALPD0 = 0.3D0 ;
ALPI0 = 3.D0 * ALPD0 / ray0 ;
* Concentration initiale de la phase continue
ALPC0 = 1 - ALPD0 ;
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;
* poussée d'Archimède multipliée par +1 car TOIMP à droite de l'égaZZlité
gravos = 'KCHT' $mt VECT SOMMET (0.D0 -9.81D0) ;
archi = gravos ;
farcc = epsilc* archi ;
farcd = epsild* archi ;
* fraction vol max
alphadm = 0.9900000001D0 ;
alphacm = MIEPS ;
* profil de fraction volumique
cory = 'COORDONNEE' 2 &mt;
alpd1 = 'TANH' (200.D0 * (cory '-' (3.*hh/10.D0))) ;
alpd1 = (0.999 '-' ALPD0) * (1.D0 '-' alpd1) '/' 2.D0;
alpd2 = 'TANH' (200.D0 * (cory '-' (9.D0 * hh/10.D0))) ;
alpd2 = ALPD0 * (1.D0 '-' alpd2) '/' 2.D0;
alpd = alpd1 '+' alpd2;
alpini = alpd;
*'TRACER' ALPD &mt;
* PAS DE TEMPS *
TERDT = 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'
*OPTI 'EF' 'IMPL' 'CONS' 'CENTREE' 'MUCONS'
OPTI 'EF' 'IMPL' 'CONS' 'CENTREE'
*ZONE $mt OPER 'KONV' 1.D0 'UC' 1.D-15 DT INCO 'JC'
ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JC'
;
RV = EQEX RV
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'
OPTI 'EF' 'IMPL' 'CENTREE'
*ZONE $mt OPER 'TOIMP' 'FARCHD' INCO 'JD'
;
RV = EQEX RV
*OPTI 'EF' 'IMPL' 'CENTREE' 'MUCONS'
OPTI 'EF' 'IMPL' 'CENTREE'
ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JC'
ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JD'
OPTI 'EF' 'IMPL' 'CONS' 'SUPG'
*ZONE $mt OPER 'KONV' 1.D0 'UD' 1.D-15 DT INCO 'JD'
;
;
RV = EQEX RV
OPTI 'EF' 'IMPL' 'CENTREP1' 'CENTREE'
*ZONE $mt OPER 'KMBT' 'GRALPS' 0.D0 'SOMM' 1.D-15 INCO 'AD' 'JD'
*ZONE $mt OPER 'KMBT' -0.001 0.D0 'SOMM' 1.D-15 INCO 'AD' 'JD'
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 'JC' 'PRES'
ZONE $mt OPER 'KMAB' 1.D0 0. INCO 'JD' 'PRES'
*OPTI 'EFMC' 'IMPL' 'CENTREE' 'INCOD' CENTREP1
*ZONE $mt OPER 'DFDT' EPSS 'PRES' DT 'JNUL' 1.D+15 INCO 'PRES'
;
*
* EQUATION DE CONSERVATION DE LA FRACTION VOLUMIQUE
*
RV = EQEX RV
OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'CENTREE'
ZONE $mt OPER 'DFDT' 1.D0 'ADP' 'ADPP' 'DT' 'UD' 1.D-15 INCO 'AD'
OPTI 'EF' 'IMPL' 'CENTREP1' 'CONS' 'SUPG' 'CMD' 0.20
ZONE $mt OPER 'KMAB' 1.D0 0.D0 'SOMM' 1.D-15 INCO 'JD' 'AD'
*OPTI 'EF' 'IMPL' 'CENTREE'
*ZONE $mt OPER 'LAPM' 'GRALPC' INCO 'AD'
;
*** FIN SECOND SYSTEME ***
**************************
*
* CONDITION AUX LIMITES
*
RV = EQEX RV
'CLIM' 'JC' 'UIMP' (&cotg 'ET' &cotd) 0.D0
'JC' 'VIMP' (&bas 'ET' &haut) 0.D0
'JD' 'UIMP' (&cotg 'ET' &cotd) 0.D0
'JD' 'VIMP' (&bas 'ET' &haut) 0.D0
'AD' 'TIMP' (&bas) 1.D0
'AD' 'TIMP' (&haut) 0.D0
;
bcp='ELEM' ('DOMA' $mt 'CENTREP1') POI1 (lect 1) ;
RV = EQEX RV
'CLIM'
'PRES' 'TIMP' bcp 0.D0 ;
*
*******************************************
* INITIALISATION DES TABLES DES INCONNUES
*******************************************
*
*
RV.INCO = TABLE INCO ;
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 * DELTATPS ;
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.'RJC' = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC' ;
RV.INCO.'RJD' = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD' ;
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.'GRALPS' = KCHT $mt VECT CENTRE (0.D0 0.D0) ;
RV.INCO.'GRALPC' = KCHT $mt VECT CENTRE (0.D0 0.D0) ;
RV.INCO.'RAY' = 'KCHT' $mt SCAL 'SOMMET' ray0 ;
RV.INCO.'ALPI' = 'KCHT' $mt SCAL 'SOMMET' ALPI0 ;
RV.INCO.'ALPIP' = 'KCHT' $mt SCAL 'SOMMET' RV.INCO.'ALPI' ;
RV.INCO.'RALPI' = 'KCHT' $mt SCAL 'SOMMET' ALPI0 ;
RV.INCO.'RAYALPI' = KCHT $mt SCAL 'SOMMET' (3.D0 * RV.INCO.'ALPI');
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' ((dx*dx) (dy*dy));
*
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.'VUD' = '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.'COAL' = 'KCHT' $mt SCAL 'SOMMET' C0 ;
*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' ;
VALPDP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADP' ;
VALPDPP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADPP';
VALPIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ALPI' ;
* 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.'JDP' = 'KCHT' $mt 'VECT' SOMMET VJDP;
RV.INCO.'JDPP' = 'KCHT' $mt 'VECT' SOMMET VJDPP;
RV.INCO.'ADP' = 'KCHT' $mt SCAL SOMMET VALPDP ;
RV.INCO.'ADPP' = 'KCHT' $mt SCAL SOMMET VALPDPP ;
RV.INCO.'ALPIP' = 'KCHT' $mt SCAL SOMMET VALPIP ;
* 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' ;
ALPIIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ALPI' ;
ALPCIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AC' ;
* Rayon moyen
ECSTA2 = RV.INCO.'RAY' MASQUE 'INFERIEURE' RAY2 ;
ECSTA2 = 'KCHT' $mt SCAL 'SOMMET' ECSTA2 ;
*RAY = 'KOPS' ALPDIP '/' ALPIIP ;
*RAY = 'KCHT' $mt SCAL 'SOMMET' RAY ;
*RAY = 'KOPS' RAY '*' ECSTA2 ;
*RAY = 'KOPS' RAY '*' 3 ;
RAYC = 'KOPS' ('KOPS' 1.D0 '-' ALPCIP) '/' ALPIIP ;
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') ;
EXEC RV ;
'MESSAGE' 'Boucle interne N°' ntempi;
* 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' UY ;
ungv = 'EXCO' rv.inco.'JC' UY ;
ungd = 'EXCO' rv.inco.'UD' UY ;
ealpini = 'EVOL' (BLEU) 'CHPO' alpini alpdco;
Ejd = 'EVOL' (JAUNE) 'CHPO' ungu alpdco ;
Ejc = 'EVOL' (JAUNE) 'CHPO' ungv alpdco ;
toto = RV.INCO.'AD';
Ealpd = 'EVOL' 'CHPO' toto alpdco ;
'DESSIN' (Ealpd 'ET' ealpini) 'TITR'
'Frac vol phase Dispersee - comparaison solu init' ;
uref = 1.D+2;
ung = vect RV.INCO.'JC' uref ux uy jaune ;
trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JC' ;
'DESSIN' Ejc 'TITRE' 'Profil vertical de Jc';
uref = 1.D+2;
ung = vect RV.INCO.'JD' uref ux uy jaune ;
trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JD' ;
'DESSIN' Ejd 'TITRE' 'Profil vertical de Jd ';
pnn = elno $mt RV.INCO.PRES CENTREP1;
Epp = 'EVOL'(TURQ) 'CHPO' pnn alpdco ;
dessin Epp 'TITRE' 'pression réduite';
'FINSI' ;
FINSI;
'FIN' BLOCKT ;
************* calcul de l'erreur ************************************
* calcul de la solution analytique.
* La vitesse volumique max de sédimentation pour JD
jdsolm = ALPD0 '*' (1.D0 '-' ALPD0);
coefsol = 2. * ray0 * ray0 * 9.81 * (rod - roc) / (9. * nu * roc);
jdsolm = coefsol * jdsolm;
* Le profile de fraction volumique au temps final NITMA * DELTATPS
deplabas = (3. * hh/10.) + (jdsolm * nitma * DELTATPS);
deplahau = hh - (hh/10.) - ((deplabas - (3.* hh/10.))/ alpd0);
cory = 'COORDONNEE' 2 &mt;
alpd1 = 'TANH' (20000.D0 * (cory '-' (deplabas))) ;
alpd1 = (0.999 '-' ALPD0) * (1.D0 '-' alpd1) '/' 2.D0;
alpd2 = 'TANH' (200.D0 * (cory '-' (deplahau))) ;
alpd2 = ALPD0 * (1.D0 '-' alpd2) '/' 2.D0;
alpd = alpd1 '+' alpd2;
alpfin = alpd;
* la vitesse de sédimenation pour JD (flux volumique), pour
* alpd numérique donné
jdsol = kops RV.INCO.'AD' '*' (kops 1.D0 '-' RV.INCO.'AD');
coefsol = -2. * ray0 * ray0 * 9.81 * (rod - roc) / (9. * nu * roc);
jdsol = coefsol * jdsol;
* la vitesse volumique verticale calculée numériquement
jdsoln = exco RV.INCO.'JD' uy;
jdsoln = kcht $mt SCAL SOMMET jdsoln;
* idem mais avec la fraction volumique analytique
alpfin = kcht $mt SCAL SOMMET alpfin;
jdsol1 = kops alpfin '*' (kops 1.D0 '-' alpfin);
jdsol1 = coefsol * jdsol1;
* calcul de l'erreur L2 sur la fraction volumique
DEBPROC CALCURR vit*'CHPOINT' vitsol*'CHPOINT' mod*'MMODEL'
mt*'MAILLAGE' coupe*'MAILLAGE';
vit = kcht mod SCAL SOMMET vit;
vitsol = kcht mod SCAL SOMMET vitsol;
udiff = kops vit - vitsol;
udif2 = kops udiff * udiff;
evudif2 = 'EVOL' 'CHPO' udif2 coupe;
ud2 = SOMME evudif2;
ud2 = MAXI ud2;
usol2 = vitsol * vitsol;
eusol2 = 'EVOL' 'CHPO' usol2 coupe;
us2 = SOMME eusol2;
us2 = (MAXI us2) + 1.e-14;
err = ud2 / us2;
err = err**0.5;
FINPROC err ;
* erreur en pression (choque inclus) en norme L2.
errorp = CALCURR RV.INCO.'AD' alpfin $mt mt alpdco;
* erreur en flux vol de vitesse JD (choque inclus)
* par rapport à la solution semi analytique (formule
* analytique faisant intervenir AD calculé numériquement).
errorv = CALCURR jdsoln jdsol $mt mt alpdco;
* erreur pour JD par rapport à la solution analytique
errorv1 = CALCURR jdsoln jdsol1 $mt mt alpdco;
SI (EGA GRAPH O);
trac mt;
Ealpf = 'EVOL' (BLEU) 'CHPO' alpfin alpdco ;
'DESSIN' (Ealpf 'ET' Ealpd) 'TITR'
'Frac vol phase Dispersee - comparaison solu analytique' ;
* on compare maintenant les flux avec le flux analytique mais
* utilisant la fraction volumique calculée numériquement.
* ceci permet d'évaluer la résolution intrinsèque de la quantité
* de mouvement. On pourrait aussi comparer avec la solution
* analytique.
Ejdsol = 'EVOL' (BLEU) 'CHPO' jdsol alpdco;
'DESSIN' (Ejd et Ejdsol) 'TITR'
'Flux vol phase Dispersee - comparaison solu semi-analytique' ;
Ejdsol1 = 'EVOL' (BLEU) 'CHPO' jdsol1 alpdco;
'DESSIN' (Ejd et Ejdsol1) 'TITR'
'Flux vol phase Dispersee - comparaison solu analytique' ;
FINSI;
*
SI (NON COMPLET);
SI (errorp > 0.07D0);
ERREUR 5;
FINSI;
SI (errorv > 1.D-3);
ERREUR 5;
FINSI;
SI (errorv1 > 0.16D0);
ERREUR 5;
FINSI;
SINON;
SI (errorp > 0.03D0);
ERREUR 5;
FINSI;
SI (errorv > 0.001D0);
ERREUR 5;
FINSI;
SI (errorv1 > 0.010);
ERREUR 5;
FINSI;
FINSI;
*OPTI SAUV '/u2/gbm/castem/diphiplan.res' ;
*REPIX RV ;
*REPIX RA ;
*SAUV ;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer