Télécharger difasyk2Dax.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : difasyk2Dax.dgibi
  2. opti trac ps ;
  3. opti dime 2 elem qua8 mode axis ;
  4. *'OPTION' 'echo' 0 ;
  5. *
  6. *
  7. **************************************************************************
  8. *
  9. * CAS TEST DIPHASIQUE LIQUIDE-LIQUIDE 2D PLAN
  10. *
  11. ***************************************************************************
  12. *
  13. *
  14. * Présentation: Ce cas test permet la simulation de la sédimention
  15. * de gouttes sphériques de liquide immergé dans
  16. * un milieu continu (autre liquide) de densité moindre.
  17. * Bien entendu les densités peuvent être modifiées.
  18. * les deux liquides sont dans un extracteur
  19. * centrifuges et sont soumis à une
  20. * force centrifuge horizontale
  21. *
  22. *
  23. * Le modèle : Le modèle simulé consiste à résoudre 2 équations de la
  24. * quantité de mouvement sur les flux volumiques des 2
  25. * phases JC et JD. Ensuite on résout une équation de la
  26. * conservation de la masse qui porte sur la fraction
  27. * volumique de la phase dispersée : AD (pour alpha_D).
  28. * La conservation de la fraction voumique de la phase
  29. * continue AC s'en déduite par AD + AC = 1.
  30. * L'équation d'incompressibilité du fluide est
  31. * Div(JC+JD) = 0.
  32. *
  33. * Des équations de mouvements sont rajoutées pour
  34. * les composantes orthoradiales pour la prise en
  35. * compte des forces de Coriolis.
  36. *
  37. * le modèle de coalescence est simplifié, ce qui
  38. * revient à jouer sur la viscosité efficace de l'émulsion
  39. * MUEFF qui tend vers la valeur de la phase pure
  40. * reconstituée lorsque la fraction volumique dépasse
  41. * un certain seuil.
  42. * Le front de sédimentation est rès raide
  43. * notre problème revient regarder évoluer le front
  44. * de sédimentation de gouttes en sédimentation dans
  45. * un extracteur perpétuellement alimenté en mélange
  46. * et pourvu de deux sorties expulsant des phases pures
  47. * séparées.
  48. *
  49. *
  50. * Auteur : Gilles Bernard-Michel
  51. *
  52. * Révision : version 1
  53. *
  54. * date : 21/12/99
  55. *
  56. * REMARQUES :
  57. *
  58. * Nous somme en diphasique et en mode axis
  59. *
  60. * différence entre mueff et hinder pas à jour pour exposant.
  61. *
  62. * AD : ALPD au Sommet
  63. * AC : ALPC au Sommet
  64. *
  65. * implicitation en alpha pour les frottements.
  66. *
  67. * attention meme mueff pour les dex phases ce qui est faux.
  68. * on doit tendre vers les bons mu homogenes.
  69. *
  70. *
  71. ************************ Options générales *********************************
  72. *
  73. *
  74. ** Sortie graphique.
  75. *
  76.  
  77. GRAPH = 'O';
  78.  
  79. *
  80. *test court uniquement.
  81. *
  82.  
  83. COMPLET = FAUX;
  84.  
  85. *
  86. *période des sorties graphiques : tous les NITMA/PERIO pas de temps
  87. *
  88.  
  89. PERIO = 1 ;
  90.  
  91. *
  92. * nombre d'itérations en temps
  93. *
  94.  
  95. NITMA = 1;
  96.  
  97. *
  98. * nombre de points fixes max.
  99. *
  100.  
  101. NITER = 2;
  102.  
  103. *
  104. * Précision du point fixe
  105. *
  106.  
  107. PREC = 1.D-4;
  108.  
  109. *
  110. * Pas de temps
  111. *
  112.  
  113. DELTATPS = 0.2 ;
  114.  
  115. *
  116. ********************* On traite les aspects géométriques ********************
  117. *
  118. ** viscosité dynamique mu/rho (utilisé uniquement pour le maillage).
  119. *
  120.  
  121. nu =1.E-6 ;
  122.  
  123. *
  124. ** La vitesse de rotation initiale. (rad/s)
  125. *
  126.  
  127. omeg = 100 ;
  128.  
  129. *
  130. ** Les paramètres de la géométrie.
  131. *
  132.  
  133. rint = 0.1 ;
  134. rext = 0.2;
  135. hh = 0.044 ;
  136. inj1 = 0.145 ;
  137. inj2 = 0.155;
  138. exi1 = 0.102;
  139. exi2 = 0.198;
  140.  
  141. *
  142. ** Coefficients affectés aux couches d'Ekman et Stewartson.
  143. *
  144.  
  145. coefekma = 8. ;
  146. coefste1 = 4. ;
  147. coefste2 = 2. ;
  148.  
  149. *
  150. ** Densité du maillage.
  151. *
  152.  
  153. npts = 12. ;
  154. d1 = 0.1 / npts;
  155. d2x1 = d1 / 10. ;
  156. d2x2 = d1 / 1.4 ;
  157. d2y = d1 / 2. ;
  158. d2y2 = d2y / 10. ;
  159.  
  160. *
  161. ** Calcul des épaisseurs des couches d'Ekman et de Stewartson.
  162. ** Ces épaisseurs sont des ordres de grandeurs sans dimension.
  163. ** on épaissit un peu les mailles en prenant nu = 10 *nu
  164.  
  165. nu = 10.*nu;
  166. ekma = nu / (omeg * (rext**2));
  167. ddek = ekma**0.5;
  168. stew = ekma**0.25;
  169.  
  170. *
  171. ** Modifications de ces épaisseurs en tenant compte des coefficents
  172. ** introduits.
  173. *
  174.  
  175. dbek = coefekma * ddek * hh ;
  176. fiek = hh * (1. - (coefekma * ddek)) ;
  177. dste = rint * coefste1 * stew ;
  178. fste = rext * coefste2 * stew ;
  179.  
  180. *
  181. ** Les sommets du bol
  182. *
  183.  
  184. p1 = rint 0. ;
  185. p2 = rext 0. ;
  186. p3 = rext hh ;
  187. p4 = rint hh ;
  188.  
  189. *
  190. ** Les milieux des cotés.
  191. *
  192.  
  193. mil1 = ((rint + rext)/2.0);
  194. mil2 = hh/2.0 ;
  195. p6 = rext mil2 ;
  196. p8 = rint mil2 ;
  197. p68 = inj1 mil2 ;
  198. *
  199. ** Les points délimitant les couches d'Ekman.
  200. *
  201.  
  202. p11 = rint dbek ;
  203. p22 = rext dbek ;
  204. p33 = rext fiek ;
  205. p44 = rint fiek ;
  206.  
  207. *
  208. ** Les points d'entrée du fluide et leurs symétriques.
  209. ** Sont rajoutées les couches de Stewartson qui les
  210. ** entourent.
  211. *
  212.  
  213. ent1 = inj1 0. ;
  214. ent2 = inj2 0. ;
  215. sym1 = inj1 hh ;
  216. sym2 = inj2 hh ;
  217. en11 = (inj1 - fste) 0. ;
  218. en22 = (inj2 + fste) 0. ;
  219. sy11 = (inj1 - fste) hh ;
  220. sy22 = (inj2 + fste) hh ;
  221.  
  222. *
  223. ents = inj1 dbek ;
  224. syms = inj1 fiek ;
  225. *
  226.  
  227. *
  228. ** Les points de sortie du fluide et leurs symétriques.
  229. ** Sont rajoutées les couches de Stewartson qui les
  230. ** entourent.
  231. *
  232.  
  233. out1 = exi1 hh ;
  234. out2 = exi2 hh ;
  235. syo1 = exi1 0. ;
  236. syo2 = exi2 0. ;
  237. ou11 = (exi1 '+' fste) hh;
  238. ou22 = (exi2 '-' fste) hh;
  239. so11 = (exi1 '+' fste) 0.D0;
  240. so22 = (exi2 '-' fste) 0.D0;
  241.  
  242.  
  243. *
  244. ** Assemblage des bords du domaine.
  245. *
  246.  
  247. bas1 = p1 'DROIT' dini d2x1 dfin d2x1 syo1
  248. 'DROIT' dini d2x1 dfin d1 so11
  249. 'DROIT' dini d1 dfin d1 en11
  250. 'DROIT' dini d1 dfin d2x2 ent1;
  251. entr = ent1 'DROIT' dini d2x2 dfin d2x2 ent2 ;
  252. bas2 = ent2 'DROIT' dini d2x2 dfin d1 en22
  253. 'DROIT' dini d1 dfin d1 so22
  254. 'DROIT' dini d1 dfin d2x1 syo2
  255. 'DROIT' dini d2x1 dfin d2x1 p2;
  256. cdro = p2 'DROIT' dini d2y dfin d2y p22
  257. 'DROIT' dini d2y dfin d1 p6
  258. 'DROIT' dini d1 dfin d2y p33
  259. 'DROIT' dini d2y dfin d2y p3;
  260. sor2 = p3 'DROIT' dini d2x1 dfin d2x1 out2;
  261. hau = out2 'DROIT' dini d2x1 dfin d1 ou22
  262. 'DROIT' dini d1 dfin d1 sy22
  263. 'DROIT' dini d1 dfin d2x2 sym2
  264. 'DROIT' dini d2x2 dfin d2x2 sym1
  265. 'DROIT' dini d2x2 dfin d1 sy11
  266. 'DROIT' dini d1 dfin d1 ou11
  267. 'DROIT' dini d1 dfin d2x1 out1;
  268. sor1 = out1 'DROIT' dini d2x1 dfin d2x1 p4;
  269. cgau = p4 'DROIT' dini d2y dfin d2y p44
  270. 'DROIT' dini d2y dfin d1 p8
  271. 'DROIT' dini d1 dfin d2y p11
  272. 'DROIT' dini d2y dfin d2y p1;
  273. cot1 = bas1 'ET' entr 'ET' bas2;
  274. cot2 = cdro;
  275. cot3 = sor2 'ET' hau 'ET' sor1;
  276. cot4 = cgau;
  277. cnt = cot1 'ET' cot2 'ET' cot3 'ET' cot4 ;
  278. *
  279. *
  280. *
  281. *borbg = p1 'DROIT' dini d2x1 dfin d2x1 syo1
  282. * 'DROIT' dini d2x1 dfin d1 ent1;
  283. *borbd = ent1 'DROIT' dini d1 dfin d2x1 syo2
  284. * 'DROIT' dini d2x1 dfin d2x1 p2;
  285. *
  286. *bord = p2 'DROIT' dini d2y2 dfin d2y2 p22
  287. * 'DROIT' dini d1 dfin d1 p6
  288. * 'DROIT' dini d1 dfin d1 p33
  289. * 'DROIT' dini d2y2 dfin d2y2 p3;
  290.  
  291. *borhd = p3 'DROIT' dini d2x1 dfin d2x1 out2
  292. * 'DROIT' dini d2x1 dfin d1 sym1;
  293. *borhg = sym1 'DROIT' dini d1 dfin d2x1 out1
  294. * 'DROIT' dini d2x1 dfin d2x1 p4;
  295. *
  296. *borg = p4 'DROIT' dini d2y2 dfin d2y2 p44
  297. * 'DROIT' dini d1 dfin d1 p8
  298. * 'DROIT' dini d1 dfin d1 p11
  299. * 'DROIT' dini d2y2 dfin d2y2 p1;
  300.  
  301. *milieud = ent1 'DROIT' dini d2y2 dfin d2y2 ents
  302. * 'DROIT' dini d1 dfin d1 p68
  303. * 'DROIT' dini d1 dfin d1 syms
  304. * 'DROIT' dini d2y2 dfin d2y2 sym1;
  305. *milieug = 'INVERSE' milieud ;
  306. *
  307. ** On maille avec des rectangles.
  308. *
  309. *mt1 = borbg milieud borhg borg 'DALLER' ;
  310. *mt2 = borbd bord borhd milieug 'DALLER' ;
  311. *mt = mt1 'ET' mt2 ;
  312. mt= cot1 cot2 cot3 cot4 'DALLER' ;
  313. 'TASSER' mt ;
  314. *'TITRE' 'Maillage de la centrifugeuse';
  315. *'TRACER' mt ;
  316.  
  317. *
  318. ** On réoriente les éléments.
  319. *
  320.  
  321. mt = 'ORIENTER' mt ;
  322.  
  323.  
  324. * définition du modèle.
  325.  
  326. &mt = CHANGER mt QUAF ;
  327. *&bord = CHANGER bord QUAF ;
  328. *&borg = CHANGER borg QUAF ;
  329. *&borhd = CHANGER borhd QUAF ;
  330. *&borhg = CHANGER borhg QUAF ;
  331. *&borbd = CHANGER borbd QUAF ;
  332. *&borbg = CHANGER borbg QUAF ;
  333. *&milieud = CHANGER milieud QUAF ;
  334.  
  335. &sor1 = CHANGER sor1 QUAF ;
  336. &sor2 = CHANGER sor2 QUAF ;
  337. &bas1 = CHANGER bas1 QUAF ;
  338. &bas2 = CHANGER bas2 QUAF ;
  339. &cdro = 'CHANGER' cdro QUAF ;
  340. &cgau = 'CHANGER' cgau QUAF ;
  341. &entr = 'CHANGER' entr QUAF ;
  342. &hau = CHANGER hau QUAF ;
  343. &cot1 = CHANGER cot1 QUAF;
  344. &cot2 = CHANGER cot2 QUAF;
  345. &cot3 = CHANGER cot3 QUAF;
  346. &cot4 = CHANGER cot4 QUAF;
  347. &cnt = CHANGER cnt QUAF ;
  348.  
  349. 'ELIMINATION' 1.e-6 (&mt 'ET' &bas1 'ET' &bas2 'ET' &cdro 'ET' &cgau
  350. 'ET' &entr 'ET' &hau 'ET' &cot1 'ET' &cot2
  351. 'ET' &cot3 'ET' &cot4 'ET' &cnt 'ET' &sor1 'ET' &sor2);
  352.  
  353. *'ELIMINATION' 1.e-6 (&mt
  354. *'ET' &borbd 'ET' &borbg 'ET' &bord 'ET' &borhd 'ET' &borhg 'ET' &borg
  355. *'ET' &milieud ) ;
  356.  
  357. $mt = MODE &mt 'NAVIER_STOKES' QUAF ;
  358. $sor1 = MODE &sor1 'NAVIER_STOKES' QUAF ;
  359. $sor2 = MODE &sor2 'NAVIER_STOKES' QUAF ;
  360.  
  361. *
  362. ** On crée le domaine associé au maillage.
  363. *
  364.  
  365. cnt = 'CONTOUR' mt;
  366.  
  367. mt = doma $mt maillage ;
  368.  
  369. lmt = 'CHANGE' mt LIGNE ;
  370. coup = lmt 'POINT' 'DROITE' p8 p6 1.D-5 ;
  371. alpdco = lmt 'ELEM' 'APPUYE' 'STRICTEMENT' coup ;
  372.  
  373.  
  374.  
  375. *
  376. ** on crée une frontière moins un noeud pour le cas
  377. ** de cavités carrés. Il ne faut pas retirer un coin
  378. ** car on peut avoir un champ sortant.
  379. *
  380.  
  381. *$cnt = 'MODELISER' &cnt 'NAVIER_STOKES' QUAF;
  382. *cntcc = 'DOMA' $cnt maillage ;
  383. *cntcc = 'CHANGER' cntcc POI1;
  384.  
  385. *retir = p8 'MANUEL' POI1;
  386. *cntcc = retir 'DIFF' cntcc;
  387.  
  388. sor2c = 'DOMA' $sor2 maillage ;
  389. sor2c = 'CHANGER' sor2c POI1;
  390.  
  391. retir = p3 'MANUEL' POI1;
  392. sor2c = retir 'DIFF' sor2c;
  393.  
  394. *
  395.  
  396.  
  397. *
  398. * *********************************
  399. * PROCEDURE CALCULANT LE MODLO
  400. * *********************************
  401. *
  402. DEBPROC MODULO i*ENTIER j*ENTIER ;
  403. k*ENTIER = i / j ;
  404. mod = i - ( k * j ) ;
  405. FINPROC mod ;
  406.  
  407. *
  408. * ***************************************
  409. * PROCEDURE DE CALCUL DES COEFFICIENTS
  410. * DES EQUATIONS DE QUANTITE DE MOUVEMENT
  411. * ***************************************
  412. *
  413. DEBPROC ACTU RV*'TABLE' $mt*'MMODEL' ;
  414.  
  415. * constantes du probleme
  416. roc = RV.INCO.'ROC' ;
  417. rod = RV.INCO.'ROD' ;
  418. muc = RV.INCO.'MUC' ;
  419. mud = RV.INCO.'MUD' ;
  420. expo = RV.INCO.'EXPOS' ;
  421. farcc = RV.INCO.'FARCC' ;
  422. farcd = RV.INCO.'FARCD' ;
  423.  
  424. * LECTURE DES INCONNUES *
  425.  
  426. JC = RV.INCO.'JC' ;
  427. JD = RV.INCO.'JD' ;
  428.  
  429. * implicitation en alpha pour les frottements.
  430. ALPM = RV.INCO.'AM' ;
  431. *ALPD = 'COPIER' RV.INCO.'AD';
  432. ALPD = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
  433. RV.INCO.'AC' = kops 1.D0 '-' RV.INCO.'AD';
  434.  
  435. MASK1D = ALPD MASQUE 'INFERIEUR' (1.D0 '-' 1.D-15) ;
  436. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  437. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  438. '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
  439. MASK1D = ALPD MASQUE 'SUPERIEUR' 1D-15 ;
  440. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  441. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  442. '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
  443.  
  444. ALPC = 'KOPS' 1.D0 '-' ALPD;
  445.  
  446. *** COUPLAGE DYNAMIQUE ***
  447.  
  448. * VISCOSITE DYNAMIQUE EFFICACE
  449. MUEFF = MUC * RV.INCO.'MUCOEF' ;
  450. RV.INCO.'MUEFF' = 'KCHT' $mt SCAL 'SOMMET' MUEFF ;
  451. RV.INCO.'SNUFFC' = 'KOPS' MUEFF '/' ROC ;
  452. RV.INCO.'SNUFFC' = 'KOPS' ALPC * RV.INCO.'SNUFFC';
  453. RV.INCO.'NUEFFC' = 'KCHT' $mt scal centre
  454. ('NOEL' $mt RV.INCO.'SNUFFC');
  455. RV.INCO.'SNUFFD' = 'KOPS' MUEFF '/' ROD ;
  456. RV.INCO.'SNUFFD' = 'KOPS' ALPD * RV.INCO.'SNUFFD';
  457. RV.INCO.'NUEFFD' = 'KCHT' $mt scal centre
  458. ('NOEL' $mt RV.INCO.'SNUFFD');
  459.  
  460.  
  461. ** UC ET UD : VITESSES INSTANTANEES **
  462.  
  463. SEUIL = RV.INCO.'AD' 'MASQUE' 'SUPERIEUR' 0.05D0;
  464. SEUIL = 'KCHT' $mt SCAL SOMMET SEUIL ;
  465.  
  466. ALPCX = NOMC 'UX' ALPC ;
  467. ALPCY = NOMC 'UY' ALPC ;
  468. ALPCV = ALPCX 'ET' ALPCY ;
  469. UC = 'KOPS' RV.INCO.'JC' '/' ALPCV ;
  470. UC = 'KOPS' SEUIL '*' UC ;
  471. RV.INCO.'UC' = 'KCHT' $mt VECT 'SOMMET' UC ;
  472.  
  473. ALPDX = NOMC 'UX' ALPD ;
  474. ALPDY = NOMC 'UY' ALPD ;
  475. ALPDV = ALPDX 'ET' ALPDY ;
  476. UD = 'KOPS' RV.INCO.'JD' '/' ALPDV ;
  477. UD = 'KOPS' SEUIL '*' UD ;
  478. RV.INCO.'UD' = 'KCHT' $mt VECT 'SOMMET' UD ;
  479.  
  480.  
  481. * COEFF DE FROTTEMENT
  482. *on calcule une fraction volumique ALPD prenant
  483. *en compte les parois verticales = max (alpd alpdm).
  484. maskp = ALPD 'MASQUE' 'SUPERIEUR' ALPM ;
  485. maskp = 'KCHT' $mt scal sommet maskp;
  486. maskpi = 'KOPS' 1.D0 '-' maskp ;
  487. ALPDF = 'KOPS' maskp * ALPD ;
  488. ALPFF = 'KOPS' maskpi * ALPM ;
  489. ALPDF = 'KOPS' ALPDF '+' ALPFF ;
  490.  
  491. *HINDER = ALPDF ;
  492. HINDER = 0.5*((KOPS 0.001D0 '+' ALPC)
  493. + ( abs (KOPS 0.001D0 '-' ALPC)));
  494. *HINDER = ALPC;
  495. HINDER = 'KCHT' $mt 'SCAL' 'SOMMET' HINDER ;
  496. HINDER = 'KOPS' HINDER '**' expo ;
  497. HINDER = 'KOPS' MUC '*' HINDER ;
  498. kft = 'KCHT' $mt SCAL 'SOMMET' (9.D0/2.D0) ;
  499. kft = 'KOPS' ('KOPS' RV.INCO.'AI' '**' 2) '*' kft;
  500.  
  501. RV.INCO.'FRCD' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' -1.D0
  502. '*' ('KOPS' kft '/' roc) ) ;
  503. RV.INCO.'FRCD1' = 'KOPS' RV.INCO.'FRCD' '*' HINDER ;
  504. RV.INCO.'FRCD2' = -1.0D0 * ('KOPS' RV.INCO.'FRCD1' '*' ALPD) ;
  505. RV.INCO.'FRCD' = 'KOPS' RV.INCO.'FRCD1' '+' RV.INCO.'FRCD2' ;
  506. RV.INCO.'MFRCD' = -0.0D0 * ('KOPS' RV.INCO.'FRCD1' '*' RV.INCO.'JD') ;
  507.  
  508. RV.INCO.'FRDC' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' -1.D0
  509. '*' ('KOPS' kft '/' rod)) ;
  510. RV.INCO.'FRDC' = 'KOPS' RV.INCO.'FRDC' '*' HINDER ;
  511. RV.INCO.'MFRDC' = 0.0D0 * ('KOPS' RV.INCO.'FRDC' '*' RV.INCO.'JC') ;
  512. RV.INCO.'FRDC' = 1.0D0 * ('KOPS' RV.INCO.'FRDC' '*' ALPD) ;
  513.  
  514. RV.INCO.'FRCC' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' 1.D0
  515. '*' ('KOPS' kft '/' roc) ) ;
  516. RV.INCO.'FRCC' = 'KOPS' RV.INCO.'FRCC' '*' HINDER ;
  517. RV.INCO.'MFRCC' = 0.0D0 * ('KOPS' RV.INCO.'FRCC' '*' RV.INCO.'JC') ;
  518. RV.INCO.'FRCC' = 1.0D0 * ('KOPS' RV.INCO.'FRCC' '*' ALPD) ;
  519.  
  520. RV.INCO.'FRDD' = 'KCHT' $mt SCAL 'SOMMET' ('KOPS' 1.D0
  521. '*' ('KOPS' kft '/' rod)) ;
  522. RV.INCO.'FRDD1' = 'KOPS' RV.INCO.'FRDD' '*' HINDER ;
  523. RV.INCO.'FRDD2' = -1.0D0 * ('KOPS' RV.INCO.'FRDD1' '*' ALPD) ;
  524. RV.INCO.'FRDD' = 'KOPS' RV.INCO.'FRDD1' '+' RV.INCO.'FRDD2' ;
  525. RV.INCO.'MFRDD' = -0.0D0 * ('KOPS' RV.INCO.'FRDD1' '*' RV.INCO.'JD') ;
  526.  
  527. RV.INCO.'KFT' = 'KCHT' $mt SCAL SOMMET kft;
  528.  
  529.  
  530. * termes de pression
  531.  
  532. PUC = 'KOPS' -1.D0 '*' RV.INCO.'AC' ;
  533. PUC = 'KOPS' PUC '/' ROC ;
  534. RV.INCO.'PUC' = kcht $mt SCAL SOMMET PUC ;
  535.  
  536. PUD = 'KOPS' -1.D0 '*' RV.INCO.'AD' ;
  537. PUD = 'KOPS' PUD '/' ROD ;
  538. RV.INCO.'PUD' = kcht $mt SCAL SOMMET PUD ;
  539. RV.INCO.'PUD' = 1.0 * RV.INCO.'PUD';
  540.  
  541. PNN = ELNO $mt RV.INCO.PRES CENTREP1;
  542. PNN = KCHT $mt SCAL SOMMET PNN;
  543. GPNN = KOPS GRAD PNN $MT;
  544. GPNN = KCHT $mt VECT CENTRE GPNN;
  545. GPNN = elno $mt GPNN;
  546. RV.INCO.'GPUD' = kcht $mt VECT SOMMET (0.5D0 * (1/ROC) * GPNN);
  547.  
  548. * terme source pesanteur signe + si TOIMP
  549. * signe - si MDIA ou dans DFDT. on met de l'ordre
  550. * 2, au lieu de alpc alpd(n+1), on a
  551. * (alpc-alpd)alpd(n+1) + alpd*alpd
  552.  
  553. farchc = 'KOPS' ('KOPS' RV.INCO.'AC' '-' RV.INCO.'AD') '*' farcc;
  554. farchci = -1.D0 * farchc ;
  555. *partie expl
  556. farchce = 'KOPS' ('KOPS' RV.INCO.'AD' * RV.INCO.'AD') '*' farcc;
  557. farchce = -1.D0 * farchce;
  558. RV.INCO.'FARCHCI' = 'KCHT' $mt VECT 'SOMMET' farchci ;
  559. RV.INCO.'FARCHCE' = 'KCHT' $mt VECT 'SOMMET' farchce ;
  560.  
  561. farchd = 'KOPS' ('KOPS' RV.INCO.'AC' '-' RV.INCO.'AD') '*' farcd;
  562. farchdi = -1.D0 * farchd ;
  563. *partie expl
  564. farchde = 'KOPS' ('KOPS' RV.INCO.'AD' * RV.INCO.'AD') '*' farcd;
  565. farchde = -1.D0 * farchde;
  566. RV.INCO.'FARCHDI' = 'KCHT' $mt VECT 'SOMMET' farchdi ;
  567. RV.INCO.'FARCHDE' = 'KCHT' $mt VECT 'SOMMET' farchde ;
  568.  
  569. * terme correction laplacien sur ortho radiale
  570.  
  571. RV.INCO.'TCORIC' = 'KOPS' RV.INCO.'SNUFFC' '*' RV.INCO.'INVCAR';
  572. RV.INCO.'TCORID' = 'KOPS' RV.INCO.'SNUFFD' '*' RV.INCO.'INVCAR';
  573.  
  574. FINPROC RV;
  575.  
  576. *************************
  577. *** DONNEES PHYSIQUES ***
  578. *************************
  579. * vitesse angulaire en rad/s
  580. omeg = 100.0D0 ;
  581. * rotation, coef de la force de Coriolis à gauche de l'égalité pour r,z
  582. domeg = ((-2.D0*omeg) 0.D0) ;
  583. * rotation, coef de la force de Coriolis à gauche de l'égalité pour theta
  584. domegt = ((2.D0*omeg) 0.D0);
  585. * champs des 1/r2
  586. invrcar = 'KCHT' $mt SCAL SOMMET ('COORDONNEE' 1 &mt);
  587. invrcar = 'KOPS' invrcar * invrcar;
  588. uni = 'KCHT' $mt SCAL SOMMET 1.D0;
  589. invrcar = 'KOPS' uni '/' invrcar ;
  590. * accélération centrifuge
  591. accelr = 'KCHT' $mt SCAL SOMMET ('COORDONNEE' 1 &mt);
  592. omeg2 = omeg * omeg;
  593. omeg2 = (omeg2 0.D0);
  594. accelr = 'KOPS' omeg2 '*' accelr;
  595. * formule Ishii trainée
  596. expo = -4.1D0;
  597. * formule viscosité
  598. expomu = -0.0D0;
  599. *coeff de viscosité (pour lisser avec mauvais maillage)
  600. MUCOEF = 10.D0 ;
  601. *Densité volumique
  602. ROC = 1000.D0 ;
  603. ROD = 1100.D0 ;
  604. epsilc = (ROC - ROD) / ROC ;
  605. epsild = (ROD - ROC) / ROD ;
  606. * Viscosite cinématique et dynamique de la phase liquide
  607. MUC = 0.001D0 ;
  608. NUC = MUC / ROC ;
  609. MUD = 0.001D0 ;
  610. NUD = MUD / ROD ;
  611. * Rayon initial des gouttes
  612. ray0 = 1.D-5 ;
  613. * Concentration initiale de la phase dispersée
  614. ALPD0 = 0.3D0 ;
  615. * Concentration initiale de la phase continue
  616. ALPC0 = 1 - ALPD0 ;
  617. *C0 = 1.D-2 ;
  618. C0 = 1.D-2;
  619. * Coeff. de Pénalisation *
  620. EPSS=1.D-9;
  621. * Valeur minimale de EPSN *
  622. MIEPS = 1.e-5 ;
  623. * Valeur maximale de EPSN *
  624. MAEPS = 0.99999 ;
  625. * Vicosité maximale et truc maximal *
  626. ALPMAX = 0.9D0 ;
  627. MUMAX = 1.D-1 ;
  628. RAY2 = 1.D-3;
  629. * fraction vol max
  630. alphadm = 0.9900000001D0 ;
  631. alphacm = MIEPS ;
  632. * profil de fraction volumique
  633. corx = 'COORDONNEE' 1 &mt;
  634. alpd1 = 'TANH' (20000.D0 * (corx '-' (7.95*rext/10.D0))) ;
  635. alpd1 = (0.998 '-' ALPD0) * (1.D0 '+' alpd1) '/' 2.D0;
  636. alpd2 = 'TANH' (20000.D0 * (corx '-' (7.05D0 * rext/10.D0))) ;
  637. alpd2 = ALPD0 * (1.D0 '+' alpd2) '/' 2.D0;
  638. alpd = 0.001D0 '+' alpd1 '+' alpd2;
  639. alpd = 'KCHT' $mt scal sommet alpd;
  640. alpini = alpd;
  641. * concentration aire interf
  642. ALPI0 = 3.D0 * alpd / ray0 ;
  643. *'TRACER' ALPD &mt;
  644. * poussée d'Archimède multipliée par +1 car TOIMP à droite de l'égaZZlité
  645. gravos = 'KCHT' $mt VECT SOMMET (0.D0 -0.D0) ;
  646. archi = 'KOPS' gravos '+' accelr ;
  647. farcc = epsilc* archi ;
  648. farcd = epsild* archi ;
  649.  
  650. * PAS DE TEMPS *
  651.  
  652. DT = DELTATPS ;
  653.  
  654. ************************************
  655. *** SYSTEMES D'EQUATIONS ET C.L. ***
  656. ************************************
  657. *
  658. *
  659. *************************************
  660. *** SYSTEME PORTANT SUR : ***
  661. *** UN, VN, P ***
  662. *************************************
  663.  
  664. * EQUATION DE CONSERVATION DE LA QUANTIE DE MVT PHASE CONT
  665. *
  666. * composantes r,z
  667. RV = 'EQEX' $mt ITMA 1
  668. OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'CENTREE'
  669. ZONE $mt OPER 'DFDT' 1.D0 'JCP' 'JCPP' 'DT' 'UC' 1.D-15 INCO 'JC'
  670. ZONE $mt OPER 'MDIA' 'FARCHCI' INCO 'AD' 'JC'
  671. ;
  672. RV = 'EQEX' RV
  673. OPTI 'EFM1' 'IMPL' 'CENTREE'
  674. ZONE $mt OPER 'MDIA' 'FRCD' INCO 'JD' 'JC'
  675. ZONE $mt OPER 'MDIA' 'FRCC' INCO 'JC' 'JC'
  676. ZONE $mt OPER 'MDIA' domeg INCO 'CT' 'JC'
  677. OPTI 'EF' 'IMPL' 'CENTREE'
  678. ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JC'
  679. ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JD'
  680. *ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JC'
  681. ;
  682.  
  683. * composante theta, otho radiale
  684. RV = 'EQEX' RV
  685. OPTI 'EFM1' 'IMPL' 'BDF2' 'CENTREE'
  686. ZONE $mt OPER 'DFDT' 1.D0 'CTP' 'CTPP' 'DT' 'CT' 1.D-15 INCO 'CT'
  687. ZONE $mt OPER 'MDIA' 'FRCD' INCO 'JDT' 'CT'
  688. ZONE $mt OPER 'MDIA' 'FRCC' INCO 'CT' 'CT'
  689. 'OPTI' 'EFM1' 'IMPL' 'CENTREE'
  690. ZONE $mt OPER 'MDIA' 'domegt' INCO 'JC' 'CT'
  691. ZONE $mt OPER 'MDIA' 'TCORIC' INCO 'CT' 'CT'
  692. ZONE $mt OPER 'MDIA' 'TCORIC' INCO 'JDT' 'CT'
  693. ;
  694.  
  695. RV = EQEX RV
  696. OPTI 'EF' 'IMPL' 'CENTREE'
  697. ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'CT'
  698. ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'JDT'
  699. *ZONE $mt OPER 'LAPN' 'NUEFFC' INCO 'CT'
  700. OPTI 'EF' 'IMPL' 'CENTREP1' 'CENTREE'
  701. ZONE $mt OPER 'KMBT' 'PUC' INCO 'PRES' 'JC'
  702. ;
  703.  
  704. *
  705. * EQUATION DE CONSERVATION DE LA QUANTITE DE MVT PHASE DISP
  706. *
  707. * composantes r, z
  708. RV = EQEX RV
  709. OPTI 'EFM1' 'IMPL' 'BDF2' 'CONS' 'CENTREE'
  710. ZONE $mt OPER 'DFDT' 1.D0 'JDP' 'JDPP' 'DT' 'UD' 1.D-15 INCO 'JD'
  711. ZONE $mt OPER 'MDIA' 'FARCHDI' INCO 'AD' 'JD'
  712. ZONE $mt OPER 'MDIA' 'FRDC' INCO 'JC' 'JD'
  713. ZONE $mt OPER 'MDIA' 'FRDD' INCO 'JD' 'JD'
  714. ZONE $mt OPER 'MDIA' domeg INCO 'JDT' 'JD'
  715. ;
  716.  
  717. * ATTENTION LES LAPLACIENS PAS BONS, LAPN JC + JD -> JC
  718. * ET RIEN POUR JD (A VERIF)
  719. * composante ortho-radiale theta
  720. RV = EQEX RV
  721. OPTI 'EFM1' 'IMPL' 'BDF2' 'CENTREE'
  722. ZONE $mt OPER 'DFDT' 1.D0 'DTP' 'DTPP' 'DT' 'CT' 1.D-15 INCO 'JDT'
  723. ZONE $mt OPER 'MDIA' 'FRDC' INCO 'CT' 'JDT'
  724. ZONE $mt OPER 'MDIA' 'FRDD' INCO 'JDT' 'JDT'
  725. ZONE $mt OPER 'MDIA' 'domegt' INCO 'JD' 'JDT'
  726. ZONE $mt OPER 'MDIA' 'TCORID' INCO 'CT' 'JDT'
  727. ZONE $mt OPER 'MDIA' 'TCORID' INCO 'JDT' 'JDT'
  728. ;
  729.  
  730.  
  731. RV = EQEX RV
  732. OPTI 'EF' 'IMPL' 'CENTREE'
  733. ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'CT'
  734. ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JDT'
  735. *ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JDT'
  736. OPTI 'EF' 'IMPL' 'CENTREE'
  737. ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JD'
  738. ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JC'
  739. *ZONE $mt OPER 'LAPN' 'NUEFFD' INCO 'JD'
  740. ;
  741.  
  742.  
  743. RV = EQEX RV
  744. OPTI 'EF' 'IMPL' 'CENTREP1' 'CENTREE'
  745. ZONE $mt OPER 'KMBT' 'PUD' INCO 'PRES' 'JD'
  746. ;
  747.  
  748. *
  749. * EQUATION DE CONSERVATION DU VOLUME (INCOMPRESSIBLE)
  750. *
  751.  
  752. RV = EQEX RV
  753. OPTI 'EFM1' 'IMPL' 'CENTREP1' 'CENTREE'
  754. ZONE $mt OPER 'KMAB' 1.D0 0. INCO 'JD' 'PRES'
  755. ZONE $mt OPER 'KMAB' 1.D0 0. INCO 'JC' 'PRES'
  756. ;
  757.  
  758.  
  759. *
  760. * EQUATION DE CONSERVATION DE LA FRACTION VOLUMIQUE
  761. *
  762.  
  763. RV = EQEX RV
  764. OPTI 'EFM1' 'IMPL' 'CONS' 'BDF2' 'SUPG' 'CMD' 0.2
  765. ZONE $mt OPER 'DFDT' 1.D0 'ADP' 'ADPP' 'DT' 'UD' 1.D-15 INCO 'AD'
  766. OPTI 'EF' 'IMPL' 'CENTREP1' 'CONS' 'SUPG' 'CMD' 0.2
  767. ZONE $mt OPER 'KMAB' 1.D0 0.D0 'SOMM' 1.D-15 INCO 'JD' 'AD'
  768. ;
  769.  
  770.  
  771. *** FIN SECOND SYSTEME ***
  772. **************************
  773.  
  774. *
  775. * CONDITION AUX LIMITES
  776. *
  777.  
  778. RV = EQEX RV
  779. 'CLIM'
  780. 'JC' 'VIMP' (&bas1 'ET' &bas2 'ET' &hau) 0.
  781. 'JC' 'UIMP' (&bas1 ET &bas2 ET &hau ET &cot2 ET &cot4) 0.D0
  782. 'JD' 'VIMP' (&bas1 'ET' &bas2 'ET' &hau) 0.
  783. 'JD' 'UIMP' (&bas1 ET &bas2 ET &hau ET &cot2 ET &cot4) 0.D0
  784. 'JD' 'VIMP' (&entr) 1.65D-4
  785. 'JC' 'VIMP' (&entr) 1.65D-4
  786. ;
  787. RV = EQEX RV
  788. 'CLIM'
  789. 'CT' 'TIMP' &cnt 0.D0
  790. 'JDT' 'TIMP' &cnt 0.D0
  791. ;
  792.  
  793. RV = EQEX RV
  794. 'CLIM'
  795. 'AD' 'TIMP' (&cot2) 0.99
  796. 'AD' 'TIMP' (&cot4) 0.
  797. 'AD' 'TIMP' &sor1 0.D0
  798. 'AD' 'TIMP' sor2c 0.99D0
  799. 'AD' 'TIMP' &entr ALPD0
  800. ;
  801.  
  802.  
  803. *
  804. *******************************************
  805. * INITIALISATION DES TABLES DES INCONNUES
  806. *******************************************
  807. *
  808.  
  809. *
  810. RV.INCO = TABLE INCO ;
  811.  
  812. * constantes du probleme
  813. RV.INCO.'ROC' = ROC ;
  814. RV.INCO.'ROD' = ROD ;
  815. RV.INCO.'MUC' = MUC ;
  816. RV.INCO.'MUD' = MUD ;
  817. RV.INCO.'EXPOS' = expo ;
  818. RV.INCO.'FARCC' = farcc ;
  819. RV.INCO.'FARCD' = farcd ;
  820. RV.INCO.'ALPHADMAX' = alphadm ;
  821.  
  822. *
  823. RV.'DT' = DELTATPS ;
  824. RV.INCO.'DT' = RV.'DT' ;
  825. RV.INCO.'GDT' = -1.D0 * epsild * 9.81D0 * RV.'DT' ;
  826.  
  827. RV.INCO.'JC' = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
  828. RV.INCO.'JD' = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
  829. RV.INCO.'JCP' = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC' ;
  830. RV.INCO.'JCPP' = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JC' ;
  831. RV.INCO.'JDP' = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD' ;
  832. RV.INCO.'JDPP' = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JD' ;
  833. RV.INCO.'CT' = 'KCHT' $mt SCAL 'SOMMET' 1.e-25 ;
  834. RV.INCO.'CTP' = 'KCHT' $mt SCAL 'SOMMET' 1.e-25 ;
  835. RV.INCO.'CTPP' = 'KCHT' $mt SCAL 'SOMMET' 1.e-25 ;
  836. RV.INCO.'JDT' = 'KCHT' $mt SCAL 'SOMMET' 1.e-25 ;
  837. RV.INCO.'DTP' = 'KCHT' $mt SCAL 'SOMMET' 1.e-25 ;
  838. RV.INCO.'DTPP' = 'KCHT' $mt SCAL 'SOMMET' 1.e-25 ;
  839.  
  840. RV.INCO.'PRES' = 'KCHT' $mt SCAL 'CENTREP1' 1.e-5 ;
  841.  
  842.  
  843. RV.INCO.'AC' = 'KCHT' $mt SCAL 'SOMMET' (1.D0 '-' ALPD) ;
  844. RV.INCO.'AD' = 'KCHT' $mt SCAL 'SOMMET' ALPD ;
  845. RV.INCO.'RALPD' = 'KCHT' $mt SCAL 'SOMMET' ALPD ;
  846. RV.INCO.'RALPC' = 'KCHT' $mt SCAL 'SOMMET' (1.D0 '-' ALPD) ;
  847. RV.INCO.'ADP'= 'KCHT' $mt SCAL 'SOMMET' RV.INCO.'AD' ;
  848. RV.INCO.'ADPP'= 'KCHT' $mt SCAL 'SOMMET' RV.INCO.'AD' ;
  849.  
  850. RV.INCO.'RAY' = 'KCHT' $mt SCAL 'SOMMET' ray0 ;
  851.  
  852. * AI est l'inverse du rayon ;
  853. RV.INCO.'AI' = 'KCHT' $mt SCAL 'SOMMET' (1.D0 '/' ray0) ;
  854.  
  855. RV.INCO.'FARCHCI' = 'KCHT' $mt vect SOMMET (0.D0 0.D0) ;
  856. RV.INCO.'FARCHDI' = 'KCHT' $mt vect SOMMET (0.D0 0.D0) ;
  857. RV.INCO.'FARCHCE' = 'KCHT' $mt vect SOMMET (0.D0 0.D0) ;
  858. RV.INCO.'FARCHDE' = 'KCHT' $mt vect SOMMET (0.D0 0.D0) ;
  859.  
  860. * INITIALISATION DES COEFFICIENTS
  861.  
  862. RV.INCO.'DR2' = kcht $mt VECT 'CENTRE' ((d1*d1) (d1*d1));
  863. RV.INCO.'MUCOEF' = kcht $mt SCAL 'SOMMET' MUCOEF;
  864. *
  865. RV.INCO.'CNUL' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  866. RV.INCO.'JNUL' = 'KCHT' $mt VECT 'SOMMET' (1.D-25 1.D-25) ;
  867. RV.INCO.'JUNI' = 'KCHT' $mt VECT 'SOMMET' (1.D0 1.D0) ;
  868.  
  869. RV.INCO.'CNUL1' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  870. RV.INCO.'CNUL2' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  871. RV.INCO.'CNUL3' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  872. RV.INCO.'CNUL4' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  873. RV.INCO.'SOUC' = 'KCHT' $mt VECT 'CENTRE' (1.e-25 1.e-25) ;
  874. RV.INCO.'SOUD' = 'KCHT' $mt VECT 'CENTRE' (1.e-25 1.e-25) ;
  875.  
  876. RV.INCO.'PUC' = 'KCHT' $mt SCAL 'SOMMET' 1.e-5 ;
  877. RV.INCO.'PUD' = 'KCHT' $mt SCAL 'SOMMET' 1.e-5 ;
  878. RV.INCO.'GPUD' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  879.  
  880. RV.INCO.'UC' = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
  881. RV.INCO.'UD' = 'KCHT' $mt VECT 'SOMMET' (1.e-25 1.e-25) ;
  882.  
  883. RV.INCO.'NUEFFC' = 'KCHT' $mt SCAL 'CENTRE' 0.0D0 ;
  884. RV.INCO.'NUEFFD' = 'KCHT' $mt SCAL 'CENTRE' 0.0D0 ;
  885. RV.INCO.'SNUFFC' = 'KCHT' $mt SCAL 'SOMMET' 0.0D0 ;
  886. RV.INCO.'SNUFFD' = 'KCHT' $mt SCAL 'SOMMET' 0.0D0 ;
  887.  
  888. RV.INCO.'MUDIFC' = 'KCHT' $mt SCAL 'CENTRE' 0.D0 ;
  889. RV.INCO.'MUDIFD' = 'KCHT' $mt SCAL 'CENTRE' 0.D0 ;
  890.  
  891. RV.INCO.'FRCC' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  892. RV.INCO.'FRCD1' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  893. RV.INCO.'FRCD2' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  894. RV.INCO.'FRCD' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  895. RV.INCO.'FRDC' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  896. RV.INCO.'FRDD1' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  897. RV.INCO.'FRDD2' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  898. RV.INCO.'FRDD' = 'KCHT' $mt SCAL 'SOMMET' 0.D0 ;
  899. RV.INCO.'MFRCC' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  900. RV.INCO.'MFRCD' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  901. RV.INCO.'MFRDC' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  902. RV.INCO.'MFRDD' = 'KCHT' $mt VECT 'SOMMET' (0.D0 0.D0) ;
  903.  
  904. RV.INCO.'TCORIC' = 'KCHT' $mt VECT SOMMET (0.D0 0.D0) ;
  905. RV.INCO.'TCORID' = 'KCHT' $mt VECT SOMMET (0.D0 0.D0) ;
  906.  
  907. RV.INCO.'INVCAR' = 'KCHT' $mt SCAL SOMMET invrcar ;
  908. RV.INCO.'domegt' = domegt;
  909.  
  910. *COEFFICIENT DE RELAXATION
  911. ALFI = 1.0D0 ;
  912. ALFID = 1.0D0 ;
  913. ALFII = 1.D0 ;
  914.  
  915. * Compteur des pas de temps
  916. ntemp = 0 ;
  917.  
  918. *** Calcul de ALPM frac vol pour la paroi ***********************
  919. *** distance à la paroi la plus proche **************************
  920.  
  921. corx = 'COOR' &mt 1;
  922. corx = 'KCHT' $mt scal sommet corx;
  923. corxi = 'KOPS' rext '-' corx;
  924. corx = 'KOPS' corx '-' rint;
  925. masko = corx 'MASQUE' 'INFERIEUR' corxi ;
  926. masko = 'KCHT' $mt scal sommet masko;
  927. maskoi = 'KOPS' 1.D0 '-' masko;
  928. corx = 'KOPS' masko * corx ;
  929. corxi = 'KOPS' maskoi * corxi ;
  930. corx = 'KOPS' corx '+' corxi ;
  931.  
  932. ratio = 'KOPS' corx '/' RV.INCO.'RAY' ;
  933. ratio = 'KOPS' ratio '+' 3.D0 ;
  934. drel = 2.D0 * (alphadm**0.34);
  935. ratio = 'KOPS' ratio '/' drel;
  936. * on fabrique une fraction volumique telle que les écarts
  937. * entre particules soient identiques à la distance à la paroi
  938. * la plus proche.
  939. RV.INCO.'AM' = 'KOPS' ratio ** (-3);
  940.  
  941. *'TRACER' RV.INCO.'AM' &mt ('CONTOUR' &mt);
  942.  
  943.  
  944. *
  945. **** BOUCLES DES ITERATIONS INTERNES ET EN TEMPS ****************
  946. *
  947. REPETER BLOCKT NITMA ;
  948. *
  949. ntemp = ntemp + 1 ;
  950. 'MESSAGE' 'Itération externe N°' ntemp;
  951.  
  952. ntempi = 0 ;
  953.  
  954. * On stocke les variables au pas de temps ntemp
  955. VJCP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JCP' ;
  956. VJCPP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JCPP' ;
  957. VJDP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JDP' ;
  958. VJDPP = 'KCHT' $mt 'VECT' 'SOMMET' RV.INCO.'JDPP' ;
  959. VJCTP = 'COPIER' RV.INCO.'CTP' ;
  960. VJCTPP = 'COPIER' RV.INCO.'CTPP' ;
  961. VJDTP = 'COPIER' RV.INCO.'DTP' ;
  962. VJDTPP = 'COPIER' RV.INCO.'DTPP' ;
  963. VALPDP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADP' ;
  964. VALPDPP = 'KCHT' $mt SCAL SOMMET RV.INCO.'ADPP';
  965.  
  966. * On effectue les itérations de point fixe
  967. *
  968. 'REPETER' BLOCKI NITER ;
  969. *
  970. ntempi = ntempi + 1 ;
  971. *
  972. *
  973. *************************************
  974. *** SYSTEME PORTANT SUR : ***
  975. *** UN, VN, P ***
  976. *************************************
  977. *
  978. * On charge les résultats du vrai pas de temps précédent
  979. RV.INCO.'JCP' = 'KCHT' $mt 'VECT' SOMMET VJCP;
  980. RV.INCO.'JCPP' = 'KCHT' $mt 'VECT' SOMMET VJCPP;
  981. RV.INCO.'CTP' = 'KCHT' $mt SCAL SOMMET VJCTP;
  982. RV.INCO.'CTPP' = 'KCHT' $mt SCAL SOMMET VJCTPP;
  983. RV.INCO.'JDP' = 'KCHT' $mt 'VECT' SOMMET VJDP;
  984. RV.INCO.'JDPP' = 'KCHT' $mt 'VECT' SOMMET VJDPP;
  985. RV.INCO.'DTP' = 'KCHT' $mt SCAL SOMMET VJDTP;
  986. RV.INCO.'DTPP' = 'KCHT' $mt SCAL SOMMET VJDTPP;
  987. RV.INCO.'ADP' = 'KCHT' $mt SCAL SOMMET VALPDP ;
  988. RV.INCO.'ADPP' = 'KCHT' $mt SCAL SOMMET VALPDPP ;
  989. * resultats de l'iter interne préc.
  990. JCIP = 'KCHT' $mt VECT SOMMET RV.INCO.'JC' ;
  991. JDIP = 'KCHT' $mt 'VECT' SOMMET RV.INCO.'JD' ;
  992. ALPDIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AD' ;
  993. ALPCIP = 'KCHT' $mt SCAL SOMMET RV.INCO.'AC' ;
  994.  
  995. *RAYC = 'KOPS' ALPCIP '/' ('KOPS' 1.D-15 '+' ALPIIP) ;
  996. *RAYC = 'KOPS' ('KOPS' ray0 '+' RAYC) '+' ('ABS' ('KOPS' ray0 '-' RAYC));
  997. *RAYC = 'KCHT' $mt SCAL 'SOMMET' RAYC ;
  998. *RAYC = 'KOPS' RAYC '*' 3 ;
  999. *RV.INCO.'RAYC' = RAYC ;
  1000.  
  1001.  
  1002. *
  1003. * On réactualise les données.
  1004. *
  1005.  
  1006. RV = ACTU RV $mt;
  1007.  
  1008. * On affecte le terme source expl de QDM par le biais de DFDT
  1009. * coef 2 car ordre 2 en temps dans dfdt. JCPP sinon la poussee
  1010. * archi est rechargéee * fois lors des reinitialisations
  1011. RV.INCO.'JCPP' = 'KOPS' RV.INCO.'JCP' '+'
  1012. ('KOPS' (2.D0*RV.INCO.'DT') * RV.INCO.'FARCHCE') ;
  1013.  
  1014. RV.INCO.'JDPP' = 'KOPS' RV.INCO.'JDP' '+'
  1015. ('KOPS' (2.D0*RV.INCO.'DT') * RV.INCO.'FARCHDE') ;
  1016.  
  1017.  
  1018. EXEC RV ;
  1019. 'MESSAGE' 'Boucle interne N°' ntempi;
  1020.  
  1021. ALPD = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
  1022. MASK1D = ALPD MASQUE 'INFERIEUR' (1.D0 '-' 1.D-15) ;
  1023. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  1024. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  1025. '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
  1026. MASK1D = ALPD MASQUE 'SUPERIEUR' 1D-15 ;
  1027. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  1028. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  1029. '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
  1030.  
  1031. * On masque AD
  1032. RV.INCO.'AD' = ALPD;
  1033.  
  1034.  
  1035.  
  1036. * CALCUL DE L'ERREUR RELATIVE MAX. *
  1037. ERR = 0.D0 ;
  1038. ERRC = 'KOPS' ('KOPS' RV.INCO.'JC' '-' JCIP)
  1039. '/' ('MAXIMUM' ('ABS' JCIP) '+' 1.D-10) ;
  1040. ERRD = 'KOPS' ('KOPS' RV.INCO.'JD' '-' JDIP)
  1041. '/' ('MAXIMUM' ('ABS' JDIP) '+' 1.D-10) ;
  1042. ERRC = 'MAXIMUM' ('ABS' ERRC) ;
  1043. ERRD = 'MAXIMUM' ('ABS' ERRD) ;
  1044. 'SI' (ERRC > ERR) ;
  1045. ERR = ERRC ;
  1046. VAR = 'MOTS' JC ;
  1047. 'FINSI' ;
  1048.  
  1049. 'SI' (ERRD > ERR) ;
  1050. ERR = ERRD ;
  1051. VAR = 'MOTS' JD ;
  1052. 'FINSI' ;
  1053. ERRA = 'KOPS' ('KOPS' RV.INCO.'AD' '-' ALPDIP)
  1054. '/' ('MAXIMUM' ('ABS' ALPDIP) '+' 1.D-10) ;
  1055. ERRA = 'MAXIMUM' ('ABS' ERRA) ;
  1056. 'SI' (ERRA > ERR) ;
  1057. ERR = ERRA ;
  1058. VAR = 'MOTS' ALPHAD ;
  1059. 'FINSI' ;
  1060.  
  1061. * CALCUL DE L'ERREUR RELATIVE MAX. *
  1062.  
  1063.  
  1064. *'LIST' VAR;
  1065. 'MESSAGE' ERR;
  1066. 'MESSAGE' ERRC ERRD ERRA;
  1067. 'SI' (ERR <EG PREC) ;
  1068. QUITTER BLOCKI ;
  1069. 'FINSI' ;
  1070.  
  1071. *
  1072. 'FIN' BLOCKI ;
  1073. *
  1074.  
  1075. 'MESSAGE' 'NBR ITER INTERNES' ntempi;
  1076. 'MESSAGE' 'ERREUR RELATIVE MAX' ERR ;
  1077. 'MESSAGE' 'Max JC' ('MAXIMUM' (RV.INCO.'JC')) ;
  1078. 'MESSAGE' 'Min JC' ('MINIMUM' (RV.INCO.'JC')) ;
  1079. 'MESSAGE' 'Max JD' ('MAXIMUM' (RV.INCO.'JD')) ;
  1080. 'MESSAGE' 'Min JD' ('MINIMUM' (RV.INCO.'JD')) ;
  1081. 'MESSAGE' 'Max AlpD' ('MAXIMUM' (RV.INCO.'AD')) ;
  1082. 'MESSAGE' 'Min AlpD' ('MINIMUM' (RV.INCO.'AD')) ;
  1083.  
  1084.  
  1085. SI (EGA GRAPH O);
  1086. 'SI' ('EGA' (ntemp MODULO (NITMA/PERIO)) 0 ) ;
  1087.  
  1088. ung = vect rv.INCO.'JD' 25. ux uy jaune;
  1089. ungu = 'EXCO' rv.inco.'JD' UX ;
  1090. ungv = 'EXCO' rv.inco.'JC' UX ;
  1091. ungd = 'EXCO' rv.inco.'UD' UX ;
  1092.  
  1093. Ealpd = 'EVOL' (BLEU) 'CHPO' rv.inco.'AD' alpdco ;
  1094. ealpini = 'EVOL' (BLEU) 'CHPO' alpini alpdco;
  1095. Ejd = 'EVOL' (JAUNE) 'CHPO' ungu alpdco ;
  1096. Ejc = 'EVOL' (JAUNE) 'CHPO' ungv alpdco ;
  1097. Eeff = 'EVOL' (VERT) 'CHPO' rv.inco.'MUEFF' alpdco ;
  1098. Etransp = 'EVOL'(TURQ) 'CHPO' ungd alpdco ;
  1099.  
  1100. 'DESSIN' Ejd 'TITRE' 'Jd';
  1101. 'DESSIN' Ejc 'TITRE' 'Jc';
  1102. 'DESSIN' Etransp 'TITRE' 'Champ transportant Ud' ;
  1103.  
  1104. toto = RV.INCO.'AD';
  1105. Ealpd = 'EVOL' 'CHPO' toto alpdco ;
  1106. 'DESSIN' (Ealpd 'ET' ealpini) 'TITR'
  1107. 'Fraction volumique de la phase Dispersee' ;
  1108.  
  1109. uref = 2.0D0;
  1110. ung = vect RV.INCO.'JD' uref ux uy jaune ;
  1111. trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JD' 'NCLK' ;
  1112.  
  1113. uref = 2.0D0;
  1114. ung = vect RV.INCO.'JC' uref ux uy jaune ;
  1115. trac ung &mt ('CONTOUR' &mt) titre 'flux volumique JC' 'NCLK' ;
  1116.  
  1117. pnn = elno $mt RV.INCO.PRES CENTREP1;
  1118. Epp = 'EVOL'(TURQ) 'CHPO' pnn alpdco ;
  1119. 'DESSIN' Epp 'TITRE' 'pression';
  1120.  
  1121. ALPD = 'KCHT' $MT SCAL SOMMET RV.INCO.'AD';
  1122. MASK1D = 'MASQUE' ALPD 'INFERIEUR' (1.D0 '-' 1.D-15) ;
  1123. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  1124. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  1125. '+' ('KOPS' (1.D0 '-' 1.D-15) '*' ('KOPS' 1.D0 '-' MASK1D));
  1126. MASK1D = 'MASQUE' ALPD 'SUPERIEUR' 1D-15 ;
  1127. MASK1D = 'KCHT' $mt SCAL 'SOMMET' MASK1D ;
  1128. ALPD = 'KOPS' ('KOPS' MASK1D '*' ALPD)
  1129. '+' ('KOPS' 1.D-15 * ('KOPS' 1.D0 '-' MASK1D));
  1130. MONTAGNE ALPD mt 0.12 ('PROG' -2. -2. 1.) 'TITRE' 'tot' 'SUPER' ;
  1131.  
  1132. 'FINSI' ;
  1133.  
  1134. 'FINSI';
  1135.  
  1136. *
  1137. 'FIN' BLOCKT ;
  1138.  
  1139. fin;
  1140.  
  1141.  
  1142.  
  1143.  
  1144.  
  1145.  
  1146.  
  1147.  
  1148.  
  1149.  
  1150.  
  1151.  
  1152.  
  1153.  

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