Télécharger la plaquette Cast3M


fr en

Télécharger difasyk2D.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  difasyk2D.dgibi
  2. opti dime 2 elem qua8 mode plan trac ps ;
  3. **************************************************************************
  4. *
  5. *        CAS TEST DIPHASIQUE LIQUIDE-LIQUIDE 2D PLAN
  6. *
  7. ***************************************************************************
  8. *
  9. *
  10. *        Présentation: Ce cas test permet la simulation de la sédimention
  11. *                      de gouttes sphériques de liquide immergé dans
  12. *                      un milieu continu (autre liquide) de densité moindre.
  13. *                      Bien entendu les densités peuvent être modifiées.
  14. *
  15. *
  16. *
  17. *        Le modèle : Le modèle simulé consiste à résoudre 2 équations de la
  18. *                    quantité de mouvement sur les flux volumiques des 2
  19. *                    phases JC et JD. Ensuite on résout une équation de la
  20. *                    conservation de la masse qui porte sur la fraction
  21. *                    volumique de la phase dispersée : AD (pour alpha_D).
  22. *                    La conservation de la fraction voumique de la phase
  23. *                    continue AC s'en déduite par AD + AC = 1.
  24. *                    L'équation d'incompressibilité du fluide est
  25. *                    Div(JC+JD) = 0.
  26. *
  27. *                    le modèle de coalescence est simplifié, ce qui
  28. *                    revient à jouer sur la viscosité efficace de l'émulsion
  29. *                    MUEFF qui tend vers la valeur de la phase pure 
  30. *                    reconstituée lorsque la fraction volumique dépasse
  31. *                    un certain seuil. Le présent cas-test, pour des raisons
  32. *                    d'obtention de résultats analytiques, laisse la
  33. *                    viscosité efficace constante quelque soit la valeur de
  34. *                    AD. Le front de sédimentation est alors très raide
  35. *                    et notre problème revient regarder évoluer le front
  36. *                    de sédimentation de billes tombant toujours à la
  37. *                    même vitesse puis s'arrêtant brutalement lorsqu'elles
  38. *                    rencontrent la zone sédimentée.
  39. *
  40. *
  41. *                    Le maillage fait une maille d'épaisseur pour accélérer
  42. *                    la simulation dans le cadre des tests du répertoire
  43. *                    castem/dgibi.
  44. *
  45. *
  46. *         Auteur     : Gilles Bernard-Michel
  47. *
  48. *         Révision   : version 1
  49. *
  50. *         date       : 21/12/99
  51. *
  52. ************************ Options générales *********************************
  53. *
  54. *
  55. ** Sortie graphique.
  56. *
  57.  
  58. GRAPH = 'O';
  59.  
  60. *
  61. *test long ou court.
  62. *
  63.  
  64. COMPLET = FAUX;
  65.  
  66.  
  67. *
  68. *période des sorties graphiques : tous les NITMA/PERIO pas de temps
  69. *
  70.  
  71. PERIO = 1 ;
  72.  
  73. *
  74. * nombre de mailles, pas de temps , intervalle de temps.
  75. *
  76.  
  77. SI (NON COMPLET) ;
  78.    nptx = 1. ;
  79.    npty = 10. ;
  80.    DELTATPS = 2.;
  81. *  nombre pas de temps
  82.    NITMA = 10 ;
  83. *  nombre points fixes
  84.    NITER = 3 ;
  85. *  convergence point fixe.
  86.    PREC = 1.D-4 ;
  87. SINON;
  88.    nptx = 1;
  89.    npty = 30;
  90.    DELTATPS = 0.5D0;
  91. *  nombre pas de temps
  92.    NITMA = 40 ;
  93. *  nombre points fixes
  94.    NITER = 5 ;
  95. *  convergence point fixe.
  96.    PREC = 1.D-4 ;
  97. FINSI;
  98.  
  99. *
  100. ************************ On définit les propriétés physiques ***************
  101. *
  102.  
  103. *
  104. ** viscosité dynamique mu/rho.
  105. *
  106.  
  107. nu =1.E-6 ;
  108.  
  109. *
  110. ********************* On traite les aspects géométriques ********************
  111. *
  112.  
  113. *
  114. ** Les paramètres de la géométrie.
  115. *
  116.  
  117. rint     = 0. ;
  118. rext     = 0.1;
  119. hh       = 0.5 ;
  120.  
  121. *
  122. ** Densité du maillage.
  123. *
  124.  
  125.  
  126. dx = (rext - rint) / nptx;
  127. dy = hh / npty ;
  128.  
  129. *
  130. ** Les sommets du bol
  131. *
  132.  
  133. p1 = rint 0. ;
  134. p2 = rext 0. ;
  135. p3 = rext hh ;
  136. p4 = rint hh ;
  137.  
  138. p6 = 1.0*dx 0.0;
  139. p8 = 1.D0*dx hh;
  140. *
  141. ** Assemblage des bords du domaine.
  142. *
  143.  
  144.  bas  = p1 'DROIT' dini dx dfin dx p2 ;
  145.  cotd = p2 'DROIT' dini dy dfin dy p3 ;
  146.  haut = p3 'DROIT' dini dx dfin dx p4 ;
  147.  cotg = p4 'DROIT' dini dy dfin dy p1 ;
  148.  
  149.  cnt  = bas 'ET' cotd 'ET' haut 'ET' cotg ;
  150.  
  151. *
  152. ** On maille avec des rectangles.
  153. *
  154.  
  155. mt= bas cotd haut cotg 'DALLER' ;
  156. 'TASSER' mt ;
  157.  
  158. *
  159. ** On réoriente les éléments.
  160. *
  161.  
  162. mt = 'ORIENTER' mt ;
  163.  
  164. &mt   = CHANGER mt   QUAF ;
  165. &bas  = CHANGER bas  QUAF ;
  166. &cotd = CHANGER cotd QUAF ;
  167. &haut = CHANGER haut QUAF ;
  168. &cotg = CHANGER cotg QUAF;
  169. &cnt  = CHANGER cnt  QUAF ;
  170.  
  171. 'ELIMINATION' 1.e-6 (&mt 'ET' &cnt 'ET' &cotg 'ET' &haut
  172.                  'ET' &cotg 'ET' &cnt);
  173.  
  174. $mt    = MODE &mt    'NAVIER_STOKES' QUAF ;
  175.  
  176. *
  177. ** On crée le domaine associé au maillage.
  178. *
  179.  
  180. mt   = doma $mt    maillage ;
  181.  
  182. * On crée une droite pour tracer des profils de fraction volumique.
  183. lmt = 'CHANGE' mt LIGNE ;
  184.  coup = lmt 'POINT' 'DROITE' p8 p6 1.D-5 ;
  185. alpdco = lmt 'ELEM' 'APPUYE' 'STRICTEMENT' coup ;
  186.  
  187. *
  188. ** on crée une frontière moins un noeud pour le cas
  189. ** de cavités carrés. Il ne faut pas retirer un coin
  190. ** car on peut avoir un champ sortant.
  191. *
  192.  
  193. $cnt =  'MODELISER' &cnt 'NAVIER_STOKES' QUAF;
  194. cntcc = 'DOMA' $cnt maillage ;
  195. cntcc = 'CHANGER' cntcc POI1;
  196.  
  197. retir = p8 'MANUEL' POI1;
  198. cntcc = retir 'DIFF' cntcc;
  199.  
  200. * *********************************
  201. * PROCEDURE CALCULANT LE MODULO 
  202. * *********************************
  203. *
  204. DEBPROC MODULO i*ENTIER j*ENTIER ;
  205. k*ENTIER = i / j ;
  206. mod = i - ( k * j ) ;
  207. FINPROC mod ;
  208.  
  209. *------------------------------------------------------------------------
  210. *-----------------------------------------------------------------------
  211. *                            - LIREXY3 -                               -
  212. * PROCEDURE POUR OBTENIR UNE LISTE CONTENANT LES VALEURS DU CHPO AUX   -
  213. * POINTS DE DONNEES : premiere étape pour le tracé d'histogramme       -
  214. * PROCEDURE EN TOUT POINT COMPARABLE A LIREXY SAUF QUE L'ON A UNE      -
  215. * TABLE EN SORTIE : POUR TRACE D'HISTOGRAMMES.                         -
  216. * (VOIR REMARQUE SUR POINT PROCHE DANS LIREXY).                        -
  217. *-----------------------------------------------------------------------
  218. * En entrée ...                                                              -
  219. *   - CHPOD : champ point dont on extrait les valeurs,                 -
  220. *   - MAILD : maillage des points centres où prendre les valeurs       -
  221. *   - K     : le nom (MOT) de la composante du champ que l'on extrait. -
  222. *-----------------------------------------------------------------------
  223. * En sortie ...                                                        -
  224. *   la liste des valeurs f(xi,yi)                                      -
  225. *-----------------------------------------------------------------------
  226. 'DEBPROC' LIREXY3 CHPOD*'CHPOINT' MAILD*'MAILLAGE' K*'MOT';
  227.   NPT   = NBEL MAILD;
  228.  
  229.   REPETER BACQ NPT;
  230.     I = &bacq;
  231.     Pt = MAILD POINT I;
  232.     Pty = 'COOR' 2 Pt ;
  233.     KI = 'EXTRAIRE'  CHPOD Pt K;
  234.     SI (I EGA 1);
  235.       LIPK = 'PROG'  KI;
  236.       LIPP = 'PROG'  Pty;
  237.     SINON;
  238.       LIPK = LIPK 'ET' ('PROG' KI);
  239.       LIPP = LIPP 'ET' ('PROG' Pty);
  240.     FINSI;
  241.   FIN BACQ;
  242.   evo = 'EVOL' 'MANU' LIPP LIPK;
  243.   dessin evo;
  244. 'FINPROC';
  245.  
  246. * Idem pour points sommet.
  247.  
  248. 'DEBPROC' LIREXY2 CHPOD*'CHPOINT' MAILD*'MAILLAGE' K*'MOT';
  249.   NPT   = NBEL MAILD;
  250.   NPT   = NPT + 1 ;
  251.  
  252.   REPETER BACQ NPT;
  253.     I = &bacq;
  254.     Pt = MAILD POINT I;
  255.     Pty = 'COOR' 2 Pt ;
  256.     KI = 'EXTRAIRE'  CHPOD Pt K;
  257.     SI (I EGA 1);
  258.       LIPK = 'PROG'  KI;
  259.       LIPP = 'PROG'  Pty;
  260.     SINON;
  261.       LIPK = LIPK 'ET' ('PROG' KI);
  262.       LIPP = LIPP 'ET' ('PROG' Pty);
  263.     FINSI;
  264.   FIN BACQ;
  265.   evo = 'EVOL' 'MANU' LIPP LIPK;
  266.   dessin evo;
  267. 'FINPROC';
  268.  
  269. *
  270. * ***************************************
  271. * PROCEDURE DE CALCUL DES COEFFICIENTS
  272. * DES EQUATIONS DE QUANTITE DE MOUVEMENT       
  273. * ***************************************
  274. *
  275. DEBPROC ACTU  RV*'TABLE' $mt*'MMODEL' ;
  276.  
  277. * constantes du probleme
  278. roc     = RV.INCO.'ROC' ;
  279. rod     = RV.INCO.'ROD' ;
  280. muc     = RV.INCO.'MUC' ;
  281. mud     = RV.INCO.'MUD' ;
  282. expo    = RV.INCO.'EXPOS' ;
  283. farcc   = RV.INCO.'FARCC' ;
  284. farcd   = RV.INCO.'FARCD' ;
  285.  
  286. * LECTURE DES INCONNUES *
  287.  
  288. RJC = RV.INCO.'RJC' ;
  289. RJD = RV.INCO.'RJD' ;
  290. JC = RV.INCO.'JC' ;
  291. JD = RV.INCO.'JD' ;
  292.  
  293. * implicitation en alpha pour les frottements.
  294. ALPM  = RV.INCO.'AM' ;
  295. *ALPD = 'COPIER' RV.INCO.'ALPD';
  296. ALPD = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
  297. RV.INCO.'AC' = kops 1.D0 '-' RV.INCO.'AD';
  298.  
  299. MASK1D = ALPD MASQUE 'INFERIEUR' (1.D0 '-' 1.D-15)  ;
  300. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D  ;
  301. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  302.           '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
  303. MASK1D = ALPD MASQUE 'SUPERIEUR' 1D-15  ;
  304. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  305. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  306.       '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
  307.  
  308. ALPC =  'KOPS' 1.D0 '-' RV.INCO.'AD';
  309.  
  310.  
  311. RAY = RV.INCO.'RAY' ;
  312. RAYC = RV.INCO.'RAYC' ;
  313. ALPI = RV.INCO.'RALPI' ;
  314. COAL = RV.INCO.'COAL' ;
  315.  
  316. *** COUPLAGE DYNAMIQUE ***
  317.  
  318. * VISCOSITE DYNAMIQUE EFFICACE
  319. MUEFF = ALPD  ;
  320. MUEFF = 'KCHT' $mt 'SCAL' 'SOMMET' MUEFF ;
  321. MUEFF = 'KOPS' 1.0D0 '-' MUEFF ;
  322. MUEFF = 'KOPS' MUEFF '**' expo ;
  323. MUEFF = 'KOPS' MUC '*' MUEFF ;
  324.  
  325. RV.INCO.'MUEFF' = 'KCHT' $mt SCAL 'SOMMET' MUEFF ;
  326. RV.INCO.'SNUFFC' = 'KOPS' MUEFF '/' ROC ;
  327. RV.INCO.'SNUFFC' = 'KOPS' ALPC * RV.INCO.'SNUFFC';
  328. RV.INCO.'NUEFFC' = 'KCHT' $mt scal centre
  329.              ('NOEL' $mt RV.INCO.'SNUFFC');
  330. RV.INCO.'SNUFFD' = 'KOPS' MUEFF '/' ROD ; 
  331. RV.INCO.'SNUFFD' = 'KOPS' ALPD * RV.INCO.'SNUFFD';
  332. RV.INCO.'NUEFFD' = 'KCHT' $mt scal centre
  333.              ('NOEL' $mt RV.INCO.'SNUFFD');
  334.  
  335.  
  336. ** UC ET UD : VITESSES INSTANTANEES relaxées **
  337. SEUIL     = (('TANH' (1.D20 *
  338.                  (ALPC '-' 1.D-3))) '+' 1)/2.D0;
  339. ALPCX     = NOMC 'UX'    ALPC                  ;
  340. ALPCY     = NOMC 'UY'    ALPC                  ;
  341. ALPCV     = ALPCX 'ET' ALPCY                    ;
  342. UC = 'KOPS' RJC '/' (ALPCV+1.e-20)             ;
  343. UC = 'KOPS' SEUIL '*' UC                       ;
  344. RV.INCO.'UC'  = 'KCHT' $mt VECT 'SOMMET' UC     ;
  345.  
  346. SEUIL     = (('TANH' (1.D20 *
  347.             (ALPD '-' 1.D-3))) '+' 1)/2.D0; 
  348. ALPDX     = NOMC 'UX'     ALPD                 ;
  349. ALPDY     = NOMC 'UY'     ALPD                 ;
  350. ALPDV     = ALPDX 'ET' ALPDY                    ;
  351. UD = 'KOPS' RJD '/' ALPDV                       ;
  352. UD = 'KOPS' SEUIL '*' UD                       ;
  353. RV.INCO.'UD'  = 'KCHT' $mt VECT 'SOMMET' UD     ;
  354.  
  355. ** VUD : VITESSES INSTANTANEES lissées**
  356. VUD = 'NOEL' $mt UD                            ;
  357. VUD = 'KCHT' $mt VECT 'CENTRE' VUD             ;
  358. VUD = 'ELNO' $mt VUD                           ;
  359. RV.INCO.'VUD'  = 'KCHT' $mt VECT 'SOMMET' VUD  ;
  360.  
  361.  
  362. * COEFF DE FROTTEMENT
  363. *on calcule une fraction volumique ALPD prenant
  364. *en compte les parois verticales = max (alpd alpdm).
  365. maskp  = ALPD 'MASQUE' 'SUPERIEUR' ALPM ;
  366. maskp  = 'KCHT' $mt scal sommet maskp;
  367. maskpi = 'KOPS' 1.D0 '-' maskp ;
  368. ALPDF  = 'KOPS' maskp * ALPD ;
  369. ALPFF  = 'KOPS' maskpi * ALPM ;
  370. ALPDF  = 'KOPS' ALPDF '+' ALPFF ;
  371.  
  372. *HINDER =  ALPDF ;
  373. HINDER  =  0.5*((KOPS 0.01D0 '+' ALPC) 
  374.            + ( abs (KOPS 0.01D0 '-' ALPC)));
  375. *HINDER = ALPC;
  376. HINDER = 'KCHT' $mt 'SCAL' 'SOMMET' HINDER ;
  377. HINDER = 'KOPS' HINDER '**' expo ;
  378. HINDER = 'KOPS' MUC '*' HINDER ;
  379. kft = 'KCHT' $mt SCAL 'SOMMET' (9.D0/2.D0) ;
  380. kft = 'KOPS' kft '*' ('KOPS' RV.INCO.'RAY' '**' -2) ;
  381.  
  382. RV.INCO.'FRCD' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' -1.D0 
  383.     '*'  ('KOPS' kft '/' roc) ) ;    
  384. RV.INCO.'FRCD1' = 'KOPS' RV.INCO.'FRCD' '*' HINDER ;
  385. RV.INCO.'FRCD2' = -1.0D0 * ('KOPS' RV.INCO.'FRCD1' '*'
  386.                  RV.INCO.'AD') ;
  387. RV.INCO.'FRCD' = 'KOPS' RV.INCO.'FRCD1' '+' RV.INCO.'FRCD2' ;
  388. RV.INCO.'MFRCD' = -0.5D0 * ('KOPS' RV.INCO.'FRCD1' '*' RV.INCO.'JD') ;
  389.  
  390. RV.INCO.'FRDC' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' -1.D0
  391.     '*' ('KOPS' kft '/' rod)) ;
  392. RV.INCO.'FRDC' = 'KOPS' RV.INCO.'FRDC' '*' HINDER ; 
  393. RV.INCO.'MFRDC' = 0.5D0 * ('KOPS' RV.INCO.'FRDC' '*' RV.INCO.'JC') ;
  394. RV.INCO.'FRDC' = 1.0D0 * ('KOPS' RV.INCO.'FRDC' '*'
  395.                   RV.INCO.'AD') ;
  396.  
  397. RV.INCO.'FRCC' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' 1.D0 
  398.     '*'  ('KOPS' kft '/' roc) ) ;    
  399. RV.INCO.'FRCC' = 'KOPS' RV.INCO.'FRCC' '*' HINDER ;
  400. RV.INCO.'MFRCC' = 0.5D0 * ('KOPS' RV.INCO.'FRCC' '*' RV.INCO.'JC') ;
  401. RV.INCO.'FRCC' = 1.0D0 * ('KOPS' RV.INCO.'FRCC' '*'
  402.                        RV.INCO.'AD') ;
  403.  
  404. RV.INCO.'FRDD' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' 1.D0
  405.     '*' ('KOPS' kft '/' rod)) ;
  406. RV.INCO.'FRDD1' = 'KOPS' RV.INCO.'FRDD' '*' HINDER ; 
  407. RV.INCO.'FRDD2' = -1.0D0 * ('KOPS' RV.INCO.'FRDD1' '*'
  408.                    RV.INCO.'AD') ;
  409. RV.INCO.'FRDD' = 'KOPS' RV.INCO.'FRDD1' '+' RV.INCO.'FRDD2' ;
  410. RV.INCO.'MFRDD' = -0.5D0 * ('KOPS' RV.INCO.'FRDD1' '*' RV.INCO.'JD') ;
  411.  
  412. RV.INCO.'KFT' = 'KCHT' $mt SCAL SOMMET kft;
  413.  
  414. * COEFFICIENT DE COALESCENCE
  415.  DIST = 'KOPS' ALPD '**' (-1.D0/3.D0);
  416.  DIST = 'KCHT' $mt SCAL 'SOMMET' DIST ;
  417.  DIST = 'KOPS' DIST '*' 2.03D0 ;
  418.  DIST = 'KOPS' ((3.D0 ** 0.5D0)/2.D0) '*' DIST ;
  419.  DIST = 'KOPS' DIST '-' 2.D0 ;
  420.  
  421.  MDIST = DIST MASQUE 'SUPERIEURE' 1.D-4 ;
  422.  MDIST = 'KCHT' $mt SCAL 'SOMMET' MDIST ;
  423.  
  424.  DIST = 'KOPS' DIST '*' MDIST ;
  425.  DIST = 'KOPS' DIST '+' ('KOPS'
  426.         ('KOPS' 1.D0 '-' MDIST) '*' 1.D-4) ;
  427.  
  428.  COAL = 'KOPS' ALPD '/' (ALPC+1.e-20) ;
  429.  COAL = 'KOPS' COAL '*' ECSTA2 ;
  430.  
  431.  COALC = 'KOPS' ALPC '/' ALPD ;
  432.  COALC = 'KCHT' $mt SCAL 'SOMMET' COALC ;
  433.  COALC = 'KOPS' COALC '*' 1.D-2 ;
  434.  COALC = 'KOPS' ('KOPS' 1.D0 - ECSTA2) '*' COALC ;
  435.  
  436.  COAL = 'KOPS' COALC '+' COAL ;
  437.  COAL = 'KOPS' COAL '*' C0 ;
  438.  
  439. RV.INCO.'COAL' = 'KCHT' $mt SCAL 'SOMMET' COAL ;
  440. RV.INCO.'C0' = 'KCHT' $mt SCAL 'SOMMET' C0 ; 
  441. RV.INCO.'DIST' = 'KCHT' $mt SCAL 'SOMMET' DIST ; 
  442.  
  443. * termes de pression 
  444.  
  445. PUC = 'KOPS' -1.D0 '*' RV.INCO.'AC'  ;
  446. PUC = 'KOPS' PUC '/' ROC ;
  447. RV.INCO.'PUC' = kcht $mt SCAL SOMMET PUC ;
  448.  
  449. PUD = 'KOPS' -1.D0 '*' RV.INCO.'AD'  ;
  450. PUD = 'KOPS' PUD '/' ROD ;
  451. RV.INCO.'PUD' = kcht $mt SCAL SOMMET PUD ;
  452. RV.INCO.'PUD' = 1.0 * RV.INCO.'PUD';
  453.  
  454. PNN = ELNO $mt RV.INCO.PRES CENTREP1;
  455. PNN = KCHT $mt SCAL SOMMET PNN;
  456. GPNN =  KOPS GRAD PNN $MT;
  457. GPNN = KCHT $mt VECT CENTRE GPNN;
  458. GPNN = elno $mt GPNN;
  459. RV.INCO.'GPUD' = kcht $mt VECT SOMMET (0.5D0 * (1/ROC) * GPNN);
  460.  
  461. * terme source pesanteur signe + si TOIMP
  462. * signe - si MDIA ou dans DFDT. on met de l'ordre
  463. * 2, au lieu de alpc alpd(n+1), on a
  464. * (alpc-alpd)alpd(n+1) + alpd*alpd
  465.  
  466. farchc = 'KOPS' ('KOPS' RV.INCO.'AC' '-' RV.INCO.'AD') '*' farcc;
  467. farchci = -1.D0 * farchc ;
  468. *partie expl
  469. farchce = 'KOPS' ('KOPS' RV.INCO.'AD' * RV.INCO.'AD') '*' farcc;
  470. farchce = -1.D0 * farchce;
  471. RV.INCO.'FARCHCI' =  'KCHT' $mt VECT 'SOMMET' farchci ;
  472. RV.INCO.'FARCHCE' =  'KCHT' $mt VECT 'SOMMET' farchce ;
  473.  
  474. farchd = 'KOPS' ('KOPS' RV.INCO.'AC' '-' RV.INCO.'AD') '*' farcd;
  475. farchdi = -1.D0 * farchd ;
  476. *partie expl
  477. farchde = 'KOPS' ('KOPS' RV.INCO.'AD' * RV.INCO.'AD') '*' farcd;
  478. farchde = -1.D0 * farchde;
  479. RV.INCO.'FARCHDI' =  'KCHT' $mt VECT 'SOMMET' farchdi ;
  480. RV.INCO.'FARCHDE' =  'KCHT' $mt VECT 'SOMMET' farchde ;
  481.  
  482. * equation du rayon
  483.  
  484. RV.INCO.'RAYALPI' = (1.D0/3.D0)* RV.INCO.'ALPI';
  485.  
  486. * terme de stabilisation QDM
  487. RV.INCO.'GRALPC' = kops GRAD RV.INCO.'AD' $mt;
  488. RV.INCO.'GRALPC' = kcht $mt VECT CENTRE RV.INCO.'GRALPC';
  489. RV.INCO.'GRALPC' = KOPS RV.INCO.'DR2' '*' RV.INCO.'GRALPC';
  490. RV.INCO.'GRALPC' = 'ABS' RV.INCO.'GRALPC';
  491. RV.INCO.'GRALPC' = .1D0 *  RV.INCO.'GRALPC';
  492. RV.INCO.'GRALPS' = kcht $mt VECT SOMMET (ELNO $mt RV.INCO.'GRALPC');
  493.  
  494. RV.INCO.'GRALJD' = 'KCHT' $mt scal sommet ('EXCO' uy RV.INCO.'JD');
  495. RV.INCO.'GRALJD' = kops GRAD RV.INCO.'GRALJD' $mt;
  496. RV.INCO.'GRALJD' = kcht $mt VECT CENTRE RV.INCO.'GRALJD';
  497. RV.INCO.'GRALJD' = KOPS RV.INCO.'DR2' '*' RV.INCO.'GRALJD';
  498. RV.INCO.'GRALJDM' = 0.1D0 *  RV.INCO.'GRALJD';
  499. RV.INCO.'GRALJD' = -10.D0 *  RV.INCO.'GRALJD';
  500. RV.INCO.'GRALJS' = kcht $mt VECT SOMMET (ELNO $mt RV.INCO.'GRALJD');
  501.  
  502. FINPROC RV;
  503.  
  504. *************************
  505. *** DONNEES PHYSIQUES ***
  506. *************************
  507. * formule Ishii
  508. *expo= -4.1;
  509. * ON NE FAIT PAS VARIER MUEFF VOIR DESCRIPTION DU MODELE
  510. expo = -0.0D0;
  511. *Densité volumique
  512. ROC = 1000.D0   ;
  513. ROD = 1100.D0   ;
  514. epsilc = (ROC - ROD) / ROC ;
  515. epsild = (ROD - ROC) / ROD ;
  516. * Viscosite cinématique et dynamique de la phase liquide
  517. MUC  = 0.001D0 ;
  518. NUC  = MUC / ROC ;
  519. MUD  = 0.001D0 ;
  520. NUD  = MUD / ROD ;
  521. * Rayon initial des gouttes
  522. ray0 = 1.D-4 ;
  523. * Concentration initiale de la phase dispersée
  524. ALPD0 = 0.3D0 ;
  525. ALPI0 = 3.D0 * ALPD0 / ray0 ;
  526. * Concentration initiale de la phase continue
  527. ALPC0 = 1 - ALPD0 ;
  528. C0 = 1.D-2 ;
  529. * Coeff. de Pénalisation *
  530. EPSS=1.D-9;
  531. * Valeur minimale de EPSN *
  532. MIEPS = 1.e-5 ;
  533. * Valeur maximale de EPSN *
  534. MAEPS = 0.99999 ;
  535. * Vicosité maximale et truc maximal *
  536. ALPMAX = 0.9D0 ;
  537. MUMAX = 1.D-1 ;
  538. RAY2  = 1.D-3;
  539. * poussée d'Archimède multipliée par +1 car TOIMP à droite de l'égaZZlité
  540. gravos = 'KCHT' $mt VECT SOMMET (0.D0 -9.81D0) ;
  541. archi  =  gravos ;
  542. farcc = epsilc* archi ;
  543. farcd = epsild* archi ;
  544. * fraction vol max
  545. alphadm = 0.9900000001D0 ;
  546. alphacm = MIEPS ;
  547. * profil de fraction volumique
  548. cory = 'COORDONNEE' 2 &mt;
  549. alpd1 = 'TANH' (200.D0 * (cory '-' (3.*hh/10.D0))) ;
  550. alpd1 = (0.999 '-' ALPD0) * (1.D0 '-' alpd1) '/' 2.D0;
  551. alpd2 = 'TANH' (200.D0 * (cory '-' (9.D0 * hh/10.D0))) ;
  552. alpd2 = ALPD0 * (1.D0 '-' alpd2) '/' 2.D0;
  553. alpd = alpd1 '+' alpd2;
  554. alpini = alpd;
  555.  
  556. *'TRACER' ALPD &mt;
  557.  
  558. * PAS DE TEMPS *
  559. TERDT = DELTATPS ;
  560.  
  561. ************************************
  562. *** SYSTEMES D'EQUATIONS ET C.L. ***
  563. ************************************
  564. *
  565. *
  566. *************************************
  567. ***    SYSTEME PORTANT SUR :      ***
  568. ***          UN, VN, P            ***
  569. *************************************
  570.  
  571. * EQUATION DE CONSERVATION DE LA QUANTIE DE MVT PHASE CONT
  572. *
  573. * composantes r,z
  574. RV = 'EQEX' $mt ITMA 1  
  575. OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'CENTREE'
  576. ZONE $mt OPER 'DFDT' 1.D0 'JCP' 'JCPP' 'DT' 'UC' 1.D-15 INCO 'JC'
  577. ZONE $mt OPER 'MDIA'  'FARCHCI'        INCO 'AD' 'JC'
  578. ;
  579. RV = 'EQEX' RV 
  580. OPTI 'EFM1' 'IMPL' 'CENTREE'
  581. ZONE $mt OPER 'MDIA' 'FRCD' INCO 'JD' 'JC'
  582. ZONE $mt OPER 'MDIA' 'FRCC' INCO 'JC' 'JC'
  583. *OPTI 'EF' 'IMPL' 'CONS' 'CENTREE' 'MUCONS'
  584. OPTI 'EF' 'IMPL' 'CONS' 'CENTREE'
  585. *ZONE $mt OPER 'KONV' 1.D0 'UC' 1.D-15 DT INCO 'JC'
  586. ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JC'
  587. ;
  588. RV = EQEX  RV
  589. OPTI 'EF'  'IMPL' 'CENTREP1' 'CENTREE'
  590. ZONE $mt OPER 'KMBT'           'PUC'        INCO 'PRES' 'JC'
  591. ;
  592.  
  593. *
  594. * EQUATION DE CONSERVATION DE LA QUANTITE DE MVT PHASE DISP
  595. *
  596. * composantes r, z
  597. RV = EQEX RV  
  598. OPTI 'EFM1' 'IMPL' 'BDF2' 'CONS' 'CENTREE'
  599. ZONE $mt OPER 'DFDT' 1.D0 'JDP' 'JDPP' 'DT' 'UD' 1.D-15 INCO 'JD'
  600. ZONE $mt OPER 'MDIA'  'FARCHDI'        INCO 'AD' 'JD'
  601. ZONE $mt OPER 'MDIA' 'FRDC' INCO 'JC' 'JD'
  602. ZONE $mt OPER 'MDIA' 'FRDD' INCO 'JD' 'JD'
  603. OPTI 'EF' 'IMPL' 'CENTREE'
  604. *ZONE $mt OPER 'TOIMP'  'FARCHD'        INCO 'JD'
  605. ;
  606.  
  607. RV = EQEX RV
  608. *OPTI 'EF' 'IMPL'  'CENTREE'   'MUCONS'
  609. OPTI 'EF' 'IMPL'  'CENTREE'
  610. ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JC'
  611. ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JD'
  612. OPTI 'EF' 'IMPL' 'CONS' 'SUPG'
  613. *ZONE $mt OPER 'KONV' 1.D0 'UD' 1.D-15 DT INCO 'JD'
  614. ;
  615.  
  616. ;
  617.  
  618. RV = EQEX  RV
  619. OPTI 'EF'  'IMPL' 'CENTREP1' 'CENTREE'
  620. *ZONE $mt OPER 'KMBT'  'GRALPS' 0.D0 'SOMM' 1.D-15  INCO 'AD' 'JD'
  621. *ZONE $mt OPER 'KMBT'   -0.001  0.D0 'SOMM' 1.D-15  INCO 'AD' 'JD'
  622. ZONE $mt OPER 'KMBT'          'PUD'         INCO 'PRES' 'JD'
  623. ;
  624.  
  625. *
  626. * EQUATION DE CONSERVATION DU VOLUME (INCOMPRESSIBLE)
  627. *
  628.  
  629. RV = EQEX  RV
  630. OPTI 'EFM1'  'IMPL' 'CENTREP1' 'CENTREE' 
  631. ZONE $mt OPER 'KMAB'        1.D0  0.           INCO 'JC' 'PRES'
  632. ZONE $mt OPER 'KMAB'        1.D0  0.           INCO 'JD' 'PRES'
  633. *OPTI 'EFMC' 'IMPL' 'CENTREE' 'INCOD' CENTREP1
  634. *ZONE $mt OPER 'DFDT'  EPSS 'PRES' DT 'JNUL' 1.D+15 INCO 'PRES'
  635. ;
  636.  
  637.  
  638. * 
  639. * EQUATION DE CONSERVATION DE LA FRACTION VOLUMIQUE
  640. *
  641.  
  642. RV = EQEX RV  
  643. OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'CENTREE' 
  644. ZONE $mt OPER 'DFDT' 1.D0 'ADP' 'ADPP' 'DT' 'UD' 1.D-15 INCO 'AD'
  645. OPTI 'EF'  'IMPL' 'CENTREP1' 'CONS' 'SUPG' 'CMD' 0.20
  646. ZONE $mt OPER 'KMAB' 1.D0 0.D0 'SOMM' 1.D-15 INCO 'JD' 'AD'
  647. *OPTI 'EF' 'IMPL'  'CENTREE'
  648. *ZONE $mt OPER 'LAPM' 'GRALPC' INCO 'AD'
  649. ;
  650.  
  651.  
  652. *** FIN SECOND SYSTEME ***
  653. **************************
  654.  
  655. *
  656. * CONDITION AUX LIMITES
  657. *
  658.  
  659. RV = EQEX  RV
  660.   'CLIM' 'JC' 'UIMP' (&cotg 'ET' &cotd)           0.D0
  661.          'JC' 'VIMP' (&bas 'ET' &haut)            0.D0
  662.          'JD' 'UIMP' (&cotg 'ET' &cotd)           0.D0
  663.          'JD' 'VIMP' (&bas 'ET' &haut)            0.D0
  664.          'AD' 'TIMP' (&bas)      1.D0
  665.          'AD' 'TIMP' (&haut)     0.D0
  666.    ;
  667. bcp='ELEM' ('DOMA' $mt 'CENTREP1')  POI1 (lect 1) ;
  668. RV = EQEX  RV
  669.    'CLIM'
  670.           'PRES' 'TIMP' bcp 0.D0  ;
  671.  
  672.  
  673. *
  674. *******************************************
  675. * INITIALISATION DES TABLES DES INCONNUES
  676. *******************************************
  677. *
  678.  
  679. *
  680. RV.INCO        = TABLE INCO                       ;
  681.  
  682. RV.INCO.'ROC'   = ROC   ;
  683. RV.INCO.'ROD'   = ROD   ;
  684. RV.INCO.'MUC'   = MUC   ;
  685. RV.INCO.'MUD'   = MUD   ;
  686. RV.INCO.'EXPOS' = expo  ;
  687. RV.INCO.'FARCC' = farcc ;
  688. RV.INCO.'FARCD' = farcd ;
  689. RV.INCO.'ALPHADMAX' = alphadm ;
  690.  
  691. *
  692. RV.'DT'       = DELTATPS ;
  693. RV.INCO.'DT'  = RV.'DT' ;
  694. RV.INCO.'GDT' = -1.D0 * epsild * 9.81D0 * DELTATPS ;
  695.  
  696. RV.INCO.'JC'   = 'KCHT' $mt VECT 'SOMMET'  (1.e-25 1.e-25) ;
  697. RV.INCO.'JD'   = 'KCHT' $mt VECT 'SOMMET'  (1.e-25 1.e-25) ;
  698. RV.INCO.'JCP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC'    ;
  699. RV.INCO.'JCPP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC'    ;
  700. RV.INCO.'JDP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD'    ;
  701. RV.INCO.'JDPP'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD'    ;
  702.  
  703. RV.INCO.'RJC'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC'    ;
  704. RV.INCO.'RJD'  = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD'    ;
  705.  
  706. RV.INCO.'PRES' = 'KCHT' $mt SCAL 'CENTREP1'   1.e-5         ;
  707.  
  708.  
  709. RV.INCO.'AC' = 'KCHT' $mt SCAL 'SOMMET'   (1.D0 '-' ALPD)  ;
  710. RV.INCO.'AD' = 'KCHT' $mt SCAL 'SOMMET'   ALPD         ;
  711. RV.INCO.'RALPD' = 'KCHT' $mt SCAL 'SOMMET'   ALPD        ;
  712. RV.INCO.'RALPC' = 'KCHT' $mt SCAL 'SOMMET'   (1.D0 '-' ALPD) ;
  713. RV.INCO.'ADP'=  'KCHT' $mt SCAL 'SOMMET' RV.INCO.'AD'    ;
  714. RV.INCO.'ADPP'= 'KCHT' $mt SCAL 'SOMMET' RV.INCO.'AD'    ;
  715. RV.INCO.'GRALPS' = KCHT $mt VECT CENTRE (0.D0 0.D0)         ;
  716. RV.INCO.'GRALPC' = KCHT $mt VECT CENTRE (0.D0 0.D0)         ;
  717.  
  718. RV.INCO.'RAY'   = 'KCHT' $mt SCAL 'SOMMET'  ray0            ; 
  719.  
  720. RV.INCO.'ALPI'  = 'KCHT' $mt SCAL 'SOMMET' ALPI0          ;
  721. RV.INCO.'ALPIP' = 'KCHT' $mt SCAL 'SOMMET'  RV.INCO.'ALPI'    ;          
  722. RV.INCO.'RALPI' = 'KCHT' $mt SCAL 'SOMMET' ALPI0          ;
  723. RV.INCO.'RAYALPI' = KCHT $mt SCAL 'SOMMET' (3.D0 * RV.INCO.'ALPI');
  724.  
  725. RV.INCO.'FARCHCI' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
  726. RV.INCO.'FARCHDI' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
  727. RV.INCO.'FARCHCE' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
  728. RV.INCO.'FARCHDE' = 'KCHT' $mt vect SOMMET (0.D0 0.D0)       ;
  729.  
  730. * INITIALISATION DES COEFFICIENTS
  731.  
  732. RV.INCO.'DR2' = kcht $mt VECT 'CENTRE' ((dx*dx) (dy*dy));
  733. *
  734. RV.INCO.'CNUL' = 'KCHT' $mt SCAL 'SOMMET' 0.D0            ;
  735. RV.INCO.'JNUL' = 'KCHT' $mt VECT 'SOMMET' (1.D-25 1.D-25) ;
  736. RV.INCO.'JUNI' = 'KCHT' $mt VECT 'SOMMET' (1.D0 1.D0)     ; 
  737.  
  738. RV.INCO.'CNUL1' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
  739. RV.INCO.'CNUL2' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
  740. RV.INCO.'CNUL3' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
  741. RV.INCO.'CNUL4' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)    ;
  742. RV.INCO.'SOUC' = 'KCHT' $mt VECT 'CENTRE' (1.e-25 1.e-25) ;
  743. RV.INCO.'SOUD' = 'KCHT' $mt VECT 'CENTRE' (1.e-25 1.e-25) ;
  744.  
  745. RV.INCO.'PUC'  = 'KCHT' $mt SCAL 'SOMMET'  1.e-5          ;
  746. RV.INCO.'PUD'  = 'KCHT' $mt SCAL 'SOMMET'  1.e-5          ;
  747. RV.INCO.'GPUD'  = 'KCHT' $mt VECT 'SOMMET'  (0.D0 0.D0)   ;
  748.  
  749. RV.INCO.'UC'   = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
  750. RV.INCO.'UD'   = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
  751.  
  752. RV.INCO.'VUD'   = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25);
  753.  
  754. RV.INCO.'NUEFFC' = 'KCHT' $mt SCAL 'CENTRE' 0.0D0         ;
  755. RV.INCO.'NUEFFD' = 'KCHT' $mt SCAL 'CENTRE' 0.0D0         ;
  756. RV.INCO.'SNUFFC' = 'KCHT' $mt SCAL 'SOMMET' 0.0D0         ;
  757. RV.INCO.'SNUFFD' = 'KCHT' $mt SCAL 'SOMMET' 0.0D0         ;
  758.  
  759. RV.INCO.'MUDIFC' = 'KCHT' $mt SCAL 'CENTRE' 0.D0          ;
  760. RV.INCO.'MUDIFD' = 'KCHT' $mt SCAL 'CENTRE' 0.D0          ;
  761.  
  762. RV.INCO.'FRCC'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  763. RV.INCO.'FRCD1'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  764. RV.INCO.'FRCD2'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  765. RV.INCO.'FRCD'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  766. RV.INCO.'FRDC'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  767. RV.INCO.'FRDD1'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  768. RV.INCO.'FRDD2'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  769. RV.INCO.'FRDD'   = 'KCHT' $mt SCAL 'SOMMET' 0.D0          ;
  770. RV.INCO.'MFRCC'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
  771. RV.INCO.'MFRCD'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
  772. RV.INCO.'MFRDC'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
  773. RV.INCO.'MFRDD'   = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0)         ;
  774.  
  775. RV.INCO.'COAL'   = 'KCHT' $mt SCAL 'SOMMET' C0            ;
  776.  
  777. *COEFFICIENT DE RELAXATION
  778. ALFI = 1.0D0 ;
  779. ALFID = 1.0D0 ;
  780. ALFII = 1.D0 ;
  781.  
  782. * Compteur des pas de temps
  783. ntemp = 0 ;
  784.  
  785. *** Calcul de ALPM frac vol pour la paroi ***********************
  786. *** distance à la paroi la plus proche **************************
  787.  
  788. corx = 'COOR' &mt 1;
  789. corx = 'KCHT' $mt scal sommet corx;
  790. corxi = 'KOPS' rext '-' corx;
  791. corx = 'KOPS'  corx '-' rint; 
  792. masko  = corx 'MASQUE' 'INFERIEUR' corxi ;
  793. masko = 'KCHT' $mt scal sommet masko;
  794. maskoi = 'KOPS' 1.D0 '-' masko;
  795. corx   = 'KOPS' masko * corx ;
  796. corxi  = 'KOPS' maskoi * corxi ;
  797. corx  = 'KOPS' corx '+' corxi ;
  798.  
  799. ratio = 'KOPS' corx '/' RV.INCO.'RAY' ;
  800. ratio = 'KOPS' ratio '+' 3.D0 ;
  801. drel  = 2.D0 * (alphadm**0.34);
  802. ratio = 'KOPS' ratio '/' drel;
  803. * on fabrique une fraction volumique telle que les écarts
  804. * entre particules soient identiques à la distance à la paroi
  805. * la plus proche.
  806. RV.INCO.'AM' = 'KOPS' ratio ** (-3);
  807.  
  808. *'TRACER' RV.INCO.'AM' &mt ('CONTOUR' &mt);
  809.  
  810.  
  811. *
  812. **** BOUCLES DES ITERATIONS INTERNES ET EN TEMPS ****************
  813. *
  814. REPETER BLOCKT NITMA ;
  815. *
  816. ntemp = ntemp + 1 ;
  817. 'MESSAGE' 'Itération externe N°' ntemp;
  818.  
  819. ntempi = 0 ;
  820.  
  821. * On stocke les variables au pas de temps ntemp
  822. VJCP =  'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JCP' ;
  823. VJCPP =  'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JCPP' ;
  824. VJDP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JDP' ;
  825. VJDPP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JDPP' ;
  826. VALPDP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADP' ;
  827. VALPDPP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADPP';
  828. VALPIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ALPI' ;
  829.  
  830.  
  831. * On effectue les itérations de point fixe
  832. *
  833. 'REPETER' BLOCKI NITER ;
  834. *
  835. ntempi = ntempi + 1 ;
  836. *
  837. *
  838. *************************************
  839. ***    SYSTEME PORTANT SUR :      ***
  840. ***          UN, VN, P            ***
  841. *************************************
  842. *
  843. * On charge les résultats du vrai pas de temps précédent
  844. RV.INCO.'JCP' = 'KCHT' $mt 'VECT' SOMMET VJCP;
  845. RV.INCO.'JCPP' = 'KCHT' $mt 'VECT' SOMMET VJCPP;
  846. RV.INCO.'JDP' = 'KCHT' $mt 'VECT' SOMMET VJDP;
  847. RV.INCO.'JDPP' = 'KCHT' $mt 'VECT' SOMMET VJDPP;
  848. RV.INCO.'ADP' = 'KCHT' $mt SCAL SOMMET VALPDP ;
  849. RV.INCO.'ADPP' = 'KCHT' $mt SCAL SOMMET VALPDPP ;
  850. RV.INCO.'ALPIP' = 'KCHT' $mt SCAL SOMMET VALPIP ;
  851. * resultats de l'iter interne préc.
  852. JCIP = 'KCHT' $mt VECT SOMMET RV.INCO.'JC' ;
  853. JDIP = 'KCHT' $mt 'VECT' SOMMET RV.INCO.'JD' ;
  854. ALPDIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AD' ;
  855. ALPIIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ALPI' ;
  856. ALPCIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AC' ;
  857. * Rayon moyen
  858. ECSTA2 = RV.INCO.'RAY' MASQUE 'INFERIEURE' RAY2 ;
  859. ECSTA2 = 'KCHT' $mt SCAL 'SOMMET' ECSTA2 ;
  860.  
  861. *RAY = 'KOPS' ALPDIP '/' ALPIIP ;
  862. *RAY = 'KCHT' $mt SCAL 'SOMMET' RAY ;
  863. *RAY = 'KOPS' RAY '*' ECSTA2 ;
  864. *RAY = 'KOPS' RAY '*' 3 ;
  865.  
  866. RAYC = 'KOPS' ('KOPS' 1.D0 '-' ALPCIP) '/' ALPIIP ;
  867. RAYC = 'KCHT' $mt SCAL 'SOMMET' RAYC ;
  868. RAYC = 'KOPS' RAYC '*' 3 ;
  869. RV.INCO.'RAYC' = RAYC ;
  870.  
  871.  
  872. *
  873. * On réactualise les données.
  874. *
  875.  
  876.  
  877.  
  878. RV = ACTU RV $mt;
  879.  
  880.  
  881. * On affecte le terme source expl de QDM par le biais de DFDT
  882. * coef 2 car ordre 2 en temps dans dfdt. JCPP sinon la poussee
  883. * archi est rechargéee * fois lors des reinitialisations
  884. RV.INCO.'JCPP'  = 'KOPS' RV.INCO.'JCP' '+'
  885.                  ('KOPS' (2.D0*RV.INCO.'DT') * RV.INCO.'FARCHCE') ;
  886.  
  887. RV.INCO.'JDPP'  = 'KOPS' RV.INCO.'JDP' '+'
  888.                  ('KOPS' (2.D0*RV.INCO.'DT') * RV.INCO.'FARCHDE') ;
  889.  
  890.  
  891. EXEC RV ;
  892. 'MESSAGE' 'Boucle interne N°' ntempi; 
  893.  
  894.  
  895. * CALCUL DE L'ERREUR RELATIVE MAX. *
  896. ERR = 0.D0 ;
  897. ERRC = 'KOPS' ('KOPS' RV.INCO.'JC' '-' JCIP)
  898.           '/' ('MAXIMUM' ('ABS' JCIP) '+' 1.D-10) ;
  899. ERRD = 'KOPS' ('KOPS' RV.INCO.'JD' '-' JDIP)
  900.           '/' ('MAXIMUM' ('ABS' JDIP) '+' 1.D-10) ;
  901. ERRC = 'MAXIMUM' ('ABS' ERRC) ;
  902. ERRD = 'MAXIMUM' ('ABS' ERRD) ;
  903. 'SI' (ERRC > ERR) ;
  904. ERR = ERRC ;
  905. VAR = 'MOTS' JC ;
  906. 'FINSI' ;
  907.  
  908. 'SI' (ERRD > ERR) ;
  909. ERR = ERRD ;
  910. VAR = 'MOTS' JD ;
  911. 'FINSI' ;
  912. ERRA = 'KOPS' ('KOPS' RV.INCO.'AD' '-' ALPDIP)
  913.           '/' ('MAXIMUM' ('ABS' ALPDIP) '+' 1.D-10) ;
  914. ERRA = 'MAXIMUM' ('ABS' ERRA) ;
  915. 'SI' (ERRA > ERR) ;
  916. ERR = ERRA ;
  917. VAR = 'MOTS' ALPHAD ;
  918. 'FINSI' ;
  919.  
  920. * CALCUL DE L'ERREUR RELATIVE MAX. *
  921.  
  922.  
  923. *'LIST' VAR;
  924. 'MESSAGE' ERR; 
  925. 'MESSAGE' ERRC ERRD ERRA;
  926. 'SI' (ERR <EG PREC) ;
  927. QUITTER BLOCKI ;
  928. 'FINSI'       ;
  929.  
  930. *
  931. 'FIN' BLOCKI ;
  932. *
  933.  
  934. 'MESSAGE' 'NBR ITER INTERNES' ntempi;
  935. 'MESSAGE' 'ERREUR RELATIVE MAX' ERR ;
  936. 'MESSAGE' 'Max JC' ('MAXIMUM' (RV.INCO.'JC'))   ;
  937. 'MESSAGE' 'Min JC' ('MINIMUM' (RV.INCO.'JC'))   ;
  938. 'MESSAGE' 'Max JD' ('MAXIMUM' (RV.INCO.'JD'))   ;
  939. 'MESSAGE' 'Min JD' ('MINIMUM' (RV.INCO.'JD'))   ;
  940. 'MESSAGE' 'Max AlpD' ('MAXIMUM' (RV.INCO.'AD')) ;
  941. 'MESSAGE' 'Min AlpD' ('MINIMUM' (RV.INCO.'AD')) ;
  942.  
  943.  
  944. SI (EGA GRAPH O);
  945.   'SI' ('EGA' (ntemp MODULO (NITMA/PERIO)) 0 ) ;
  946.  
  947.    ung = vect rv.INCO.'JD' 25. ux uy jaune;
  948.    ungu = 'EXCO' rv.inco.'JD' UY ;
  949.    ungv = 'EXCO' rv.inco.'JC' UY ;
  950.    ungd = 'EXCO' rv.inco.'UD' UY ;
  951.  
  952.    ealpini = 'EVOL' (BLEU) 'CHPO' alpini alpdco;
  953.    Ejd    = 'EVOL' (JAUNE) 'CHPO' ungu alpdco ;
  954.    Ejc    = 'EVOL' (JAUNE) 'CHPO' ungv alpdco ;
  955.  
  956.  
  957.  
  958.     toto = RV.INCO.'AD';
  959.     Ealpd = 'EVOL' 'CHPO' toto alpdco ;
  960.     'DESSIN' (Ealpd  'ET' ealpini) 'TITR'
  961.     'Frac vol phase Dispersee - comparaison solu init' ;
  962.  
  963.    uref = 1.D+2;
  964.    ung = vect RV.INCO.'JC' uref ux uy jaune ;
  965.    trac ung &mt ('CONTOUR' &mt)  titre 'flux volumique JC' ;
  966.  
  967.    'DESSIN' Ejc 'TITRE' 'Profil vertical de Jc';
  968.  
  969.    uref = 1.D+2;
  970.    ung = vect RV.INCO.'JD' uref ux uy jaune ;
  971.    trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JD' ;
  972.  
  973.    'DESSIN' Ejd 'TITRE' 'Profil vertical de Jd ';
  974.  
  975.    pnn = elno $mt RV.INCO.PRES CENTREP1;
  976.    Epp = 'EVOL'(TURQ) 'CHPO' pnn alpdco ;
  977.    dessin Epp 'TITRE' 'pression réduite';
  978.  
  979.  
  980.    'FINSI' ;
  981. FINSI;
  982.  
  983. 'FIN' BLOCKT ;
  984.  
  985.  
  986. ************* calcul de l'erreur ************************************
  987.  
  988. * calcul de la solution analytique.
  989.  
  990. * La vitesse volumique max de sédimentation pour JD
  991.  
  992. jdsolm = ALPD0 '*' (1.D0 '-' ALPD0);
  993. coefsol = 2. * ray0 * ray0 * 9.81 * (rod - roc) / (9. * nu * roc);
  994. jdsolm = coefsol * jdsolm;
  995.  
  996. * Le profile de fraction volumique au temps final NITMA * DELTATPS
  997.  
  998. deplabas = (3. * hh/10.) + (jdsolm * nitma * DELTATPS);
  999. deplahau = hh - (hh/10.) - ((deplabas - (3.* hh/10.))/ alpd0);
  1000.  
  1001. cory = 'COORDONNEE' 2 &mt;
  1002. alpd1 = 'TANH' (20000.D0 * (cory '-' (deplabas))) ;
  1003. alpd1 = (0.999 '-' ALPD0) * (1.D0 '-' alpd1) '/' 2.D0;
  1004. alpd2 = 'TANH' (200.D0 * (cory '-' (deplahau))) ;
  1005. alpd2 = ALPD0 * (1.D0 '-' alpd2) '/' 2.D0;
  1006. alpd = alpd1 '+' alpd2;
  1007. alpfin = alpd;
  1008.  
  1009.  
  1010. * la vitesse de sédimenation pour JD (flux volumique), pour
  1011. * alpd numérique donné
  1012.  
  1013. jdsol = kops RV.INCO.'AD' '*' (kops 1.D0 '-' RV.INCO.'AD');
  1014. coefsol = -2. * ray0 * ray0 * 9.81 * (rod - roc) / (9. * nu * roc);
  1015. jdsol = coefsol * jdsol;
  1016.  
  1017. * la vitesse volumique verticale calculée numériquement
  1018.  
  1019. jdsoln = exco RV.INCO.'JD' uy;
  1020. jdsoln = kcht $mt SCAL SOMMET jdsoln;
  1021.  
  1022. * idem mais avec la fraction volumique analytique
  1023.  
  1024. alpfin = kcht $mt SCAL SOMMET alpfin;
  1025. jdsol1 = kops alpfin '*' (kops 1.D0 '-' alpfin);
  1026. jdsol1 = coefsol * jdsol1;
  1027.  
  1028. * calcul de l'erreur L2 sur la fraction volumique
  1029.  
  1030. DEBPROC CALCURR vit*'CHPOINT' vitsol*'CHPOINT' mod*'MMODEL'
  1031.   mt*'MAILLAGE' coupe*'MAILLAGE';
  1032.   vit = kcht mod SCAL SOMMET vit;
  1033.   vitsol = kcht mod SCAL SOMMET vitsol;
  1034.   udiff = kops vit - vitsol;
  1035.   udif2 = kops udiff * udiff;
  1036.   evudif2 = 'EVOL' 'CHPO' udif2 coupe;
  1037.   ud2 = INTG evudif2;
  1038. *  ud2 = MAXI ud2;
  1039.   usol2 = vitsol * vitsol;
  1040.   eusol2 = 'EVOL' 'CHPO' usol2 coupe;
  1041.   us2 = INTG eusol2 + 1.e-14 ;
  1042. *  us2 = (MAXI us2) + 1.e-14;
  1043.   err = ud2 / us2;
  1044.   err = err**0.5;
  1045. FINPROC err ;
  1046.  
  1047. * erreur en pression (choque inclus) en norme L2.
  1048.   errorp = CALCURR RV.INCO.'AD' alpfin $mt mt alpdco;
  1049.  
  1050. * erreur en flux vol de vitesse JD (choque inclus) 
  1051. * par rapport à la solution semi analytique (formule
  1052. * analytique faisant intervenir AD calculé numériquement).
  1053.   errorv = CALCURR jdsoln jdsol $mt mt alpdco;
  1054.  
  1055. * erreur pour JD par rapport à la solution analytique  
  1056.   errorv1 =  CALCURR jdsoln jdsol1 $mt mt alpdco;
  1057.  
  1058. SI (EGA GRAPH O);
  1059.   trac mt;
  1060. Ealpf = 'EVOL' (BLEU) 'CHPO' alpfin alpdco ;
  1061. 'DESSIN' (Ealpf 'ET' Ealpd) 'TITR'
  1062.  'Frac vol phase Dispersee - comparaison solu analytique' ;
  1063. * on compare maintenant les flux avec le flux analytique mais
  1064. * utilisant la fraction volumique calculée numériquement.
  1065. * ceci permet d'évaluer la résolution intrinsèque de la quantité
  1066. * de mouvement. On pourrait aussi comparer avec la solution
  1067. * analytique.
  1068. Ejdsol = 'EVOL' (BLEU) 'CHPO' jdsol alpdco;
  1069. 'DESSIN' (Ejd et Ejdsol) 'TITR' 
  1070.  'Flux vol phase Dispersee - comparaison solu semi-analytique' ;
  1071. Ejdsol1 = 'EVOL' (BLEU) 'CHPO' jdsol1 alpdco;
  1072. 'DESSIN' (Ejd et Ejdsol1) 'TITR'
  1073.  'Flux vol phase Dispersee - comparaison solu analytique' ;
  1074. FINSI;
  1075.  
  1076. *
  1077. SI (NON COMPLET);
  1078.    SI (errorp > 0.07D0);
  1079.      ERREUR 5;
  1080.    FINSI;
  1081.    SI (errorv > 1.D-3);
  1082.      ERREUR 5;
  1083.    FINSI;
  1084.    SI (errorv1 > 0.16D0);
  1085.      ERREUR 5;
  1086.    FINSI;
  1087. SINON;
  1088.    SI (errorp > 0.03D0);
  1089.      ERREUR 5;
  1090.    FINSI;
  1091.    SI (errorv > 0.001D0);
  1092.      ERREUR 5;
  1093.    FINSI;
  1094.    SI (errorv1 > 0.010);
  1095.      ERREUR 5;
  1096.    FINSI;
  1097. FINSI;
  1098.  
  1099.  
  1100. *OPTI SAUV '/u2/gbm/castem/diphiplan.res' ;
  1101. *REPIX RV ;
  1102. *REPIX RA ;
  1103. *SAUV ;
  1104. fin;
  1105.  
  1106.  
  1107.  
  1108.  
  1109.  
  1110.  
  1111.  
  1112.  
  1113.  
  1114.  
  1115.  
  1116.  
  1117.  
  1118.  
  1119.  
  1120.  
  1121.  
  1122.  
  1123.  
  1124.  

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