Télécharger transport2EFMH.dgibi

Retour à la liste

Numérotation des lignes :

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

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