Télécharger transport6.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : transport6.dgibi
  2. ******************** CAS TEST : transport6.dgibi ************************
  3. *
  4. GRAPH = faux ;
  5. CRIT1 = 1.D-3 ;
  6. 'SAUTER' 'PAGE';
  7.  
  8. ************************************************************************
  9. *
  10. * Petit programme testant le transport en milieu poreux avec
  11. * - Retard non linéaire suivant une isotherme de type Langmuir.
  12. *
  13. * pour un cas test 2D simple :
  14. * Une petite source baignant dans un champ de vitesse uniforme.
  15. *
  16. * P4 *---------------* P3
  17. * | |
  18. * | |
  19. * | |
  20. * P5 |---------------| P2
  21. * | |
  22. * | |
  23. * | |
  24. * | |
  25. * P6 *----* P9 - - - * P11
  26. * | | |
  27. * P7 *----* P8 - - - * P10
  28. * | |
  29. * | |
  30. * P0 *---------------* P1
  31. *
  32. ************************************************************************
  33. 'OPTION' 'ECHO' 0 ;
  34. 'SAUTER' 'PAGE';
  35. 'TITRE' 'Validation du transport';
  36. * ==============
  37. * | PARAMETRES |
  38. * ==============
  39. *
  40. * Paramètres internes :
  41. * =====================
  42. 'OPTION' 'ECHO' 0 'DIME' 2 'ELEM' 'QUA4' ;
  43. 'OPTION' ISOV 'SURFACE' ;
  44. 'OPTION' 'TRAC' PSC ;
  45.  
  46. * Paramètres de convivialité :
  47. * ============================
  48. * Trace les limites de solubilité ?
  49. TrLimSol= Vrai ;
  50. * Trace les profils de concentration suivant certaines lignes ?
  51. TrLEvol = Vrai ;
  52. * Trace les évolutions de concentration en certaines points ?
  53. TrPEvol = Vrai ;
  54. * Trace les champs de concentration en soluté, précipité ?
  55. TrConc = vrai ;
  56. TrPrec = vrai ;
  57. * Trace la quantité de précipité dissout ?
  58. TrLib = vrai ;
  59. * Trace le flux à travers une ligne ?
  60. TrFlux = Vrai ;
  61.  
  62. SI ('NON' GRAPH) ;
  63. TrLEvol = faux ;
  64. TrPEvol = faux ;
  65. TrConc = faux ;
  66. TrPrec = faux ;
  67. TrLib = faux ;
  68. TrFlux = faux ;
  69. TrLimSol= faux ;
  70. FINSI;
  71.  
  72. * Grandeurs physiques :
  73. * =====================
  74. * VITESSE DE DARCY
  75. UXImp = 0.;
  76. UYImp = 1.;
  77.  
  78. * POROSITE
  79. ValPor = 1.;
  80.  
  81. *-- Facteur de forme :
  82. FacForme = 1.;
  83.  
  84. * RETARD
  85. *-- Coefficient de retard (éventuellement scalaire) :
  86. * Isotherme de Langmuir : F(C) = [(R-1) C] / [1 + (R-1)C/Fsat]
  87. * F étant exprimé en mol/dm3 de milieu.
  88. * où R est la pente à l'origine
  89. * et Fsat est la limite d'adsorption
  90. * Isotherme de Freundlich : F(C) = (R-1) C^(1/n)
  91. * F étant exprimé en mol/dm3 de milieu.
  92. * Attention, R n'est pas un coefficient de retard car il n'y a
  93. * aucune zone de linéarité, et ce n'est pas non plus un nombre
  94. * pur. Sa dimension est (mol/l)^(1 - 1/n)
  95. * L'abus de langage est caractérisé dans ce cas.
  96. * On pose alpharet = 1/n
  97.  
  98. * Cas Langmuir :
  99. CofRetar = 2. ;
  100. Fsat = .5 ;
  101. * Cas Freundlich :
  102. KFreund = 2. ;
  103. AFreund = .5 ;
  104.  
  105. * LIMITE DE SOLUBILITE
  106. *-- Limite de solubilité en mol/l (global et fond):
  107. * Cette valeur ne sert que si la dissolution est instantanée ou
  108. * d'ordre 1.
  109. CSat1 = 1. ;
  110. CSat2 = 0.01 ;
  111. Distrib = 'expo';
  112.  
  113. * DIFFUSIVITE
  114. *-- Diffusivité moléculaire du traceur en m²/s
  115. DifMolec = 2.d0;
  116.  
  117. * NATURE DISSOLUTION
  118. *-- Nature de la dissolution (si existe) et logiques associés :
  119. NaturDis = 'Pas de dissolution';
  120. *NaturDis = 'Imposée';
  121. *NaturDis = 'Ordre 1';
  122. *NaturDis = 'Instant_Pénal';
  123. *NaturDis = 'Instant_PFDis';
  124.  
  125. Solub = 'NEG' NaturDis 'Pas de dissolution';
  126. SoluA = 'EGA' NaturDis 'Imposée';
  127. SoluP = 'EGA' NaturDis 'Ordre 1';
  128. SoluI = ('EGA' NaturDis 'Instant_PFDis')
  129. OU ('EGA' NaturDis 'Instant_Pénal');
  130. SoluL = SoluP OU SoluI;
  131.  
  132. * PARAM. NUM. DISSOLUTION
  133. *-- Grandeur(s) associée(s) à la dissolution :
  134. Methode = 'Point fixe Dissolution';
  135. * Critère d'arrêt (relatif) du point fixe
  136. EpsSolub = 4.D-2;
  137. * Nombre maxi d'itérations
  138. NbItMax = 100;
  139.  
  140. * QUANTITES INITIALES
  141. *-- Concentration en marqueur dans la source à t=0 (mol/l d'eau):
  142. CActif0 = CSat1;
  143.  
  144. *-- Précipité de marqueur dans la source à t=0 (mol/dm² de milieu):
  145. SI Solub;
  146. PActif0 = 10. * CSat1 ;
  147. FINSI;
  148.  
  149. * TEMPS DE CALCUL
  150. * initial, final, et pas de temps :
  151. TIni = 0.;
  152. TFin = 10.;
  153. DelTaT = 1. ;
  154.  
  155. * METHODES NUMERIQUES
  156. *-- Coef. d'implicitation (1 = implicite, 0 = explicite)
  157. * Diffusion :
  158. TetaDif = 1. ;
  159. * Convection :
  160. TetaCnv = 1. ;
  161.  
  162. *-- Coefficient d'élimination :
  163. EpsElim = .1 ;
  164.  
  165. *-- Précision de la machine :
  166. EpsPrec = 1.D-14;
  167.  
  168. * ============
  169. * | MAILLAGE |
  170. * ============
  171. * Dimensions (m) :
  172. * ----------------
  173. LX1 = 1.;
  174. LXM = 5.;
  175. LY1 = 4.;
  176. LY2 = 5.;
  177. LY3 = 8.;
  178. LYM = 12.;
  179.  
  180. * Découpages :
  181. * ------------
  182. N = 1;
  183. NX = 'ENTIER' (N * LXM);
  184. NY = 'ENTIER' (N * LYM + 1);
  185. NX1 = 'MAXIMUM' ('LECT' ('ENTIER' (NX * LX1 / LXM)) 1);
  186. NX2 = NX - NX1;
  187. NY1 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * LY1 / LYM)) 1);
  188. NY2 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LY2 - LY1) / LYM)) 1);
  189. NY3 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LY3 - LY2) / LYM)) 1);
  190. NY4 = 'MAXIMUM' ('LECT' ('ENTIER' (NY * (LYM - LY3) / LYM)) 1);
  191. NY = NY1 + NY2 + NY3 + NY4;
  192.  
  193. * Points :
  194. * --------
  195. p0 = 0. 0.;
  196. p1 = LXM 0.;
  197. p2 = LXM LY3;
  198. p3 = LXM LYM;
  199. p4 = 0. LYM;
  200. p5 = 0. LY3;
  201. p6 = 0. LY2;
  202. p7 = 0. LY1;
  203. p8 = LX1 LY1;
  204. p9 = LX1 LY2;
  205. p10 = LXM LY1;
  206. p11 = LXM LY2;
  207.  
  208. * Droites :
  209. * ---------
  210. 'OPTION' 'COULEUR' 'BLAN';
  211. D4D5 = 'DROITE' P4 NY4 P5;
  212. D5D6 = 'DROITE' P5 NY3 P6;
  213. D6D7 = 'DROITE' P6 NY2 P7;
  214. D7D0 = 'DROITE' P7 NY1 P0;
  215. D2D5 = 'DROITE' P2 NX P5;
  216. 'OPTION' 'COULEUR' 'TURQ';
  217. D0D1 = 'DROITE' P0 NX P1;
  218. D1D10 = 'DROITE' P1 NY1 P10;
  219. D10D11= 'DROITE' P10 NY2 P11;
  220. D11D2 = 'DROITE' P11 NY3 P2;
  221. D2D3 = 'DROITE' P2 NY4 P3;
  222. D3D4 = 'DROITE' P3 NX P4;
  223. D7D8 = 'DROITE' P7 NX1 P8;
  224. D8D9 = 'DROITE' P8 NY2 P9;
  225. D9D6 = 'DROITE' P9 NX1 P6;
  226. D10D8 = 'DROITE' P10 NX2 P8;
  227. D11D9 = 'DROITE' P11 NX2 P9;
  228. D2D5 = 'DROITE' P2 NX P5;
  229. D8D7 = 'INVERSE' D7D8;
  230. D8D10 = 'INVERSE' D10D8;
  231. D9D8 = 'INVERSE' D8D9;
  232. D6D9 = 'INVERSE' D9D6;
  233. D9D11 = 'INVERSE' D11D9;
  234. D5D2 = 'INVERSE' D2D5;
  235. D4D3 = 'INVERSE' D3D4;
  236.  
  237. *-- Pour les conditions aux limites :
  238. LInf = D0D1;
  239. Lsup = D3D4;
  240. LDroit = D1D10 ET D10D11 ET D11D2 ET D2D3;
  241. LGau = 'INVERSE' (D4D5 ET D5D6 ET D6D7 ET D7D0);
  242.  
  243. *-- Pour faire la somme des flux :
  244. * pour le tranfert d'activité (loin des conditions aux limites) :
  245. LFlux = D5D2;
  246.  
  247. * Surfaces maillées :
  248. * -------------------
  249. SActif = 'ORIENTER' ('COULEUR' ('DALLER' D7D8 D8D9 D9D6 D6D7) 'ROUG') ;
  250. SFond = 'ORIENTER' ('COULEUR' ('DALLER' D2D3 D3D4 D4D5 D5D2) 'VERT') ;
  251. SCoeur = 'ORIENTER' ('COULEUR' (
  252. ('DALLER' D0D1 D1D10 (D10D8 ET D8D7) D7D0)
  253. ET ('DALLER' D8D10 D10D11 D11D9 D9D8)
  254. ET ('DALLER' (D6D9 ET D9D11) D11D2 D2D5 D5D6)
  255. ) 'TURQ' ) ;
  256. SGran = 'COULEUR' (SCoeur et SFond) 'TURQ' ;
  257. Mtot = SActif ET SGran;
  258.  
  259. LTrame = 'COULEUR' ((CONT SActif) ET (CONT SGran)) 'BLAN';
  260.  
  261. *
  262. * Elements contenant tous les points de calcul
  263.  
  264. QMtot = 'CHAN' Mtot 'QUAF' ;
  265. QFond = 'CHAN' SFond 'QUAF' ;
  266. QActif = 'CHAN' SActif 'QUAF' ;
  267. QGran = 'CHAN' SGran 'QUAF' ;
  268. QInf = 'CHAN' LInf 'QUAF' ;
  269. QDro = 'CHAN' LDroit 'QUAF' ;
  270. QSup = 'CHAN' Lsup 'QUAF' ;
  271. QGau = 'CHAN' LGau 'QUAF' ;
  272. QLFlux = 'CHAN' LFlux 'QUAF' ;
  273. ELIM 0.001 ( QMtot ET QActif ET QFond ET QGran ET QLFlux ET
  274. QGau ET QInf ET QDro ET QSup) ;
  275.  
  276. * Domaines :
  277. * ----------
  278. * Modèle (Mo) :
  279. MoHydro = 'MODELE' QMtot DARCY ISOTROPE ;
  280. MoFond = 'MODELE' QFond DARCY ISOTROPE ;
  281. MoActif = 'MODELE' QActif DARCY ISOTROPE ;
  282. MoGran = 'MODELE' QGran DARCY ISOTROPE ;
  283. MoInf = 'MODELE' QInf DARCY ISOTROPE ;
  284. MoDro = 'MODELE' QDro DARCY ISOTROPE ;
  285. MoSup = 'MODELE' QSup DARCY ISOTROPE ;
  286. MoGau = 'MODELE' QGau DARCY ISOTROPE ;
  287. MoFlux = 'MODELE' QLFlux DARCY ISOTROPE ;
  288.  
  289. * Longueurs, Normales et Orientations des arêtes du domaine ;
  290. Chyb1 = 'DOMA' MoHydro 'SURFACE';
  291. Chyb2 = 'DOMA' MoHydro 'NORMALE';
  292. VolChPro = 'DOMA' MoHydro 'VOLUME';
  293. CeHydro = 'DOMA' MoHydro 'CENTRE';
  294. CeFlux = 'DOMA' MoFlux 'CENTRE';
  295.  
  296. * Création des normales à la ligne sur laquelle on mesure les flux :
  297. NleFlux = 'MANU' 'CHPO' CeFlux 2 'UX' 0. 'UY' +1. NATURE DIFFUS;
  298. DFluxTIT = 'CHAINE' 'la ligne y=' 'FORMAT' '(F5.0)'
  299. ('COOR' 2 P2) 'm';
  300.  
  301. * Points où l'on suivra la concentration :
  302. * ----------------------------------------
  303. PEvol = 'TABLE';
  304. * Centre du combustible :
  305. PEvol.1 = 'TABLE';
  306. PEvol.1 .'T' = 'CHAINE' 'dans la source' ;
  307. PEvol.1 .'P' = 'CENTRE';
  308. PEvol.1 .'M' = ( (.5 * LX1) ((LY2 + LY1) / 2.) );
  309. PEvol.1 .'M' = CeHydro 'POINT' 'PROC' (PEvol.1 .'M');
  310. * Juste après la transition aval :
  311. PEvol.2 = 'TABLE';
  312. PEvol.2 .'T' = 'CHAINE' 'post-transition aval' ;
  313. PEvol.2 .'P' = 'CENTRE';
  314. PEvol.2 .'M' = CeHydro 'POINT'
  315. 'PROC' (P5 'PLUS' (0.05 0.05));
  316.  
  317. * Lignes où capturer les évolutions des concentrations :
  318. * ------------------------------------------------------
  319. LEvol = 'TABLE';
  320.  
  321. 'REPETER' Bcl8 ('NBEL' LGau);
  322. pt = CeHydro 'POINT' 'PROC'
  323. ((LGau 'POINT' &Bcl8) 'PLUS' (0 0.01));
  324. SI (&Bcl8 > 2);
  325. LGauC = LGauC ET ('MANU' 'SEG2' pt0 pt);
  326. FINSI;
  327. SI (&Bcl8 'EGA' 2);
  328. LGauC = ('MANU' 'SEG2' pt0 pt);
  329. FINSI;
  330. pt0 = pt;
  331. 'FIN' Bcl8;
  332. 'OUBLIER' pt0;
  333. 'OUBLIER' pt;
  334.  
  335. LEvol.1 = 'TABLE';
  336. LEvol.1 .M = 'COULEUR' LGauC 'BLAN';
  337. LEvol.1 .P = 'CENTRE';
  338. LEvol.1 .L = 'Y (m)';
  339. LEvol.1 .T = 'CHAINE' 'longitudinale' ;
  340.  
  341. * ===============
  342. * | HYDRAULIQUE |
  343. * ===============
  344. * La donnée est U. Pour obtenir VDarcy et QF, on simule une résolution
  345. * de Darcy. On prend pour cela K = 1 et grad h = U :
  346.  
  347. * Matériau associé (Ma) :
  348. MaHydro = 'MATERIAU' MoHydro 'K' 1.;
  349.  
  350. * Matrice hybride élémentaire :
  351. MeHydro = 'MHYB' MoHydro MaHydro;
  352.  
  353. * Trace de Concentration :
  354. FAHydro= 'DOMA' MoHydro 'FACE' ;
  355. XX YY = 'COOR' FAHydro;
  356. Un = 'MANU' 'CHPO' FAHydro 1 'SCAL' 1.;
  357. THHydro = (-1.) * (UXImp * XX) - (UYImp * YY) + Un;
  358. THHydro = 'NOMC' THHydro 'TH';
  359.  
  360. * Concentration :
  361. XX YY = 'COOR' CeHydro;
  362. Un = 'MANU' 'CHPO' CeHydro 1 'SCAL' 1.;
  363. HHydro = (-1.) * (UXImp * XX) - (UYImp * YY) + Un;
  364. HHydro = 'NOMC' HHydro 'H';
  365.  
  366. * Flux et vitesses :
  367. QFHydro = 'HDEB' MoHydro MeHydro HHydro THHydro ;
  368. VCDarcy = 'HVIT' MoHydro QFHydro ;
  369.  
  370. * =============
  371. * | TRANSPORT |
  372. * =============
  373. * Modèle propre au transport :
  374. * ----------------------------
  375. MoTransp = 'MODELE' QMtot DARCY ISOTROPE ;
  376.  
  377. * Matrice Porosité cinématique (Type ChamElem) :
  378. * ----------------------------------------------
  379. CeTransp= 'DOMA' MoTransp 'CENTRE' ;
  380. MaPor1 = 'MANU' 'CHPO' CeTransp 1 'CK' ValPor NATURE DIFFUS;
  381. MaPor = 'KCHA' MoTransp MaPor1 'CHAM';
  382.  
  383. * Matrice coef. de retard (Type chpo centre ou scalaire) :
  384. * ---------------------------------------------------------
  385. Maret1 = CofRetar ;
  386.  
  387. * Porosité * volume de maille
  388. * ----------------------------
  389. MaPorVol= VolChPro * ('EXCO' 'CK' MaPor1 'SCAL');
  390.  
  391. * Temps :
  392. * -------
  393. * Temps à calculer
  394. LiTCalc = 'PROG' TIni 'PAS' DeltaT Tfin;
  395. NTCalc = 'DIME' LiTCalc;
  396. 'REMPLACER' LiTCalc NTCalc Tfin;
  397.  
  398. * Ecriture de ces temps sous forme de table ;
  399. TabTCalc = 'TABLE';
  400. 'REPETER' Bcl26 NTCalc;
  401. Id = &Bcl26 - 1;
  402. Tps = 'EXTR' LiTCalc &Bcl26;
  403. TabTCalc.Id = Tps;
  404. 'FIN' Bcl26;
  405.  
  406. * Temps sauvegardés
  407. LiTSauv = LiTCalc ;
  408. NTSauv = NTCalc ;
  409.  
  410. * Temps de tracage des champs
  411. LiIndGph = 'LECT' 0 5 10 ;
  412. NTGrph = 'DIME' LiIndGph;
  413.  
  414. * Temps de tracage des évolutions
  415. LiTEvol = litsauv ;
  416. liindEvl = 'LECT' 0 pas 1 (ntsauv - 1) ;
  417. NTEvol = 'DIME' litevol;
  418.  
  419. * Temps servant à la création des chargements
  420. TSup = 1.2 * TFin;
  421. DelTDef = 0.001 * TSup;
  422.  
  423. * ====================================
  424. * | CONSTRUCTION DES DONNEES DE BASE |
  425. * ====================================
  426. * Dissolution / Libération ;
  427. * --------------------------
  428. *-- Matrice Limite de solubilité ;
  429. * Variation exponentielle décroissante de la source vers l'aval :
  430. XX YY = 'COOR' CeHydro;
  431. Sup2 = 'POINT' YY 'EGINF' LY2;
  432. MaSol2 = 'MANU' 'CHPO' Sup2 1 'H' CSat1 NATURE DIFFUS;
  433. Sup3 = 'POINT' YY 'SUPERIEUR' LY2;
  434. YY3 = 'REDU' (YY - LY2) Sup3;
  435. a = ('LOG' (CSat2 / CSat1)) / (LYM - LY2);
  436. MaSol3 = CSat1 * ('EXP' (a * YY3));
  437. MaSol3 = 'NOMC' MaSol3 'H' NATURE DIFFUS;
  438. MaSol1 = MaSol2 ET MaSol3;
  439.  
  440. * Matériau associé au transport :
  441. * -------------------------------
  442. * Terme diffusif pur (diffusivité effective) = Porosité * Diffusivité :
  443. DifEff = DifMolec * FacForme * MaPor;
  444. MaTransp = 'MATERIAU' MoTransp 'K' DifEff ;
  445.  
  446. * Matrice hybride élémentaire :
  447. * -----------------------------
  448. MeTransp = 'MHYB' MoTransp MaTransp;
  449.  
  450. * Conditions Initiales :
  451. * ----------------------
  452. *-- Traces de concentration (quelconque car le transport est implicité).
  453. TC0 = 'MANU' 'CHPO' FAHydro 1 'TH' 0. 'NATURE' 'DISCRET';
  454.  
  455. *-- Concentration aux centres des éléments (sert à définir la quantité
  456. * initiale d'élément) :
  457. Conc01 = 'MANU' 'CHPO' CEHydro 1 'H' 0. 'NATURE' 'DISCRET';
  458. Conc02 = 'MANU' 'CHPO' (DOMA MoActif 'CENTRE' ) 1 'H' CActif0
  459. 'NATURE' 'DISCRET';
  460. Conc0 = Conc01 ET Conc02;
  461. Conc0 = 'NOMC'
  462. ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Conc0 'SCAL')) 'H';
  463.  
  464. *-- Flux diffusif initial (aux centres des faces) - décoratif puisque la
  465. * diffusion est implicite :
  466. QFDif0 = 'MANU' 'CHPO' FAHydro 1 'FLUX' 0. 'NATURE' 'DISCRET';
  467.  
  468. *-- Précipité :
  469. SI Solub;
  470. Prec01 = 'MANU' 'CHPO' ( 'DOMA' MoActif 'CENTRE') 1 'H' PActif0
  471. 'NATURE' 'DISCRET';
  472. Prec02 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET';
  473. Prec0 = Prec01 ET Prec02;
  474. Prec0 = 'NOMC'
  475. ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Prec0 'SCAL')) 'H';
  476. FINSI;
  477.  
  478. *-- Libération de soluté à partir du précipité (décoratif aussi):
  479. SI Solub;
  480. Libe0 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET';
  481. FINSI;
  482.  
  483. * Conditions aux limites (type CHARGEMENT) :
  484. * ------------------------------------------
  485. * Concentration imposée en amont et en aval, flux nul ailleurs.
  486. Ceinf= DOMA MoInf 'CENTRE' ;
  487. Cesup= DOMA Mosup 'CENTRE' ;
  488. CCL = 'MANU' 'CHPO' (Ceinf et Cesup) 1 'TH' 0.
  489. 'NATURE' 'DISCRET';
  490. RigCL = 'BLOQUE' (Ceinf et Cesup) 'TH';
  491. CharCL = 'DEPIMPOSE' RigCL CCL ;
  492. ChgCL = 'CHARG' CharCL
  493. ('EVOL' 'MANU' ('PROG' 0. TSup) ('PROG' 1. 1.));
  494.  
  495. * ==============
  496. * | Résolution |
  497. * ==============
  498. *-- Table de transport :
  499. Transp = 'TABLE';
  500. Transp.'TEMPS' = 'TABLE';
  501. Transp.'TRACE_CONC' = 'TABLE';
  502. Transp.'CONCENTRATION' = 'TABLE';
  503. Transp.'FLUX' = 'TABLE';
  504. SI Solub;
  505. Transp.'PRECIPITE' = 'TABLE';
  506. Transp.'DISSOLUTION' = 'TABLE';
  507. FINSI;
  508.  
  509. Transp.'SOUSTYPE' = 'DARCY_TRANSITOIRE';
  510. Transp.'MODELE' = MoTransp;
  511. Transp.'CARACTERISTIQUES' = MaTransp;
  512. Transp.'POROSITE' = MaPor;
  513. Transp.'CONVECTION' = QFHydro;
  514. Transp.'COEF_RETARD' = MaRet1;
  515. Transp.'LANGMUIR' = FSat;
  516. SI SoluL;
  517. Transp.'LIMITE_SOLUBILITE' = MaSol1;
  518. FINSI;
  519. SI SoluP;
  520. Transp.'COEF_DISSOLUTION' = MaCoDis1;
  521. FINSI;
  522.  
  523. * Conditions initiales :
  524. Transp.'TEMPS'. 0 = TIni;
  525. Transp.'TRACE_CONC'. 0 = TC0;
  526. Transp.'CONCENTRATION'. 0 = Conc0;
  527. Transp.'FLUX'. 0 = QFDif0;
  528. SI Solub;
  529. Transp.'PRECIPITE'. 0 = Prec0;
  530. Transp.'DISSOLUTION'. 0 = Libe0;
  531. FINSI;
  532.  
  533. * Conditions aux limites :
  534. Transp.'BLOCAGE' = RigCL;
  535. Transp.'TRACE_IMPOSE' = ChgCL;
  536.  
  537. * Paramètres numériques :
  538. Transp.'THETA_DIFF' = TetaDif;
  539. Transp.'THETA_CONVECTION' = TetaCnv;
  540. SI SoluP;
  541. Transp.'THETA_DISS' = TetaDis;
  542. FINSI;
  543. SI ('EGA' Methode 'Point fixe Dissolution');
  544. Transp.'EPSILON' = EpsSolub;
  545. Transp.'ITMAX' = NbItMax;
  546. FINSI;
  547. Transp.'TEMPS_CALCULES' = LiTCalc;
  548. Transp.'TEMPS_SAUVES' = LiTSauv;
  549.  
  550. * ==========
  551. * | CALCUL |
  552. * ==========
  553. DARCYTRA Transp;
  554. * ===================
  555. * | POST-TRAITEMENT |
  556. * ===================
  557. MESS ' ';
  558. MESS 'Post-Traitement :';
  559. MESS '-----------------';
  560. MESS ' ';
  561.  
  562. IndTSauv = 'INDEX' (Transp.temps);
  563.  
  564. * Préparation tracé :
  565. * -------------------
  566. * Calcul du flux à travers une certaine ligne :
  567. * ---------------------------------------------
  568. * Les propriétés de la ligne sont stockées dans DFlux et NleFlux.
  569. * Le nom de la ligne doit etre mis dans DFlux.'TITRE' .
  570. NORHydro = DOMA MoTransp 'NORMALE';
  571. Ori2 = PSCAL ('REDU' NORHydro CeFlux ) NleFlux
  572. ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY');
  573. FlInt = 0.;
  574. liabs2 = 'ENLEVER' litcalc 1;
  575. 'REPETER' Bcl21 (NTSauv - 1);
  576. i = &Bcl21 + 1;
  577. Id = IndTSauv.i;
  578. Tps = Transp.Temps.Id;
  579. Id1 = IndTSauv.(i - 1);
  580. Tp1 = Transp.Temps.(id - 1);
  581. dt = Tps - Tp1;
  582.  
  583. *-- Flux diffusif et convectif à l'instant d'indice Id
  584. * (diffusion implicite TetaDif et convection implicite TetaCnv) :
  585. * méthode par produit scalaire flux * normale sur 1 ligne
  586.  
  587. * Débit diffusif sortant :
  588. DebD1 = ((1. - TetaDif) * Transp.Flux.Id1)
  589. + (TetaDif * Transp.Flux.Id);
  590. DebD2 = 'REDU' DebD1 CeFlux;
  591. DebD3 = 'RESULT' (DebD2 * Ori2);
  592. FlD = ('MAXIMUM' DebD3);
  593. Fl = FlD;
  594.  
  595. * Débit convectif sortant :
  596. THCnv = ((1. - TetaCnv) * Transp.Trace_Conc.Id1)
  597. + (TetaCnv * Transp.Trace_Conc.Id);
  598. DebC1 = ('EXCO' 'TH' THCnv 'SCAL') * Transp.Convection;
  599. DebC2 = 'REDU' DebC1 CeFlux;
  600. DebC3 = 'RESULT' (DebC2 * Ori2);
  601. FlC = ('MAXIMUM' DebC3);
  602. Fl = FlC + FlD;
  603.  
  604. FlInt = FlInt + (Fl * dt);
  605.  
  606. SI (i 'EGA' 2);
  607. LiFlC = 'PROG' FlC;
  608. LiFlD = 'PROG' FlD;
  609. LiFl = 'PROG' Fl;
  610. LiFlI = 'PROG' FlInt;
  611. SINON;
  612. LiFlC = 'INSERER' LiFlC (i - 1) FlC;
  613. LiFlD = 'INSERER' LiFlD (i - 1) FlD;
  614. LiFl = 'INSERER' LiFl (i - 1) Fl;
  615. LiFlI = 'INSERER' LiFlI (i - 1) FlInt;
  616. FINSI;
  617. 'FIN' Bcl21;
  618.  
  619. EvFlC = 'EVOL' 'JAUN' 'MANU' 't (ans)' LiAbs2 'Debit_Conv' LiFlC;
  620. EvFlD = 'EVOL' 'BLEU' 'MANU' 't (ans)' LiAbs2 'Debit_Diff' LiFlD;
  621. EvFl = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit' LiFl;
  622. EvFlI = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit_int.' LiFlI;
  623.  
  624. * Tracé :
  625. * ------
  626. * Décalage des champs les uns par rapport aux autres :
  627. Decal1 = (1.2 * LXM) 0.;
  628. Decal2 = (2.4 * LXM) 0.;
  629.  
  630. * Tracé du champ de limite de solubilité :
  631. * ----------------------------------------
  632. SI (TrLimSol ET SoluL);
  633. 'TITRE' 'Limite de solubilite';
  634. MaSol2 = 'KCHA' MoTransp MaSol1 'CHAM';
  635. 'TRAC' MoTransp MaSol2 LTrame;
  636. FINSI;
  637.  
  638. * Tracé des profils de conc, prec, diss le long des lignes choisies :
  639. * -------------------------------------------------------------------
  640. SI TrLEvol;
  641.  
  642. IndLEvol = 'INDEX' LEvol;
  643. TabEvlC = table;
  644. SI Solub;
  645. TabEvlP = table;
  646. TabEvlL = table;
  647. Finsi;
  648.  
  649. 'REPETER' Bcl ('DIME' LEvol);
  650. * pour chaque ligne i où suivre les valeurs :
  651. i = IndLEvol.&Bcl;
  652.  
  653. * détermination du nom de la composante à suivre et évolution de
  654. * la limite de solubilité :
  655. si ('EGA' (LEvol.i.'P') 'CENTRE');
  656. mot1 = 'H';
  657. finsi;
  658. si ('EGA' (LEvol.i.'P') 'FACE');
  659. mot1 = 'TH';
  660. finsi;
  661. si ('EGA' (LEvol.i.'P') 'SOMMET');
  662. mot1 = 'SCAL';
  663. finsi;
  664.  
  665. * Evolution des concentrations :
  666. TitDesC = 'CHAINE' 'Concentration ' (LEvol.i.'T') ;
  667. TabEvlC.i = TRACHIS Transp 'CONCENTRATION' ' ' LiIndevl
  668. ('MOTS' mot1) ('MOTS' ' ') (LEvol.i.'M');
  669.  
  670. si (SoluL et ('NEG' (LEvol.i.'P') 'FACE'));
  671. * Adjonction de la limite de solubilité :
  672. si ('EGA' (LEvol.i.'P') 'CENTRE');
  673. EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSol1 'H' (LEvol.I.'M');
  674. finsi;
  675. si ('EGA' (LEvol.i.'P') 'SOMMET');
  676. EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSolS 'SCAL' (LEvol.I.'M');
  677. finsi;
  678. EvSatL = 'EXTR' EvSatL 'COUR' 1;
  679. nn = 'DIME' TabEvlC.i ;
  680. 'REPETER' bcl2 nn ;
  681. j = nn - &bcl2 + 1;
  682. TabEvlC.i.(j + 1) = TabEvlC.i.j;
  683. 'FIN' bcl2;
  684. TabEvlC.i.1 = table;
  685. TabEvlC.i.1 .'VALEUR' = EvSatL ;
  686. TabEvlC.i.1 .'LEGEND1' = 'Csat' ;
  687. TabEvlC.i.1 .'LEGEND2' = ' ';
  688. finsi;
  689.  
  690. DESTRA (TabEvlC.i) 'MIMA' 'AXES'
  691. 'TITR' TitDesC 'TITX' (LEvol.i.'L') 'TITY' 'C (mol/l)'
  692. 'YBOR' 0. CSat1 ;
  693.  
  694. SI (Solub et ('NEG' (LEvol.i.'P') 'FACE'));;
  695. * Evolution des précipités :
  696. TitDesP = 'CHAINE' 'Precipite ' (LEvol.i.'T') ;
  697. TabEvlP.i = TRACHIS Transp 'PRECIPITE' ' ' LiIndevl
  698. ('MOTS' mot1) ('MOTS' ' ') (LEvol.i.'M');
  699. DESTRA (TabEvlP.i) 'MIMA' 'AXES' 'TITR' TitDesP
  700. 'TITX' (LEvol.i.'L') 'TITY' 'P (mol/dm3)' ;
  701. * Evolution des dissolutions :
  702. TitDesL = 'CHAINE' 'Dissolution ' (LEvol.i.'T') ;
  703. TabEvlL.i = TRACHIS Transp 'DISSOLUTION' ' ' LiIndevl
  704. ('MOTS' mot1) ('MOTS' ' ') (LEvol.i.'M');
  705. DESTRA (TabEvlL.i) 'MIMA' 'AXES' 'TITR' TitDesL
  706. 'TITX' (LEvol.i.'L') 'TITY' 'L (mol/dm3/an)' ;
  707. Finsi;
  708.  
  709. 'FIN' bcl;
  710.  
  711. FINSI;
  712.  
  713. * Tracé des évolutions de conc, prec, diss aux points choisis :
  714. * -------------------------------------------------------------
  715. * (L'indice 0 est réservé à un point de la source)
  716. SI (TrPEvol ET (('DIME' ('INDEX' PEVol)) > 0));
  717.  
  718. IndPEvol = 'INDEX' PEvol;
  719. TabEvpC = table;
  720. SI Solub;
  721. TabEvpP = table;
  722. TabEvpL = table;
  723. Finsi;
  724.  
  725. tabrien = table;
  726. tabrien.1 = ' ' ;
  727.  
  728. 'REPETER' Bcl ('DIME' PEvol);
  729. * pour chaque point i où suivre les valeurs :
  730. i = IndPEvol.&Bcl;
  731.  
  732. * détermination du nom de la composante à suivre et évolution de
  733. * la limite de solubilité :
  734. si ('EGA' (PEvol.i.'P') 'CENTRE');
  735. mot1 = 'H';
  736. finsi;
  737. si ('EGA' (PEvol.i.'P') 'FACE');
  738. mot1 = 'TH';
  739. finsi;
  740. si ('EGA' (PEvol.i.'P') 'SOMMET');
  741. mot1 = 'SCAL';
  742. finsi;
  743.  
  744. * Evolution des concentrations :
  745. TitDesC = 'CHAINE' 'Concentration ' (PEvol.i.'T') ;
  746. TabEvpC.i = TRACHIT Transp 'CONCENTRATION' ' ' LiIndEvl
  747. ('MOTS' mot1) tabrien ('MANU' 'POI1' (PEvol.i.'M'));
  748.  
  749. si (SoluL et ('NEG' (PEvol.i.'P') 'FACE'));
  750. * Adjonction de la limite de solubilité :
  751. si ('EGA' (PEvol.i.'P') 'CENTRE');
  752. ValCSat = MaSol1 'EXTR' 'H' (PEvol.i.'M');
  753. finsi;
  754. si ('EGA' (PEvol.i.'P') 'SOMMET');
  755. ValCSat = MaSolS 'EXTR' 'SCAL' (PEvol.i.'M');
  756. finsi;
  757. EvSatP = 'EVOL' 'MANU' ('PROG' 0. TFin) ('PROG' 2*ValCSat);
  758. 'REPETER' bcl2 ('DIME' TabEvpC.i);
  759. j = ('DIME' TabEvpC.i) - &bcl2 + 1;
  760. TabEvpC.i.(j + 1) = TabEvpC.i.j;
  761. 'FIN' bcl2;
  762. TabEvpC.i.1 = table;
  763. TabEvpC.i.1 .'VALEUR' = EvSatP ;
  764. TabEvpC.i.1 .'LEGEND1' = 'Csat' ;
  765. TabEvpC.i.1 .'LEGEND2' = ' ';
  766. FINSI;
  767.  
  768. DESTRA (TabEvpC.i) 'MIMA' 'AXES'
  769. 'TITR' TitDesC 'TITX' 't (an)' 'TITY' 'C (mol/l)'
  770. 'YBOR' 0. CSat1 'XBOR'Tini TFin ;
  771.  
  772. SI (Solub et ('NEG' (PEvol.i.'P') 'FACE'));;
  773. * Evolution des précipités :
  774. TitDesP = 'CHAINE' 'Precipite ' (PEvol.i.'T') ;
  775. TabEvpP.i = TRACHIT Transp 'PRECIPITE' ' ' LiIndevl
  776. ('MOTS' mot1) tabrien ('MANU' 'POI1' (PEvol.i.'M'));
  777. DESTRA (TabEvpP.i) 'MIMA' 'AXES' 'TITR' TitDesP
  778. 'TITX' 't (an)' 'TITY' 'P (mol/dm3)'
  779. 'XBOR' Tini TFin ;
  780. * Evolution des dissolutions :
  781. TitDesL = 'CHAINE' 'Dissolution ' (PEvol.i.'T') ;
  782. TabEvpL.i = TRACHIT Transp 'DISSOLUTION' ' ' LiIndevl
  783. ('MOTS' mot1) tabrien ('MANU' 'POI1' (PEvol.i.'M'));
  784. DESTRA (TabEvpL.i) 'MIMA' 'AXES' 'TITR' TitDesL
  785. 'TITX' 't (an)' 'TITY' 'L (mol/dm3/an)'
  786. 'XBOR' Tini TFin ;
  787. Finsi;
  788.  
  789. 'FIN' bcl;
  790.  
  791. FINSI; COMM 'SI TrPEvol';
  792.  
  793. * Tracé des champs de conc, prec, aux différents instants :
  794. * ---------------------------------------------------------
  795. * On déplace les supports des chpos pour tracer les fractions de
  796. * précipité à côté des concentrations :
  797. TrPrec = TrPrec ET Solub;
  798.  
  799. *-- On trace pour chaque temps choisi :
  800. SI (TrConc OU TrPrec);
  801.  
  802. *-- Echelles :
  803. SI TrConc;
  804. SI SoluL;
  805. CMinD = 'MINIMUM' MaSol1;
  806. CMaxD = 'MAXIMUM' MaSol1;
  807. SINON;
  808. CMinD = 0.;
  809. CMaxD = 0.;
  810. FINSI;
  811. FINSI;
  812.  
  813. 'REPETER' Bcl17 NTSauv;
  814. i = &Bcl17;
  815. Id = IndTSauv.i;
  816. SI TrConc;
  817. CMinD = 'MINIMUM'
  818. ('PROG' CMinD ('MINIMUM' Transp.Concentration.Id));
  819. CMaxD = 'MAXIMUM'
  820. ('PROG' CMaxD ('MAXIMUM' Transp.Concentration.Id));
  821. FINSI;
  822. 'FIN' Bcl17;
  823.  
  824. SI TrConc;
  825. CMinD = 'MINIMUM' ('PROG' CMinD (CMaxD * 1.e-2));
  826. CMinD = 'MAXIMUM' ('PROG' CMinD 0.);
  827. EchC0 = 'PROG' 0. 'PAS' (CMaxD / 13.) CMaxD;
  828. FINSI;
  829.  
  830. 'REPETER' Bcl3 NTGrph;
  831. Id = 'EXTR' LiIndgph &Bcl3;
  832. Tps = Transp.Temps.Id;
  833. TitGraph = ' ' ;
  834. SI TrConc;
  835. TitGraph = 'CHAINE' TitGraph 'Solute ';
  836. FINSI;
  837. SI (TrConc ET TrPrec);
  838. TitGraph = 'CHAINE' TitGraph 'et ';
  839. FINSI;
  840. SI TrPrec;
  841. TitGraph = 'CHAINE' TitGraph 'Precipite ';
  842. FINSI;
  843. TitGraph = 'CHAINE' TitGraph 'a t=' Tps ' ans ';
  844. 'TITRE' TitGraph;
  845.  
  846. SI TrConc;
  847. Conc2 = 'KCHA' MoTransp Transp.Concentration.Id 'CHAM';
  848. MoTot = MoTransp;
  849. Tot2 = Conc2;
  850. FINSI;
  851. SI TrPrec;
  852. Prec1 = 'KCHA' MoTransp Transp.Precipite.Id 'CHAM';
  853. MoTra2 Prec2 = MoTransp Prec1 'PLUS' Decal1;
  854. LTra2 = LTrame 'PLUS' Decal1;
  855. 'ELIMINATION' (ltra2 'ET' ('EXTR' Prec2 'MAIL')) epselim;
  856. SI (NON TrConc);
  857. Tot2 = Prec2;
  858. MoTot = MoTra2;
  859. SINON;
  860. Tot2 = Tot2 ET Prec2;
  861. MoTot = MoTot ET MoTra2;
  862. FINSI;
  863. FINSI;
  864.  
  865. SI (('MINIMUM' Tot2) < (('MAXIMUM' Tot2) - EpsPrec));
  866. SI (TrConc ET (NON TrPrec));
  867. SI (('MINIMUM' Conc2) < (('MAXIMUM' Conc2) - EpsPrec));
  868. 'TRACER' Motot conc2 echc0 LTrame;
  869. FINSI;
  870. SINON;
  871. SI TrConc;
  872. 'TRACER' MoTot Tot2 echc0;
  873. SINON;
  874. 'TRACER' MoTot Tot2 echc0 ltra2;
  875. FINSI;
  876. FINSI;
  877. FINSI;
  878. 'FIN' Bcl3;
  879. FINSI;
  880.  
  881. * Tracé des champs de dissolution :
  882. * ---------------------------------
  883. TrLib = TrLib ET Solub;
  884. SI TrLib;
  885. 'REPETER' Bcl22 (NTgrph - 1);
  886. Id = ('EXTR' LiIndgph (&Bcl22 + 1)) - 1 ;
  887. Tps = Transp.Temps.Id;
  888. Id2 = Id + 1;
  889. Tp2 = Transp.Temps.Id2;
  890. dt = Tp2 - Tps;
  891. TitGraph = 'CHAINE' 'Quantite dissoute entre '
  892. Tps ' et ' Tp2 ' annees (mol/m3/an)';
  893. 'TITRE' TitGraph;
  894. Libe2 = Transp.Dissolution.Id2;
  895. Libe3 = 'KCHA' MoTransp Libe2 'CHAM';;
  896. SI ((('MAXIMUM' Libe3) - ('MINIMUM' Libe3)) > EpsPrec);
  897. 'TRAC' MoTransp Libe3 Ltrame ;
  898. FINSI;
  899. 'FIN' Bcl22;
  900. FINSI;
  901.  
  902. * Tracé du flux traversant la ligne choisie :
  903. * ------------------------------------------
  904. SI TrFlux ;
  905.  
  906. TitDess = 'CHAINE' 'Flux traversant ' DFlux.'TITRE' ;
  907. TabDess = 'TABLE';
  908. TabDess.1 = 'MARQ CARR' ;
  909. TabDess.2 = 'MARQ TRIA' ;
  910. TabDess.3 = 'MARQ PLUS' ;
  911. TabDess.'TITRE' = 'TABLE';
  912. TabDess.'TITRE' . 1 = 'Debit_Convectif';
  913. TabDess.'TITRE' . 2 = 'Debit_Diffusif';
  914. TabDess.'TITRE' . 3 = 'Debit_Total';
  915. EvTot = EvFlC ET EvFlD ET EvFl;
  916. 'DESSIN' EvTot 'MIMA' TabDess 'TITR' TitDess 'TITX' 't (an)' 'TITY'
  917. 'mol/l/an' 'LEGE' 'AXES' 'XBOR' TIni TFin;
  918. FINSI;
  919.  
  920. * ===========================
  921. * | TEST DE BON DEROULEMENT |
  922. * ===========================
  923. * Le test se fait en comparant le flux total traversant la ligne avec
  924. * une liste obtenue par un calcul préalable.
  925. *liflsav = 'PROG' .68457 1.1205 1.3157 1.3279 1.2973 1.1393 .79427
  926. * .48185 .34038 .24040 ;
  927.  
  928. liflsav = 'PROG' 5.69132E-02 .10321 .11576 .10992 9.835E-02 8.6297E-02
  929. 7.53200E-02 6.57577E-02 5.75179E-02 5.04195E-02 ;
  930.  
  931. *'LISTE' lifl;
  932.  
  933. err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav);
  934. mess 'Erreur = ' err1 ;
  935.  
  936. si (err1 > CRIT1);
  937. 'ERREUR' 5;
  938. 'OPTION' donn 5 ;
  939. sinon;
  940. 'ERREUR' 0;
  941. finsi;
  942.  
  943. Fin;
  944.  
  945.  
  946.  
  947.  
  948.  
  949.  
  950.  
  951.  
  952.  
  953.  
  954.  
  955.  
  956.  

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