Télécharger la plaquette Cast3M


fr en

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.     'OPTION' ISOV SURF;
  809.  
  810.     SI (('MINIMUM' Tot2) < (('MAXIMUM' Tot2) - EpsPrec));
  811.       SI (TrConc ET (NON TrPrec));
  812.         SI (('MINIMUM' Conc2) < (('MAXIMUM' Conc2) - EpsPrec));
  813.           'TRACER' Motot conc2 echc0 LTrame;
  814.         FINSI;
  815.       SINON;
  816.         SI TrConc;
  817.           'TRACER' MoTot Tot2 echc0;
  818.         SINON;
  819.           'TRACER' MoTot Tot2 echc0 ltra2 ;
  820.         FINSI;
  821.       FINSI;
  822.     FINSI;
  823.   'FIN' Bcl3;
  824. FINSI;
  825.  
  826. * Tracé des champs de dissolution :
  827. * ---------------------------------
  828. TrLib = TrLib;
  829. SI TrLib;
  830.   'REPETER' Bcl22 (NTgrph - 1);
  831.     Id  = ('EXTR' LiIndgph (&Bcl22 + 1)) - 1 ;
  832.     Tps = Transp.Temps.Id;
  833.     Id2 = Id + 1;
  834.     Tp2 = Transp.Temps.Id2;
  835.     dt  = Tp2 -  Tps;
  836.     TitGraph = 'CHAINE' 'Quantite dissoute entre ' 
  837.                Tps ' et ' Tp2 ' annees (mol/m3/an)';
  838.     'TITRE' TitGraph;
  839.   'FIN' Bcl22;
  840. FINSI;
  841.  
  842. * Tracé du flux traversant la ligne choisie :
  843. * ------------------------------------------
  844. SI TrFlux ;
  845.  
  846.   TitDess = 'CHAINE' 'Flux traversant ' DFluxTIT ;
  847.   TabDess   = 'TABLE';
  848.   TabDess.1 = 'MARQ CARR' ;
  849.   TabDess.2 = 'MARQ TRIA' ;
  850.   TabDess.3 = 'MARQ PLUS' ;
  851.   TabDess.'TITRE'     = 'TABLE';
  852.   TabDess.'TITRE' . 1 = 'Debit_Convectif';
  853.   TabDess.'TITRE' . 2 = 'Debit_Diffusif';
  854.   TabDess.'TITRE' . 3 = 'Debit_Total';
  855.   EvTot = EvFlC ET EvFlD ET EvFl;
  856.   'DESSIN' EvTot 'MIMA' TabDess 'TITR' TitDess 'TITX' 't (an)' 'TITY' 
  857.                  'mol/l/an' 'LEGE' 'AXES' 'XBOR' TIni TFin;
  858. FINSI;
  859.  
  860. *                                           ===========================
  861. *                                           | TEST DE BON DEROULEMENT |
  862. *                                           ===========================
  863. * Le test se fait en comparant le flux total traversant la ligne avec
  864. * une liste obtenue par un calcul préalable.
  865. liflsav = 'PROG' 0.56061 1.0216 1.2444 1.3013 1.2365 1.0595 0.93058
  866.                  0.64840 0.41501 0.29818;
  867.  
  868. err1 = 'MAXIMUM' ('ABS' (lifl - liflsav)) / ('MAXIMUM' liflsav);
  869. 'MESSAGE' 'difference avec sol de ref' err1;
  870.  
  871. si (err1 > CRIT1);
  872.   'ERREUR' 5;
  873. sinon;
  874.   'ERREUR' 0;
  875. finsi;
  876.  
  877. 'LISTE' crit1;
  878.  
  879. Fin;
  880.  
  881.  
  882.  
  883.  
  884.  
  885.  
  886.  
  887.  
  888.  
  889.  
  890.  
  891.  
  892.  
  893.  
  894.  

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