Télécharger chimsour1d.dgibi

Retour à la liste

Numérotation des lignes :

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

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