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' MoHydro QFHydro ;
  272.  
  273.  
  274. * Points où l'on suivra la concentration :
  275. * ----------------------------------------
  276. PEvol = 'TABLE';
  277. * Centre du combustible :
  278. PEvol.1 = 'TABLE';
  279. PEvol.1 .'T' = 'CHAINE' 'dans la source' ;
  280. PEvol.1 .'P' = 'CENTRE';
  281. PEvol.1 .'M' = ( (.5 * LX1) ((LY2 + LY1) / 2.) );
  282. PEvol.1 .'M' = CeHydro 'POINT' 'PROC' (PEvol.1 .'M');
  283. * Juste après la transition aval :
  284. PEvol.2 = 'TABLE';
  285. PEvol.2 .'T' = 'CHAINE' 'post-transition aval' ;
  286. PEvol.2 .'P' = 'CENTRE';
  287. PEvol.2 .'M' = CeHydro 'POINT'
  288. 'PROC' (P5 'PLUS' (0.05 0.05));
  289.  
  290. * Lignes où capturer les évolutions des concentrations :
  291. * ------------------------------------------------------
  292. LEvol = 'TABLE';
  293.  
  294. 'REPETER' Bcl8 ('NBEL' LGau);
  295. pt = CeHydro 'POINT' 'PROC'
  296. ((LGau 'POINT' &Bcl8) 'PLUS' (0 0.01));
  297. SI (&Bcl8 > 2);
  298. LGauC = LGauC ET ('MANU' 'SEG2' pt0 pt);
  299. FINSI;
  300. SI (&Bcl8 'EGA' 2);
  301. LGauC = ('MANU' 'SEG2' pt0 pt);
  302. FINSI;
  303. pt0 = pt;
  304. 'FIN' Bcl8;
  305. 'OUBLIER' pt0;
  306. 'OUBLIER' pt;
  307.  
  308. LEvol.1 = 'TABLE';
  309. LEvol.1 .M = 'COULEUR' LGauC 'BLAN';
  310. LEvol.1 .P = 'CENTRE';
  311. LEvol.1 .L = 'Y (m)';
  312. LEvol.1 .T = 'CHAINE' 'longitudinale' ;
  313.  
  314. * =============
  315. * | TRANSPORT |
  316. * =============
  317. * Modèle propre au transport :
  318. * ----------------------------
  319. MoTransp = 'MODELE' QMtot DARCY ANISOTROPE ;
  320.  
  321. * Matrice Porosité cinématique (Type ChamElem) :
  322. * ----------------------------------------------
  323. CeTransp = 'DOMA' MoTransp 'CENTRE' ;
  324. MaPor1 = 'MANU' 'CHPO' CeTransp 1 'CK' ValPor NATURE DIFFUS;
  325. MaPor = MaPor1;
  326.  
  327. * Porosité * volume de maille
  328. * ----------------------------
  329. MaPorVol= VolChPro * ('EXCO' 'CK' MaPor1 'SCAL');
  330.  
  331. * Temps :
  332. * -------
  333. * Temps à calculer
  334. LiTCalc = 'PROG' TIni 'PAS' DeltaT Tfin;
  335. NTCalc = 'DIME' LiTCalc;
  336. 'REMPLACER' LiTCalc NTCalc Tfin;
  337.  
  338. * Ecriture de ces temps sous forme de table ;
  339. TabTCalc = 'TABLE';
  340. 'REPETER' Bcl26 NTCalc;
  341. Id = &Bcl26 - 1;
  342. Tps = 'EXTR' LiTCalc &Bcl26;
  343. TabTCalc.Id = Tps;
  344. 'FIN' Bcl26;
  345.  
  346. * Temps sauvegardés
  347. LiTSauv = LiTCalc ;
  348. NTSauv = NTCalc ;
  349.  
  350. * Temps de tracage des champs
  351. LiIndGph = 'LECT' 0 5 10 ;
  352. NTGrph = 'DIME' LiIndGph;
  353.  
  354. * Temps de tracage des évolutions
  355. LiTEvol = litsauv ;
  356. liindEvl = 'LECT' 0 pas 1 (ntsauv - 1) ;
  357. NTEvol = 'DIME' litevol;
  358.  
  359. * Temps servant à la création des chargements
  360. TSup = 1.2 * TFin;
  361. DelTDef = 0.001 * TSup;
  362.  
  363. * ====================================
  364. * | CONSTRUCTION DES DONNEES DE BASE |
  365. * ====================================
  366. * Dissolution / Libération ;
  367. * --------------------------
  368. *-- Matrice Limite de solubilité ;
  369. * Variation exponentielle décroissante de la source vers l'aval :
  370. XX YY = 'COOR' CeTransp ;
  371. Sup2 = 'POINT' YY 'EGINF' LY2;
  372. MaSol2 = 'MANU' 'CHPO' Sup2 1 'H' CSat1 NATURE DIFFUS;
  373. Sup3 = 'POINT' YY 'SUPERIEUR' LY2;
  374. YY3 = 'REDU' (YY - LY2) Sup3;
  375. a = ('LOG' (CSat2 / CSat1)) / (LYM - LY2);
  376. MaSol3 = CSat1 * ('EXP' (a * YY3));
  377. MaSol3 = 'NOMC' MaSol3 'H' NATURE DIFFUS;
  378. MaSol1 = MaSol2 ET MaSol3;
  379.  
  380. * Matériau associé au transport :
  381. * -------------------------------
  382. * Terme diffusif pur (diffusivité effective) = Porosité * Diffusivité :
  383. DifEff = DifMolec * FacForme * MaPor1;
  384. MaTransp = NOMC 'K' DifEff ;
  385.  
  386. * Conditions Initiales :
  387. * ----------------------
  388.  
  389. *-- Concentration aux centres des éléments (sert à définir la quantité
  390. * initiale d'élément) :
  391. Conc01 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET';
  392. Conc02 = 'MANU' 'CHPO'( 'DOMA' MoActif 'CENTRE') 1 'H' CActif0
  393. 'NATURE' 'DISCRET';
  394. Conc0 = Conc01 ET Conc02;
  395. Conc0 = 'NOMC'
  396. ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Conc0 'SCAL')) 'H';
  397.  
  398. *-- Flux diffusif initial (aux centres des faces) - décoratif puisque la
  399. * diffusion est implicite :
  400. QFDif0 = 'MANU' 'CHPO' FaHydro 1 'FLUX' 0. 'NATURE' 'DISCRET';
  401.  
  402. *-- Précipité :
  403. Prec01 = 'MANU' 'CHPO' ( 'DOMA' MoActif 'CENTRE') 1 'H'
  404. PActif0 'NATURE' 'DISCRET';
  405. Prec02 = 'MANU' 'CHPO' CeHydro 1 'H' 0. 'NATURE' 'DISCRET';
  406. Prec0 = Prec01 ET Prec02;
  407. Prec0 = 'NOMC'
  408. ('KCHT' MoHydro 'SCAL' 'CENTRE' ('NOMC' Prec0 'SCAL')) 'H';
  409.  
  410.  
  411. * Conditions aux limites (type CHARGEMENT) :
  412. * ------------------------------------------
  413. * Concentration imposée en amont et en aval, flux nul ailleurs.
  414. Ceinf = 'DOMA' Moinf 'CENTRE' ;
  415. Cesup = 'DOMA' Mosup 'CENTRE' ;
  416. CCL = 'MANU' 'CHPO' (Ceinf et Cesup) 1 'H' 0. 'NATURE' 'DISCRET';
  417. CharCL = CCL ;
  418. ChgCL = 'CHARG' CharCL ('EVOL' 'MANU' ('PROG' 0. TSup) ('PROG' 1. 1.));
  419.  
  420. * ==============
  421. * | Résolution |
  422. * ==============
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445. *-- Table de transport :
  446. Transp = 'TABLE';
  447. Transp.'TEMPS' = 'TABLE';
  448. Transp.'TRACE_CONC' = 'TABLE';
  449. Transp.'CONCENTRATION' = 'TABLE';
  450. Transp.'FLUX' = 'TABLE';
  451. Transp.'PRECIPITE' = 'TABLE';
  452. Transp.'DISSOLUTION' = 'TABLE';
  453.  
  454.  
  455. Transp.'SOUSTYPE' = 'DARCY_TRANSITOIRE';
  456. Transp.'MODELE' = MoTransp;
  457. Transp.'CARACTERISTIQUES' = MaTransp;
  458. Transp.'POROSITE' = MaPor;
  459. Transp.'COEF_RETARD' = 10.D0;
  460. *Transp.'CONVECTION' = QFHydro;
  461. Transp.'CONVECTION' = (nomc SCAL (1.D0 * QFHydro)) /
  462. (doma MOTRANSP SURFACE) ;
  463. Transp . 'CONVECTION' = Transp . CONVECTION *
  464. (doma MoTRANSP NORMALE);
  465. Transp . 'VITELEM' = 'HVIT' MOtransp QFHYDRO ;
  466.  
  467.  
  468.  
  469. Transp.'LIMITE_SOLUBILITE' = MaSol1;
  470. Transp.'COEF_DISSOLUTION' = 1.D8;
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  
  477. * Conditions initiales :
  478. Transp.'TEMPS'. 0 = TIni;
  479. Transp.'CONCENTRATION'. 0 = Conc0;
  480. Transp . 'FLUXDIFF' = table;
  481. Transp.'FLUXDIFF' . 0 = 0.D0 * QFDif0;
  482. Transp . 'FLUXCONV' = table;
  483. Transp.'FLUXCONV' . 0 = 0.D0 * QFDif0;
  484.  
  485.  
  486. Transp.'PRECIPITE'. 0 = Prec0;
  487. * Transp.'DISSOLUTION'. 0 = Libe0;
  488.  
  489.  
  490.  
  491.  
  492.  
  493. * Conditions aux limites :
  494. Transp.'TRACE_IMPOSE' = ChgCL;
  495.  
  496. * Paramètres numériques :
  497. Transp.'THETA_DIFF' = TetaDif;
  498. Transp.'THETA_CONVECTION' = TetaCnv;
  499. Transp.'LUMP' = VRAI;
  500. Transp.'TYPDISCRETISATION' = 'EFMH';
  501. Transp.'DECENTR' = VRAI;
  502.  
  503. TABRES = table METHINV;
  504. TABRES . 'TYPINV' = 1;
  505. TABRES . 'PRECOND' = 3;
  506.  
  507. Transp . 'METHINV' = TABRES;
  508. *ransp.'THETA_DISS' = TetaDis;
  509.  
  510.  
  511.  
  512. SI ('EGA' Methode 'Point fixe Dissolution');
  513. Transp.'EPSI_LIM' = EpsSolub;
  514. Transp.'ITMAX_LIM' = NbItMax;
  515. FINSI;
  516. Transp.'TEMPS_CALCULES' = LiTCalc;
  517. Transp.'TEMPS_SAUVES' = LiTSauv;
  518.  
  519. * ==========
  520. * | CALCUL |
  521. * ==========
  522.  
  523. TRANSGEN Transp;
  524.  
  525. * ===================
  526. * | POST-TRAITEMENT |
  527. * ===================
  528.  
  529. MESS ' ';
  530. MESS 'Post-Traitement :';
  531. MESS '-----------------';
  532. MESS ' ';
  533.  
  534. IndTSauv = 'INDEX' (Transp.temps);
  535.  
  536. * Préparation tracé :
  537. * -------------------
  538. * Calcul du flux à travers une certaine ligne :
  539. * ---------------------------------------------
  540. * Les propriétés de la ligne sont stockées dans DFlux et NleFlux.
  541. * Le nom de la ligne doit etre mis dans DFluxTIT .
  542. NORHydro = 'DOMA' MoTransp 'NORMALE' ;
  543. Ori2 = 'PSCAL' ('REDU' NORHydro CeFlux ) NleFlux
  544. ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY');
  545. FlInt = 0.;
  546. liabs2 = 'ENLEVER' litcalc 1;
  547. 'REPETER' Bcl21 (NTSauv - 1);
  548. i = &Bcl21 + 1;
  549. Id = IndTSauv.i;
  550. Tps = Transp.Temps.Id;
  551. Id1 = IndTSauv.(i - 1);
  552. Tp1 = Transp.Temps.(id - 1);
  553. dt = Tps - Tp1;
  554.  
  555. *-- Flux diffusif et convectif à l'instant d'indice Id
  556. * (diffusion implicite TetaDif et convection implicite TetaCnv) :
  557. * méthode par produit scalaire flux * normale sur 1 ligne
  558.  
  559. * Débit diffusif sortant :
  560. DebD1 = ((1. - TetaDif) * Transp.FLUXDIFF.Id1)
  561. + (TetaDif * Transp.FLUXDIFF.Id);
  562. DebD2 = 'REDU' DebD1 CeFlux ;
  563. DebD3 = 'RESULT' (DebD2 * Ori2);
  564. FlD = ('MAXIMUM' DebD3);
  565. Fl = FlD;
  566.  
  567. * Débit convectif sortant :
  568. THCnv = ((1. - TetaCnv) * Transp.FLUXCONV.Id1)
  569. + (TetaCnv * Transp.FLUXCONV.Id);
  570. DebC1 = THCnv ;
  571. DebC2 = 'REDU' DebC1 CeFlux ;
  572. DebC3 = 'RESULT' (DebC2 * Ori2);
  573. FlC = ('MAXIMUM' DebC3);
  574. Fl = FlC + FlD;
  575.  
  576. FlInt = FlInt + (Fl * dt);
  577.  
  578. SI (i 'EGA' 2);
  579. LiFlC = 'PROG' FlC;
  580. LiFlD = 'PROG' FlD;
  581. LiFl = 'PROG' Fl;
  582. LiFlI = 'PROG' FlInt;
  583. SINON;
  584. LiFlC = 'INSERER' LiFlC (i - 1) FlC;
  585. LiFlD = 'INSERER' LiFlD (i - 1) FlD;
  586. LiFl = 'INSERER' LiFl (i - 1) Fl;
  587. LiFlI = 'INSERER' LiFlI (i - 1) FlInt;
  588. FINSI;
  589. 'FIN' Bcl21;
  590.  
  591. EvFlC = 'EVOL' 'JAUN' 'MANU' 't (ans)' LiAbs2 'Debit_Conv' LiFlC;
  592. EvFlD = 'EVOL' 'BLEU' 'MANU' 't (ans)' LiAbs2 'Debit_Diff' LiFlD;
  593. EvFl = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit' LiFl;
  594. EvFlI = 'EVOL' 'VERT' 'MANU' 't (ans)' LiAbs2 'Debit_int.' LiFlI;
  595.  
  596. * Tracé :
  597. * ------
  598. * Décalage des champs les uns par rapport aux autres :
  599. Decal1 = (1.2 * LXM) 0.;
  600. Decal2 = (2.4 * LXM) 0.;
  601.  
  602.  
  603.  
  604. * Tracé des profils de conc, prec, diss le long des lignes choisies :
  605. * -------------------------------------------------------------------
  606. SI TrLEvol;
  607.  
  608. IndLEvol = 'INDEX' LEvol;
  609. TabEvlC = table;
  610. TabEvlP = table;
  611. TabEvlL = table;
  612.  
  613. 'REPETER' Bcl ('DIME' LEvol);
  614. * pour chaque ligne i où suivre les valeurs :
  615. i = IndLEvol.&Bcl;
  616.  
  617. * détermination du nom de la composante à suivre et évolution de
  618. * la limite de solubilité :
  619. si ('EGA' (LEvol.i.'P') 'CENTRE');
  620. mot1 = 'H';
  621. finsi;
  622. si ('EGA' (LEvol.i.'P') 'FACE');
  623. mot1 = 'TH';
  624. finsi;
  625. si ('EGA' (LEvol.i.'P') 'SOMMET');
  626. mot1 = 'SCAL';
  627. finsi;
  628.  
  629. * Evolution des concentrations :
  630. TitDesC = 'CHAINE' 'Concentration ' (LEvol.i.'T') ;
  631. TabEvlC.i = TRACHIS Transp 'CONCENTRATION' (LEvol.i.'M') LiIndevl
  632. ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  633.  
  634. si ( ('NEG' (LEvol.i.'P') 'FACE'));
  635. * Adjonction de la limite de solubilité :
  636. si ('EGA' (LEvol.i.'P') 'CENTRE');
  637. EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSol1 'H' (LEvol.I.'M');
  638. finsi;
  639. si ('EGA' (LEvol.i.'P') 'SOMMET');
  640. EvSatL = 'EVOL' 'BLAN' 'CHPO' MaSolS 'SCAL' (LEvol.I.'M');
  641. finsi;
  642. EvSatL = 'EXTR' EvSatL 'COUR' 1;
  643. 'REPETER' bcl2 ('DIME' TabEvlC.i);
  644. j = ('DIME' TabEvlC.i) - &bcl2 + 1;
  645. TabEvlC.i.(j + 1) = TabEvlC.i.j;
  646. 'FIN' bcl2;
  647. TabEvlC.i.1 = table;
  648. TabEvlC.i.1 .'VALEUR' = EvSatL ;
  649. TabEvlC.i.1 .'LEGEND1' = 'Csat' ;
  650. TabEvlC.i.1 .'LEGEND2' = ' ';
  651. finsi;
  652.  
  653. * DESTRA (TabEvlC.i) 'MIMA' 'AXES'
  654. * 'TITR' TitDesC 'TITX' (LEvol.i.'L') 'TITY' 'C (mol/l)'
  655. * 'YBOR' 0. CSat1 ;
  656.  
  657. SI ( ('NEG' (LEvol.i.'P') 'FACE'));;
  658. * Evolution des précipités :
  659. TitDesP = 'CHAINE' 'Precipite ' (LEvol.i.'T') ;
  660. TabEvlP.i = TRACHIS Transp 'PRECIPITE' (LEvol.i.'M') LiIndevl
  661. ('MOTS' mot1) ('MOTS' ' ') 'PREF' ' ' ;
  662. DESTRA (TabEvlP.i) 'MIMA' 'AXES' 'TITR' TitDesP
  663. 'TITX' (LEvol.i.'L') 'TITY' 'P (mol/dm3)' ;
  664. Finsi;
  665.  
  666. 'FIN' bcl;
  667.  
  668. FINSI;
  669.  
  670. * Tracé des évolutions de conc, prec, diss aux points choisis :
  671. * -------------------------------------------------------------
  672. * (L'indice 0 est réservé à un point de la source)
  673. SI (TrPEvol ET (('DIME' ('INDEX' PEVol)) > 0));
  674.  
  675. IndPEvol = 'INDEX' PEvol;
  676. TabEvpC = table;
  677. TabEvpP = table;
  678. TabEvpL = table;
  679.  
  680. 'REPETER' Bcl ('DIME' PEvol);
  681. * pour chaque point i où suivre les valeurs :
  682. i = IndPEvol.&Bcl;
  683.  
  684. * détermination du nom de la composante à suivre et évolution de
  685. * la limite de solubilité :
  686. si ('EGA' (PEvol.i.'P') 'CENTRE');
  687. mot1 = 'H';
  688. finsi;
  689. si ('EGA' (PEvol.i.'P') 'FACE');
  690. mot1 = 'TH';
  691. finsi;
  692. si ('EGA' (PEvol.i.'P') 'SOMMET');
  693. mot1 = 'SCAL';
  694. finsi;
  695.  
  696. * Evolution des concentrations :
  697. TitDesC = 'CHAINE' 'Concentration ' (PEvol.i.'T') ;
  698. TabEvpC.i = TRACHIT Transp 'CONCENTRATION'
  699. ('MANU' 'POI1' (PEvol.i.'M')) LiIndEvl
  700. ('MOTS' mot1) 'PREF' ' ' ;
  701.  
  702. si (('NEG' (PEvol.i.'P') 'FACE'));
  703. * Adjonction de la limite de solubilité :
  704. si ('EGA' (PEvol.i.'P') 'CENTRE');
  705. ValCSat = MaSol1 'EXTR' 'H' (PEvol.i.'M');
  706. finsi;
  707. si ('EGA' (PEvol.i.'P') 'SOMMET');
  708. ValCSat = MaSolS 'EXTR' 'SCAL' (PEvol.i.'M');
  709. finsi;
  710. EvSatP = 'EVOL' 'MANU' ('PROG' 0. TFin) ('PROG' 2*ValCSat);
  711. 'REPETER' bcl2 ('DIME' TabEvpC.i);
  712. j = ('DIME' TabEvpC.i) - &bcl2 + 1;
  713. TabEvpC.i.(j + 1) = TabEvpC.i.j;
  714. 'FIN' bcl2;
  715. TabEvpC.i.1 = table;
  716. TabEvpC.i.1 .'VALEUR' = EvSatP ;
  717. TabEvpC.i.1 .'LEGEND1' = 'Csat' ;
  718. TabEvpC.i.1 .'LEGEND2' = ' ';
  719. FINSI;
  720.  
  721. * DESTRA (TabEvpC.i) 'MIMA' 'AXES'
  722. * 'TITR' TitDesC 'TITX' 't (an)' 'TITY' 'C (mol/l)'
  723. * 'YBOR' 0. CSat1 'XBOR' Tini TFin ;
  724.  
  725. SI ( ('NEG' (PEvol.i.'P') 'FACE'));;
  726. * Evolution des précipités :
  727. TitDesP = 'CHAINE' 'Precipite ' (PEvol.i.'T') ;
  728. TabEvpP.i = TRACHIT Transp 'PRECIPITE'
  729. ('MANU' 'POI1' (PEvol.i.'M')) LiIndevl
  730. ('MOTS' mot1) 'PREF' ' ' ;
  731. * DESTRA (TabEvpP.i) 'MIMA' 'AXES' 'TITR' TitDesP
  732. * 'TITX' 't (an)' 'TITY' 'P (mol/dm3)'
  733. * 'XBOR' Tini TFin ;
  734. Finsi;
  735.  
  736. 'FIN' bcl;
  737.  
  738. FINSI; COMM 'SI TrPEvol';
  739.  
  740. * Tracé des champs de conc, prec, aux différents instants :
  741. * ---------------------------------------------------------
  742. * On déplace les supports des chpos pour tracer les fractions de
  743. * précipité à côté des concentrations :
  744. TrPrec = TrPrec;
  745.  
  746. *-- On trace pour chaque temps choisi :
  747. SI (TrConc OU TrPrec);
  748.  
  749. *-- Echelles :
  750. SI TrConc;
  751. CMinD = 'MINIMUM' MaSol1;
  752. CMaxD = 'MAXIMUM' MaSol1;
  753. SINON;
  754. CMinD = 0.;
  755. CMaxD = 0.;
  756. FINSI;
  757. FINSI;
  758.  
  759. 'REPETER' Bcl17 NTSauv;
  760. i = &Bcl17;
  761. Id = IndTSauv.i;
  762. SI TrConc;
  763. CMinD = 'MINIMUM'
  764. ('PROG' CMinD ('MINIMUM' Transp.Concentration.Id));
  765. CMaxD = 'MAXIMUM'
  766. ('PROG' CMaxD ('MAXIMUM' Transp.Concentration.Id));
  767. FINSI;
  768. 'FIN' Bcl17;
  769.  
  770. SI TrConc;
  771. CMinD = 'MINIMUM' ('PROG' CMinD (CMaxD * 1.e-2));
  772. CMinD = 'MAXIMUM' ('PROG' CMinD 0.);
  773. EchC0 = 'PROG' 0. 'PAS' (CMaxD / 13.) CMaxD;
  774. FINSI;
  775.  
  776. 'REPETER' Bcl3 NTGrph;
  777. Id = 'EXTR' LiIndgph &Bcl3;
  778. Tps = Transp.Temps.Id;
  779. TitGraph = ' ' ;
  780. SI TrConc;
  781. TitGraph = 'CHAINE' TitGraph 'Solute ';
  782. FINSI;
  783. SI (TrConc ET TrPrec);
  784. TitGraph = 'CHAINE' TitGraph 'et ';
  785. FINSI;
  786. SI TrPrec;
  787. TitGraph = 'CHAINE' TitGraph 'Precipite ';
  788. FINSI;
  789. TitGraph = 'CHAINE' TitGraph 'a t=' Tps ' ans ';
  790. 'TITRE' TitGraph;
  791.  
  792. Conc2 = 'KCHA' MoTransp Transp.Concentration.Id 'CHAM';
  793. MoTot = MoTransp;
  794. Tot2 = Conc2;
  795. SI TrPrec;
  796. Prec1 = 'KCHA' MoTransp Transp.Precipite.Id 'CHAM';
  797. MoTra2 Prec2 = MoTransp Prec1 'PLUS' Decal1;
  798. LTra2 = LTrame 'PLUS' Decal1;
  799. 'ELIMINATION' (ltra2 'ET' ('EXTR' Prec2 'MAIL')) epselim;
  800. SI (NON TrConc);
  801. Tot2 = Prec2;
  802. MoTot = MoTra2;
  803. SINON;
  804. Tot2 = Tot2 ET Prec2;
  805. MoTot = MoTot ET MoTra2;
  806. FINSI;
  807. FINSI;
  808.  
  809. SI (('MINIMUM' Tot2) < (('MAXIMUM' Tot2) - EpsPrec));
  810. SI (TrConc ET (NON TrPrec));
  811. SI (('MINIMUM' Conc2) < (('MAXIMUM' Conc2) - EpsPrec));
  812. 'TRACER' Motot conc2 echc0 LTrame;
  813. FINSI;
  814. SINON;
  815. SI TrConc;
  816. 'TRACER' MoTot Tot2 echc0;
  817. SINON;
  818. 'TRACER' MoTot Tot2 echc0 ltra2 ;
  819. FINSI;
  820. FINSI;
  821. FINSI;
  822. 'FIN' Bcl3;
  823. FINSI;
  824.  
  825. * Tracé des champs de dissolution :
  826. * ---------------------------------
  827. TrLib = TrLib;
  828. SI TrLib;
  829. 'REPETER' Bcl22 (NTgrph - 1);
  830. Id = ('EXTR' LiIndgph (&Bcl22 + 1)) - 1 ;
  831. Tps = Transp.Temps.Id;
  832. Id2 = Id + 1;
  833. Tp2 = Transp.Temps.Id2;
  834. dt = Tp2 - Tps;
  835. TitGraph = 'CHAINE' 'Quantite dissoute entre '
  836. Tps ' et ' Tp2 ' annees (mol/m3/an)';
  837. 'TITRE' TitGraph;
  838. 'FIN' Bcl22;
  839. FINSI;
  840.  
  841. * Tracé du flux traversant la ligne choisie :
  842. * ------------------------------------------
  843. SI TrFlux ;
  844.  
  845. TitDess = 'CHAINE' 'Flux traversant ' DFluxTIT ;
  846. TabDess = 'TABLE';
  847. TabDess.1 = 'MARQ CARR' ;
  848. TabDess.2 = 'MARQ TRIA' ;
  849. TabDess.3 = 'MARQ PLUS' ;
  850. TabDess.'TITRE' = 'TABLE';
  851. TabDess.'TITRE' . 1 = 'Debit_Convectif';
  852. TabDess.'TITRE' . 2 = 'Debit_Diffusif';
  853. TabDess.'TITRE' . 3 = 'Debit_Total';
  854. EvTot = EvFlC ET EvFlD ET EvFl;
  855. 'DESSIN' EvTot 'MIMA' TabDess 'TITR' TitDess 'TITX' 't (an)' 'TITY'
  856. 'mol/l/an' 'LEGE' 'AXES' 'XBOR' TIni TFin;
  857. FINSI;
  858.  
  859. * ===========================
  860. * | TEST DE BON DEROULEMENT |
  861. * ===========================
  862. * Le test se fait en comparant le flux total traversant la ligne avec
  863. * une liste obtenue par un calcul préalable.
  864. liflsav = 'PROG' 0.56061 1.0216 1.2444 1.3013 1.2365 1.0595 0.93058
  865. 0.64840 0.41501 0.29818;
  866.  
  867. err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav);
  868. 'MESSAGE' 'difference avec sol de ref' err1;
  869.  
  870. si (err1 > CRIT1);
  871. 'ERREUR' 5;
  872. sinon;
  873. 'ERREUR' 0;
  874. finsi;
  875.  
  876. 'LISTE' crit1;
  877.  
  878. Fin;
  879.  
  880.  
  881.  
  882.  
  883.  
  884.  
  885.  
  886.  
  887.  
  888.  
  889.  
  890.  
  891.  

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