Télécharger la plaquette Cast3M


fr en

Télécharger chimsour1d.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  chimsour1d.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. GRAPH = faux ;
  5. * repertoire des fichiers "divers"
  6. DIVERS = VENV 'CASTEM_DIVERS';
  7. *
  8. ******************* CAS TEST : chimsour1d.dgibi **********************
  9.  
  10. * --------------------------------------------------------------------
  11. *
  12. *   TRANSPORT GEOCHIMIE AVEC SOURCE ( CAS 1D)
  13. *    la partie post-traitement avec graphique  contient des exemples
  14. *       d'utilisation des procedures TRACHIS TRACHIT DESTRA
  15. *    avec des recherches d'identificateurs par NOESPCHI et NOCOMCHI
  16. *
  17. * --------------------------------------------------------------------
  18.  
  19. * emplacement de COMPOM
  20. emp1 = 'CHAINE' DIVERS '/COMPOM' ;
  21. *emp1 = 'MOT' 'COMPOM' ;
  22. *emp1 = 'MOT' '/export/home/castem2001/DGIBI/COMPOM' ;
  23. 'OPTION' 'ECHO' 0 ;
  24. *
  25. *        Génération du maillage
  26. *        ----------------------
  27. *
  28. 'OPTION' 'DIME' 2 'ELEM' 'QUA4' 'TRACER' 'PSC' ;
  29.  
  30. * Axe
  31. ORIG   = 0.0 0.0 ;
  32. R10    = 1. 0.0;
  33. XAXE   = 'DROIT' 1 ORIG R10;
  34. 'ELIMINATION' XAXE 0.01;
  35.  
  36. * Points
  37. Z90    = 0.0   90.0;
  38. Z5     = 0.0    5.0;
  39. Z10    = 0.0   10.0;
  40. Z40    = 0.0   40.0 ;
  41. Z390   = 0.0  390.0;
  42. ZNAD   = 0.0 -100.0;
  43. ZM40   = 0.0  -40.0 ;
  44. ZBOT   = 0.0  -10.0;
  45. ZB2    = 0.0   -5.0;
  46. ZT2    = 0.0   5.0;
  47. ZTOP   = 0.0  10.0;
  48. ZZEN   = 0.0 400.0;
  49. *
  50. XAXE2  = 'PLUS' XAXE ZNAD;
  51. *
  52. ZONE1  = 'TRANSLATION' XAXE2 6 ('MOIN' ZM40 ZNAD) ;
  53. XAXM40 = 'INVERSE' ('COTE' 3 ZONE1) ;
  54. YN     = 'INVERSE' ('COTE' 4 ZONE1) ;
  55. LIGM40 = 'ELEM' ZONE1 'APPUYE' 'LARGEMENT'  XAXM40 ;
  56. ZONE2  = 'TRANSLATION' XAXM40 6 ('MOIN' ZBOT ZM40) ;
  57. XAXBOT = 'INVERSE' ('COTE' 3 ZONE2) ;
  58. YK     = 'INVERSE' ('COTE' 4 ZONE2)  ;
  59. LIGBOT = 'ELEM' ZONE2 'APPUYE' 'LARGEMENT'  XAXBOT ;
  60. ZONE3  = 'TRANSLATION' XAXBOT 8 ('MOIN' ZTOP ZBOT) ;
  61. XAXTOP = 'INVERSE' (COTE 3 ZONE3) ;
  62. LIGTOP = 'ELEM' ZONE3 'APPUYE' 'LARGEMENT'  XAXTOP ;
  63. YB     = 'INVERSE' (COTE 4 ZONE3)  ;
  64. ZONE4  = 'TRANSLATION' XAXTOP 6 ('MOIN' Z40 ZTOP) ;
  65. XAXP40 = 'INVERSE' (COTE 3 ZONE4) ;
  66. YS     = 'INVERSE' (COTE 4 ZONE4)  ;
  67. ZONE5  = 'TRANSLATION' XAXP40 36 ('MOIN' ZZEN Z40) ;
  68. LIGP40 = 'ELEM' ZONE5 'APPUYE' 'LARGEMENT'  XAXP40 ;
  69. XAXE1  = 'INVERSE' (COTE 3 ZONE5) ;
  70. YZ     = 'INVERSE' (COTE 4 ZONE5)  ;
  71. LIGAX1 = 'ELEM' ZONE5 'APPUYE' 'LARGEMENT'  XAXE1 ;
  72. *
  73. YAXE   = YN 'ET' YK 'ET' YB 'ET' YS 'ET' YZ;
  74. X80    = 80.0 0.0;
  75. *
  76. MT     = ZONE1 'ET' ZONE2 'ET' ZONE3 'ET' ZONE4 'ET' ZONE5  ;
  77. *
  78. 'ELIMINATION' (MT 'ET' YAXE) 0.01;
  79. 'TASSER' MT;
  80. SDOM   = 'ELEM' MT 'CONTENANT' ORIG ;
  81. SDOM   = 'COULEUR' SDOM 'ROUGE' ;
  82. AXELEM = 'ELEM' MT 'APPUYE' 'LARGEMENT' YAXE ;
  83. *
  84. * Calcul CHI1
  85. * -----------
  86. *
  87. TABDON              = 'TABLE' ;
  88. TABDON.IDEN         = 'LECT' 50 5 112 ;
  89. TABDON.CLIM         = 'TABLE' ;
  90. TABDON.CLIM . TYP3  = 'LECT' 21442  ;
  91. TABDON.CLIM . COMP3 = 'LECT' 112 ;
  92. TABDON.CLIM . TYP6  = 'LECT' 21440 21441 12721 12722 ;
  93.  
  94. TB1 = 'CHI1' TABDON 'COMP' emp1 'LOGK' emp1 ;
  95. *
  96. *- Création des maillages hybrides
  97. *
  98. MFTOT  = 'CHANER' MT     'QUAF'  ;
  99. MFSOUR = 'CHANER' SDOM   'QUAF'  ;
  100. MFBAS  = 'CHANER' XAXE2  'QUAF'  ;
  101. MFAXE  = 'CHANER' AXELEM 'QUAF'  ;
  102. MFP40  = 'CHANER' LIGP40 'QUAF'  ;
  103. MFTOP  = 'CHANER' LIGTOP 'QUAF'  ;
  104. 'ELIMINATION' 0.001
  105. (MFTOT 'ET' MFSOUR 'ET' MFBAS 'ET' MFAXE 'ET' MFP40 'ET' MFTOP) ;
  106. *
  107. *- Modèle
  108. *
  109. MODHYB = 'MODELE' MFTOT 'DARCY' ANISOTROPE  ;
  110. MMSOUR = 'MODELE' MFSOUR 'DARCY' ANISOTROPE  ;
  111. MMBAS  = 'MODELE' MFBAS 'DARCY' ANISOTROPE  ;
  112. MMAXE  = 'MODELE' MFAXE 'DARCY' ANISOTROPE  ;
  113. MMP40  = 'MODELE' MFP40 'DARCY' ANISOTROPE  ;
  114. MMTOP  = 'MODELE' MFTOP 'DARCY' ANISOTROPE  ;
  115. CHYB1  = 'DOMA' MODHYB 'SURFACE' ;
  116. CHYB2  = 'DOMA' MODHYB 'NORMALE' ;
  117. XVOLU  = 'DOMA' MODHYB 'VOLUME' ;
  118. CETOT  = 'DOMA' MODHYB 'CENTRE' ;
  119. MATOT  = 'DOMA' MODHYB 'MAILLAGE ';
  120. FATOT  = 'DOMA' MODHYB 'FACE' ;
  121. *
  122. *- On génère la ligne SEGAXE pour le post-traitement
  123. *
  124. CEAXE  = 'DOMA' MMAXE 'CENTRE' ;
  125. NBCC   = 'NBNO' CEAXE ;
  126. PI1    = 'POINT' CEAXE 1 ;
  127. PI0    = PI1 ;
  128. I      = 2 ;
  129. 'SI' (NBCC > 1) ;
  130.   PI2    = 'POINT' CEAXE I ;
  131.   SEGAXE = 'QUELCONQUE' 'SEG2' PI1 PI2 ;
  132.   SI (NBCC > 2) ;
  133.     NBCC2 = NBCC - 2 ;
  134.     'REPETER' BLOC6 nbcc2 ;
  135.       PI1    = PI2 ;
  136.       I      = I + 1 ;
  137.       PI2    = 'POINT' CEAXE I ;
  138.       LILI   = 'QUELCONQUE' 'SEG2' PI1 PI2 ;
  139.       SEGAXE = SEGAXE 'ET' LILI ;
  140.     'FIN' BLOC6 ;
  141.   'FINSI' ;
  142. 'FINSI' ;
  143. *
  144. * Données physiques
  145. * -----------------
  146. *
  147. *- on entre la vitesse on en deduit le flux aux faces
  148. *
  149. V     = 'MANU' CHPO  FATOT  'NATURE' 'DISCRET' 2 'UX' 0. 'UY' 0.1;
  150. VNCH  = 'VECTEUR' V 1. 'UX' 'UY' 'ROUGE' ;
  151. MOT1  = 'MOTS' 'UX' 'UY' ;
  152. VAVN  = 'PSCAL' V CHYB2 MOT1 MOT1 ;
  153. VAVN  = 'NOMC' 'SCAL' VAVN ;
  154. QFACE = VAVN * CHYB1 ;
  155. QFACE = 'NOMC' 'FLUX' QFACE ;
  156. *
  157. *- matériau pour le transport
  158. *
  159. DIR1  = 1. 0. ;
  160. PORO  = 0.5 ;
  161. D11C  = 1.D0 ;
  162. D22C  = 5.D0 ;
  163. D11   = 'MANU' CHML MATOT  SCAL D11C ;
  164. D22   = 'MANU' CHML MATOT  SCAL (D22C* poro)  ;
  165. D21   = 'MANU' CHML MATOT  SCAL 1.D-15 ;
  166. MAT1  = 'MATERIAU' MODHYB 'DIRECTION' DIR1 K11 D11 K21 D21 K22 D22 ;
  167. *
  168. *- pas de temps
  169. *
  170. *NBPAST nombre de pas de temps  (200 pour un calcul significatif)
  171. DELTAT= 1.25 ;
  172. NBPAST= 10 ;
  173. TETA  = 1. ;
  174. TFINAL                 = NBPAST * DELTAT ;
  175. *
  176. * Initialisation du système
  177. * -------------------------
  178. *
  179. * CTOTC  concentrations au centre
  180. * CTOT   concentrations aux faces
  181. *
  182. CTOT   = 'MANU' 'CHPO' FATOT 3 X050 2.D-9 X005 1.D-9 X112 0.D0 ;
  183. TT1    = 'EXTR' CTOT 'COMP' ;
  184. NOCOMP = 'EXTR' TT1 1 ;
  185. CTOTC  = 'MANU' 'CHPO' CETOT 1 NOCOMP 0. ;
  186. NBCOMP = 'DIME' TT1 ;
  187. MASHYB = 'MHYB' MODHYB MAT1 ;
  188. 'REPETER' BOUC1 NBCOMP ;
  189.   NOCOMP = 'EXTR' TT1 &BOUC1 ;
  190.   TPR    = 'EXCO' CTOT NOCOMP 'TH' ;
  191.   TPC    = 'HYBP' MODHYB MASHYB TPR ;
  192.   TPC    = 'NOMC' NOCOMP TPC ;
  193.   CTOTC  = CTOTC + TPC ;
  194. 'FIN' BOUC1 ;
  195. CLOGC = 'MANU' 'CHPO' CETOT 3 X050 -7. X005 -9. X112 -7. ;
  196. FLOGC = 'MANU' 'CHPO' FATOT 3 X050 -7. X005 -9. X112 -7. ;
  197. *
  198. * Source
  199. *
  200. QELEM0 = 1.D-3;
  201. QELEM  = QELEM0 / DELTAT  ;
  202. SOURC0 = 'MANU' 'CHPO' CETOT 3  X050 0.  X005 0. X112 0. ;
  203. SOURC1 = 'MANU' 'CHPO' ('DOMA' MMSOUR 'CENTRE') 3
  204.                 X050 (-1.* QELEM) X005 QELEM X112 0. ;
  205. SOURC1 = SOURC1 + SOURC0 ;
  206. SOURC  = 'CHARGEMENT' SOURC1
  207.            ('EVOL' 'MANU' ('PROG' 0. 'PAS' DELTAT (200. * DELTAT))
  208.                           (PROG 1. 1. 199. * 0.) )  ;
  209. *
  210. * On initialise les concentrations des aqueux aux faces
  211. * par un premier calcul CHI2
  212. *
  213. TBPAR2            = 'TABLE' ;
  214. TBPARM            = 'TABLE' ;
  215. TBPARM.'SOUSTYPE' = 'DONNEES_CHIMIQUES' ;
  216. TBPARM.TOT        = CTOT ;
  217. TBPARM.LOGC       = FLOGC ;
  218. TBPAR2.ITMAX      = 25 ;
  219. TBPAR2.EPS        = 1.D-4 ;
  220. TBPAR2.NFI        = 2 ;
  221. *
  222. TB4 = 'CHI2' TB1 TBPAR2 TBPARM ;
  223. *
  224. TAQU = TB4.AQUE ;
  225. *
  226. * conditions aux limites
  227. *
  228. LIMBAS = 'REDU' TAQU ('DOMA' MMBAS 'CENTRE') ;
  229. BBBAS  = 'BLOQ' ('DOMA' MMBAS 'CENTRE') 'TH' ;
  230. *
  231. * Table de données de TRANSPORT
  232. * -----------------------------
  233. *
  234. TABTRAN                = 'TABLE' ;
  235. TABTRAN.SOUSTYPE       = 'MOT' 'GEOCHIMIE' ;
  236. TABTRAN.'MODELE'       = MODHYB ;
  237. TABTRAN.'DIFFUSION'    = MAT1 ;
  238. TABTRAN.'POROSITE'     = 'MANU' 'CHPO' CETOT 1 'CK' PORO ;
  239. TABTRAN.'CONVECTION'   = QFACE ;
  240. TABTRAN.'CHIMI1'       = TB1 ;
  241. TABTRAN.'TAQU'         = 'TABLE' ;
  242. TABTRAN.'TAQU'. 0      = TAQU ;
  243. TABTRAN.ITMAX          = 25;
  244. TABTRAN.EPS            = 1.D-4 ;
  245. TABTRAN.NFI            = 1 ;
  246. TABTRAN.LOGC           = 'TABLE' ;
  247. TABTRAN.LOGC. 0        = CLOGC ;
  248. TABTRAN.TOT            = 'TABLE' ;
  249. TABTRAN.TOT. 0         = CTOTC ;
  250. TABTRAN.PRECISION      = 5.E-3 ;
  251. TABTRAN.SORTIE         = 'MOTS' 'FION' 'SOLU'  ;
  252. TABTRAN.'BLOCAGE'      = BBBAS  ;
  253. TABTRAN.'TRACE_IMPOSE' = 'CHARGEMENT' LIMBAS
  254.                 ('EVOL' 'MANU' ('PROG' 0. 450.) ('PROG' 1. 1.))  ;
  255. TABTRAN.'SOURCE'       =  SOURC  ;
  256. TABTRAN.'PAS_DE_TEMPS' = DELTAT ;
  257. TABTRAN.'TEMPS_FINAL'  = TFINAL ;
  258. TABTRAN.'THETA'        = 1.D0 ;
  259. *
  260. *  Calcul couplé transport/géochimie
  261. *  ---------------------------------
  262. *
  263. CHITRNSP TABTRAN ;
  264. *
  265. * Controle des résultats
  266. * ----------------------
  267. *
  268. * Calcul théorique  1D   (pour le NA+) 
  269. *
  270. VITCENT = 'HVIT' MODHYB QFACE  ;
  271. UY      = 'EXCO' VY VITCENT SCAL ;
  272. XX YY   = 'COOR' CETOT ;
  273. YY      = 'NOMC' YY SCAL ;
  274. PETITT  = TFINAL ;
  275. DENO1   = (4. * PI * PETITT * D22C) ** 0.5 ;
  276. DENO2   = (4. * PETITT * D22C) ** 0.5 ;
  277. CXYT1   = QELEM0 / PORO / DENO1 ;
  278. CXYT21  = ((YY -((UY / PORO) * PETITT)) / DENO2) ** 2 ;
  279. CXYT2   = CXYT1 * ('EXP' (-1. * CXYT21)) ;
  280. SC1     = CXYT2 * XVOLU ;
  281. SCXYT   = 'RESULT' SC1 ;
  282. tttt    = 'CHAINE' 'Concentration en NA+ le long de l axe au temps '
  283.                    petitt ;  
  284. 'TITRE' TTTT ;
  285. EV100   = 'EVOL' CHPOI CXYT2 'SCAL' SEGAXE ;
  286. *
  287. *- Calcul de l'intégrale de la concentration de NA+ et comparaison
  288. * à la valeur théorique  (on compare l'intégrale et la valeur maximale)
  289. *
  290. ECXN    = 'EXCO' W002 TABTRAN.'SOLU'. NBPAST 'SCAL' ;
  291. SECX1   = ECXN * XVOLU ;
  292. SECXN   = 'RESULT' SECX1 ;
  293. SS1     = 'EXTR' SCXYT 'SCAL' ('POINT' ('EXTR' SCXYT 'MAIL') 1)  ;
  294. SS2     = 'EXTR' SECXN 'SCAL' ('POINT' ('EXTR' SECXN 'MAIL') 1)  ;
  295. DIFS    = ('ABS' (SS1 - SS2) ) / SS1;
  296. 'SI' (DIFS  < 1.D-2)  ;
  297.    'ERREUR' 0 ;
  298. 'SINON' ;
  299.    'ERREUR' 5 ;
  300. 'FINSI' ;
  301.  
  302. MAECXN  = 'MAXIMUM' ECXN ;
  303. MACXYT2 = 'MAXIMUM' CXYT2 ;
  304. DIFFNA  = 'ABS'  (ECXN - CXYT2) ;
  305. MADIFF  = 'MAXIMUM' DIFFNA ;
  306. ERRDIF  = MADIFF / MACXYT2 ;
  307. 'SI'  (ERRDIF > 1.D-1)  ;
  308.    'ERREUR' 5 ;
  309. 'FINSI' ;
  310. *
  311. * Post-traitement graphique
  312. * -------------------------
  313. *
  314. 'SI' GRAPH ;
  315.   EV51    = 'EVOL' 'VERT' CHPOI ECXN 'SCAL' SEGAXE ;
  316.   'DESSIN' (ev51 'ET' ev100) ;
  317.   MOTIA   = 'CHAINE' 'Concentration le long de l axe'
  318.                      ' pour différents temps ';
  319.   LLTMP   = 'LECT' 0 'PAS' 10 NBPAST  ;
  320.   LISMD   = 'MOTS' 'W002' ;
  321.   NE1 NE2 = NOESPCHI TB1  W002 ;
  322.   NO1 NO2 NO3 = NOCOMCHI TB1 'NUMCOMP' NE2 ;
  323.   TBDES   = TRACHIS TABTRAN 'SOLU' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
  324.   DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
  325.  
  326.   LISMD   = 'MOTS' 'X050' ;
  327.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NOMINT' X050 ;
  328.   TBDES   = TRACHIS TABTRAN 'AQUE' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
  329.   DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
  330.  
  331.   LISMD   = 'MOTS' 'X005' ;
  332.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NOMINT' X005 ;
  333.   LISMD   = 'MOTS' 'X112' ;
  334.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NOMINT' X112 ;
  335.   LISMD   = 'MOTS' 'W001' ;
  336.   NE1 NE2 =  NOESPCHI TB1  W001 ;
  337.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NUMCOMP' NE2 ;
  338.   TBDES   = TRACHIS TABTRAN 'SOLU' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
  339.   DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
  340.  
  341.   LISMD   = 'MOTS' 'W002' ;
  342.   NE1 NE2 =  NOESPCHI TB1  W002 ;
  343.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NUMCOMP' NE2 ;
  344.   TBDES   = TRACHIS TABTRAN 'SOLU' LLTMP LISMD ('MOTS' NO1) SEGAXE ;
  345.   DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' MOTIA ;
  346.  
  347.   TBN     = 'TABLE' ;
  348.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NUMCOMP' 50 ;
  349.   TBN.1   = NO1;
  350.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NUMCOMP' 5  ;
  351.   TBN.2   = NO1;
  352.   NO1 NO2 NO3 =  NOCOMCHI TB1 'NUMCOMP' 112 ;
  353.   TBN.3   = NO1;
  354.   MO1 NO1 = NOESPCHI TB1 50 ;
  355.   MO2 NO2 = NOESPCHI TB1 5 ;
  356.   MO3 NO3 = NOESPCHI TB1 112 ;
  357.   LISMD   = 'MOTS' MO1 MO2 MO3 ;
  358.   CETOP   = 'DOMA' MMTOP 'CENTRE' ;
  359.   XXTOP YYTOP= 'COOR' CETOP ;
  360.   YCTOP   = 'EXTR' YYTOP 'SCAL' ('POINT' ('EXTR' YYTOP 'MAIL') 1)  ;
  361.   TITOP   = CHAIN 'Evolution en fonction du temps a '
  362.                   'FORMAT' '(F6.3)' YCTOP 'metres' ;
  363.   TBDES   = TRACHIT TABTRAN 'SOLU' LISMD TBN CETOP ;
  364.   DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' TITOP ;
  365.  
  366.   LISMD   = 'MOTS' 'X050' 'X005' 'X112' ;
  367.   TBDES   = TRACHIT TABTRAN 'AQUE' LISMD TBN CETOP ;
  368.   DESTRA TBDES 'MIMA' 'LOGO' 'DATE' 'TITRE' TITOP ;
  369. 'FINSI' ;
  370.  
  371. 'FIN' ;
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  

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