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;
  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.  

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