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