Télécharger la plaquette Cast3M


fr en

Télécharger transport5.dgibi

Retour à la liste

Numérotation des lignes : 

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

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