Télécharger transgen.procedur

Retour à la liste

Numérotation des lignes :

  1. * TRANSGEN PROCEDUR JC220346 12/09/12 21:15:09 7501
  2. *TRANSGEN PROCEDUR GBM 04/01/03
  3. 'DEBP' TRANSGEN TRANSI*'TABLE';
  4.  
  5.  
  6. * SUPPRIMER TOUS LES 'H', stocker composante de CHRG 'ET'
  7. * REMPLACER LES 'H' par cela. A faire.
  8.  
  9.  
  10.  
  11. *---------------------------------------------------------------------
  12. * Résolution de l'équation de darcy pour un problème d'écoulement ou
  13. * de transport par une méthode d'éléments
  14. * finis mixtes hybrides ou VF. Les inconnues du problème sont
  15. * - en EFMH, la concentration (H), la trace de
  16. * concentration (TH) et le débit diffusif (FLUX).
  17. * - en VF, la concentration (H)
  18. *---------------------------------------------------------------------
  19. *
  20. *------------------------------
  21. * Phrase d'appel (en GIBIANE) :
  22. *------------------------------
  23. *
  24. * TRANSGEN TABLE ;
  25. *
  26. *----------------------------------
  27. * Opérandes (à mettre dans TABLE) :
  28. *----------------------------------
  29. *
  30. * ___________________________________________________________________
  31. * | |
  32. * | Indice Contenu |
  33. * | |
  34. * -------------------------------------------------------------------
  35. * | |
  36. * |------------------------------------------------ |
  37. * |Données physiques, géométriques et materielles : |
  38. * |------------------------------------------------ |
  39. * | |
  40. * |'MODELE' Objet modèle (MMODEL créé par MODE) |
  41. * | |
  42. * |'CARACTERISTIQUES' Données physiques et materielles : |
  43. * | diffusivité effective (CHAMELEM créé par MATE) |
  44. * | |
  45. * |'POROSITE' Valeur de la porosité (Type Champoint, Comp |
  46. * | 'SCAL', ou FLOTTANT) - Défaut 1. |
  47. * | |
  48. * |'DECROISSANCE' Valeur du terme de décroissance (Type FLOTTANT) |
  49. * | Tel que dC/dt = - Lambda * C - Défaut 0. |
  50. * | |
  51. * |'COEF_RETARD' Coefficient de retard linéaire dans le cas simple, |
  52. * | ou Pente à l'origine de la fonction F(C) dans le |
  53. * | cas d'isotherme non linéaire de Langmuir |
  54. * | ou Coefficient K de l'isotherme de Freundlich |
  55. * | (Type CHPO Centre 'SCAL', ou FLOTTANT) |
  56. * | |
  57. * |'LANGMUIR' Quantité maximale 'Fsat' adsorbée sur le solide |
  58. * | rapportée à l'unité de volume du fluide et exprimée|
  59. * | dans la meme unité que le soluté. |
  60. * | (Type CHPO Centre 'SCAL', ou FLOTTANT). |
  61. * | F = (R-1) C / [1 + ((R-1) C / Fsat)] |
  62. * | Si cet indice et le suivant sont absents, |
  63. * | l'équilibre d'adsorption est linéaire. Cet indice a|
  64. * | priorité sur l'indice FREUNDLICH. |
  65. * | |
  66. * |'FREUNDLICH' Exposant de la loi de Freundlich F = K (C ^ 1/n) |
  67. * | (Type FLOTTANT). |
  68. * | Dans ce cas (et si l'indice LANGMUIR n'existe pas),|
  69. * | l'indice 'COEF_RETARD' contient le coefficient |
  70. * | K ramené à une unité de volume de fluide. |
  71. * | - Non disponible pour l'instant - |
  72. * | |
  73. * |'LIMITE_SOLUBILITE' Limite de solubilité (Type chpoin, Comp 'H') |
  74. * | |
  75. * |'COEF_DISSOLUTION' Coef. de dissolution (Type CHPO Centre, Comp |
  76. * | 'SCAL'). Tel que dC/dt = Coef * (Csat - C) - Par |
  77. * | défaut, la dissolution est instantanée |
  78. * | |
  79. * |'CONVECTION' Débit intégré de la vitesse convective à travers |
  80. * | chaque face des éléments (Type CHPO Face, comp. |
  81. * | 'FLUX') |
  82. * | |
  83. * |---------------------- |
  84. * |Conditions initiales : |
  85. * |---------------------- |
  86. * | |
  87. * |'TEMPS' TABLE contenant à l'indice 0 la valeur du temps |
  88. * | initial (FLOTTANT) |
  89. * | |
  90. * | |
  91. * |'CONCENTRATION' TABLE contenant à l'indice 0 la concentration |
  92. * | (quantité d'élément par unité de volume d'eau) |
  93. * | (Type CHPO Centre, Comp 'H') |
  94. * | |
  95. * |'FLUX' TABLE contenant à l'indice 0 le flux diffusif |
  96. * | initial intégré sur chaque face (Type CHPO Face, |
  97. * | comp. 'FLUX') |
  98. * | |
  99. * |'PRECIPITE' TABLE contenant à l'indice 0 la quantité initiale |
  100. * | de précipité par unité de volume de milieu solide |
  101. * | (Type CHPO Centre, Comp 'H') |
  102. * | |
  103. * |'DISSOLUTION' TABLE contenant à l'indice 0 la quantité initiale |
  104. * | pour estimer la dissolution au premier pas de temps|
  105. * | (Type CHPO, Comp 'H'), voir plus loin. |
  106. * | |
  107. * |-------------------------------------- |
  108. * |Conditions aux limites / chargements : |
  109. * |-------------------------------------- |
  110. * | |
  111. * |'BLOCAGE' Contient les matrices de blocage (RIGIDITE) |
  112. * | |
  113. * |'TRACE_IMPOSE' Valeurs des traces imposées (charge ou concentra- |
  114. * | -tion) (CHARGEMENT 'TH' - Obligatoire si BLOCAGE) |
  115. * | |
  116. * |'FLUX_IMPOSE' Valeurs des flux imposés intégrés par face |
  117. * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- |
  118. * | défaut 0.) |
  119. * | |
  120. * |'FLUXTOT_IMP' Valeurs des flux totaux imposés intégrés par face |
  121. * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- |
  122. * | défaut 0.). Il s'agit du flux diffusif + convectif |
  123. * | |
  124. * |'MIXTES' Table : indice C contient les valeurs des flux |
  125. * | mixtes imposés intégrés par face |
  126. * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- |
  127. * | défaut 0.) il est égal à A * flux diffusif + |
  128. * | B * Concentration |
  129. * | Indice A et B sont les réels A et B |
  130. * | |
  131. * |'SOURCE' Valeurs du terme source par maille et par unité de |
  132. * | temps (ex : puits, filiation) |
  133. * | Les valeurs à l'indice i sont les valeurs entre |
  134. * | les temps i-1 et i. |
  135. * | (CHARGEMENT de CHPO Centre, comp 'SOUR'- défaut 0.)|
  136. * | |
  137. * | |
  138. * |'DISSOLUTION_IMPOSEE' Valeurs des dissolutions imposées par unité|
  139. * | de temps et par maille. (Type CHARGEMENT de CHPO, |
  140. * | Comp 'H'). Les valeurs à l'indice i sont les |
  141. * | valeurs moyennes de dissolution par unité de temps |
  142. * | entre les temps i-1 et i. |
  143. * | Priorité de la dissolution imposée sur les |
  144. * | cinétiques. |
  145. * | |
  146. * |-------------------- |
  147. * |Données numériques : |
  148. * |-------------------- |
  149. * | |
  150. * | |
  151. * |'TEMPS_CALCULES' Valeur des temps calculés (LISTREEL) |
  152. * | Contient obligatoirement le temps final. |
  153. * | |
  154. * |'TEMPS_SAUVES' Valeur des temps sauvegardés (LISTREEL - défaut : |
  155. * | on sauve tous les pas de temps) |
  156. * | |
  157. * | |
  158. * | |
  159. * |'THETA_DIFF' Coefficient de relaxation compris entre 0. et 1. |
  160. * | (theta-méthode diffusion) ('FLOTTANT' - défaut 1.) |
  161. * | |
  162. * |'THETA_CONV' Idem pour la convection |
  163. * | ('FLOTTANT', Défaut = THETA_DIFF) |
  164. * |'THETA_DEC' Idem mais pour la décroissance |
  165. * | ('FLOTTANT' - défaut 1/2) |
  166. * | |
  167. * |'THETA_DISS' Idem mais pour la dissolution |
  168. * | ('FLOTTANT' - défaut 1.) |
  169. * | |
  170. * |'PENALISATION' Coefficient de pénalisation pour la prise en |
  171. * | compte de la limite de solubilité. La présence de |
  172. * | cet indice ou du suivant indique quel schéma a été |
  173. * | choisi. |
  174. * | (Type 'FLOTTANT') - Valeur conseillée 1.D7 |
  175. * | |
  176. * |'EPSI_LIM' Précision relative d'arrêt pour le shéma limite de |
  177. * | solubilité prédicteur-correcteur itératif |
  178. * | (Type FLOTTANT) - Valeur conseillée 5.D-3 |
  179. * | |
  180. * |'ITMAX_LIM' Nombre maxi d'itérations correspondant aux modules |
  181. * | de dissolution avant d'abandonner |
  182. * | (Type 'ENTIER') - Défaut 50 |
  183. * | |
  184. * |'EPSI_RET' Précision relative d'arrêt pour la résolution |
  185. * | itérative (Picard) de l'adsorption non linéaire |
  186. * | (Type FLOTTANT) - Défaut 1.D-4 |
  187. * | |
  188. * |'EPSI_COR' Petit saut de concentration pour calculer le coef. |
  189. * | de retard par la méthode de la corde lorsque le |
  190. * | retard est non-linéaire. |
  191. * | (Type FLOTTANT) - Défaut 1.D-4 |
  192. * | |
  193. * |'ITMAX_RET' Nombre maxi d'itérations correspondant au retard |
  194. * | non linéaire avant d'abandonner. |
  195. * | (Type 'ENTIER') - Défaut 20 |
  196. * |_________________________________________________________________|
  197. *
  198. *
  199. *
  200. *---------------------------------
  201. * Résultats (stockés dans TABLE) :
  202. *---------------------------------
  203. *
  204. * ___________________________________________________________________
  205. * | |
  206. * | Indice Contenu |
  207. * | |
  208. * -------------------------------------------------------------------
  209. * | |
  210. * | |
  211. * |'TEMPS' TABLE contenant les temps sauvegardés (FLOTTANT) |
  212. * | |
  213. * |'CONCENTRATION' TABLE contenant les concentrations |
  214. * | (Type CHPO Centre, Comp 'H') |
  215. * | |
  216. * |'FLUX' TABLE contenant les débits diffusifs intégrés |
  217. * | par face (Type CHPO Face, comp. 'FLUX') |
  218. * | |
  219. * |'PRECIPITE' TABLE contenant la quantité de précipité par maille|
  220. * | (Type CHPO Centre, Comp 'H') |
  221. * | |
  222. * |'DISSOLUTION' TABLE contenant la quantité de précipité dissoute |
  223. * | entre deux pas de temps par unité de volume et par |
  224. * | unité de temps. La valeur stockée à l'indice i, |
  225. * | est valable entre les temps i-1 et i |
  226. * | (Type CHPO, Comp 'H'). |
  227. * | ATTENTION, les valeurs de cette table résultat |
  228. * | n'ont aucun sens lorsque les temps sauvegardés ne |
  229. * | sont pas les memes que les temps calculés. Toute |
  230. * | tentative d'exploitation donnera alors des |
  231. * | résultats incohérents (erreurs de bilan) |
  232. * | |
  233. * |'RETARD' Si cet indice a été préalablement défini comme une |
  234. * | TABLE, alors il contient les valeurs du coefficient|
  235. * | de retard (Type 'CHPO' centre, Comp 'SCAL'). Sinon,|
  236. * | les valeurs du coefficient de retard ne sont pas |
  237. * | sauvegardées. |
  238. * |_________________________________________________________________|
  239. *
  240. *
  241. * ___________________________________________________________________
  242. * | |
  243. * | Les tables résultats sont indicées par des entiers variant de 0 |
  244. * | à N . |
  245. * | A l'indice 0 on stocke les valeurs initiales, aux indices |
  246. * | suivants les champs correspondant au temps de sortie TEMPS.I . |
  247. * | Les champs servant en cas de reprise sont ceux correpondant au |
  248. * | dernier indice. |
  249. * |_________________________________________________________________|
  250. *
  251. *
  252. * TYPE DE RESOLUTION
  253.  
  254. NOMINC = 'MOT' 'CONCENTRATION' ;
  255. NOMTETA = 'MOT' 'THETA_DIFF' ;
  256. NOMDDT = 'MOT' 'POROSITE' ;
  257.  
  258. *
  259. *--------------------------------------------------------------------*
  260. * RECUPERATION DES DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES *
  261. *--------------------------------------------------------------------*
  262. *
  263. * MODELE
  264. 'SI' ( 'EXISTE' TRANSI 'MODELE' ) ;
  265. MODDARCY = TRANSI . 'MODELE' ;
  266. * METTRE LE PRECOND à 1 pour ne pas recalculer les préconditionnements
  267. * sur le maillage
  268. ('DOMA' MODDARCY TABLE) . 'PRECONDI' = 1 ;
  269. 'SINON' ;
  270. 'ERREUR' 'Il manque le modele.' ;
  271. 'QUITTER' TRANSGEN ;
  272. 'FINSI' ;
  273.  
  274. * VOLUME
  275. VOLU1 = 'DOMA' MODDARCY 'VOLUME' ;
  276. * ORIENTATION
  277. MCHYB = 'DOMA' MODDARCY 'ORIENTAT' ;
  278. *
  279. MCENT = 'DOMA' MODDARCY 'CENTRE' ;
  280. *
  281. MFACE = 'DOMA' MODDARCY 'FACE' ;
  282. *
  283. MAILL = 'DOMA' MODDARCY 'MAILLAGE' ;
  284. * FILIATION
  285. 'SI' ('NON' ('EXISTE' TRANSI 'PERE')) ;
  286. TRANSI . 'PERE' = FAUX ;
  287. 'FINSI' ;
  288. * CARACTERISTIQUES
  289. 'SI' ( 'EXISTE' TRANSI 'CARACTERISTIQUES' ) ;
  290. 'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHPOINT') ;
  291. MAT1 = TRANSI . 'CARACTERISTIQUES' ;
  292. 'FINSI' ;
  293. * Si la diffusion est un chargement, le traitement à lieu dans le
  294. * corps du programme
  295. 'SINON' ;
  296. 'ERREUR'
  297. 'Indice CARACTERISTIQUES absent de la table de données.' ;
  298. 'QUITTER' TRANSGEN ;
  299. 'FINSI' ;
  300. * EMMAGASINEMENT OU POROSITE
  301. 'SI' ( 'EXISTE' TRANSI NOMDDT ) ;
  302. POROS = TRANSI . NOMDDT ;
  303. 'SINON' ;
  304. POROS = 1.D0 ;
  305. 'FINSI' ;
  306. 'SI' ('EGA' ('TYPE' POROS) 'FLOTTANT') ;
  307. POROS = 'MANU' 'CHPO' MCENT 'SCAL' POROS ;
  308. POROS = 'CHAN' 'ATTRIBUT' POROS 'NATURE' 'DISCRET' ;
  309. * VERRUE
  310. * POROS = 'MANU' 'CHPO' MCENT 'SCAL' 1.D0 NATURE DISCRET ;
  311. *
  312. 'FINSI' ;
  313. 'SI' ('NEG' ('TYPE' POROS) 'CHPOINT ') ;
  314. 'ERREUR' ('CHAINE' 'L indice ' NOMDDT ' doit etre de type '
  315. 'FLOTTANT ou CHAMPOINT') ;
  316. 'SINON' ;
  317. POROS = 'NOMC' 'SCAL' POROS ;
  318. 'FINSI' ;
  319.  
  320.  
  321. *'SI' ((('MINI' POROS) < 0.D0) 'OU' (('MAXI' POROS) > 1.D0)) ;
  322. * 'ERREUR' 'La porosité doit etre comprise entre 0 et 1' ;
  323. *'FINSI' ;
  324.  
  325. * RETARD
  326. 'SI' ( 'EXISTE' TRANSI 'COEF_RETARD' ) ;
  327. RETAR1 = TRANSI . 'COEF_RETARD' ;
  328. 'SI' ('EGA' ('TYPE' RETAR1) 'CHPOINT ') ;
  329. RETAR1 = NOMC 'SCAL' RETAR1 ;
  330. MINRET = ('MINI' RETAR1) ;
  331. 'SINON' ;
  332. MINRET = RETAR1 ;
  333. RETAR1 = 'MANU' 'CHPO' MCENT 'SCAL' MINRET ;
  334. RETAR1 = 'CHAN' 'ATTRIBUT' RETAR1 'NATURE' 'DISCRET' ;
  335. 'FINSI' ;
  336. 'SI' (MINRET < 1.D0) ;
  337. 'ERREUR' 'Le coefficient de retard doit etre supérieur à 1';
  338. 'FINSI' ;
  339. 'SINON' ;
  340. * RETAR1 = kcht modhyb scal centre 1.D0 ;
  341. RETAR1 = 'MANU' 'CHPO' MCENT 'SCAL' 1.D0 ;
  342. RETAR1 = 'CHAN' 'ATTRIBUT' RETAR1 'NATURE' 'DISCRET' ;
  343. 'FINSI' ;
  344. RETAR1M1 = RETAR1 - 1. ;
  345. SAUVRET = 'EXISTE' TRANSI 'RETARD' ;
  346. * RETARD NON LINEAIRE
  347. RNONLINL = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
  348. ('EXISTE' TRANSI 'LANGMUIR') ;
  349. RNONLINF = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
  350. ('EXISTE' TRANSI 'FREUNDLICH') 'ET' ('NON' RNONLINL) ;
  351. RNONLIN = RNONLINL 'OU' RNONLINF ;
  352. * Cas Langmuir :
  353. 'SI' RNONLINL ;
  354. FSAT = TRANSI . 'LANGMUIR' ;
  355. 'SI' ('EGA' ('TYPE' FSAT) 'CHPOINT ') ;
  356. 'SI' (('MINI' FSAT) < 0.D0) ;
  357. 'ERREUR'
  358. 'La masse adsorbee a saturation doit etre positive' ;
  359. 'QUITTER' TRANSGEN ;
  360. 'FINSI' ;
  361. 'SINON' ;
  362. 'SI' (FSAT < 0.D0) ;
  363. 'ERREUR'
  364. 'La masse adsorbee a saturation doit etre positive' ;
  365. 'QUITTER' TRANSGEN ;
  366. 'FINSI' ;
  367. 'FINSI' ;
  368. RM1SURF = RETAR1M1 / FSAT ;
  369. 'FINSI' ;
  370. * Cas Freundlich :
  371. 'SI' RNONLINF ;
  372. UNSURN = TRANSI . 'FREUNDLICH' ;
  373. 'SI' ((UNSURN &lt;EG 0.D0) 'OU' (UNSURN > 1.D0)) ;
  374. 'ERREUR'
  375. 'L exposant de Freundlich doit etre compris entre 0'
  376. ' (exclu) et 1' ;
  377. 'QUITTER' TRANSGEN ;
  378. 'FINSI' ;
  379. 'FINSI' ;
  380. * DECROISSANCE
  381. DECRO = ( 'EXISTE' TRANSI 'DECROISSANCE' ) ;
  382. 'SI' DECRO ;
  383. LAMBD0 = TRANSI . 'DECROISSANCE' ;
  384. 'SI' (LAMBD0 < 0.D0) ;
  385. 'ERREUR' 'Le coefficient de décroissance doit etre positif.';
  386. 'FINSI' ;
  387. 'SINON' ;
  388. LAMBD0 = 0. ;
  389. 'FINSI' ;
  390. * DISSOLUTION_IMPOSEE
  391. SOLUA = ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE' ) ;
  392. 'SI' SOLUA ;
  393. DISIMP = TRANSI . 'DISSOLUTION_IMPOSEE' ;
  394. 'FINSI' ;
  395. * LIMITE_SOLUBILITE
  396. * (Priorité de la dissolution imposée sur les autres processus)
  397. SOLUL = (('EXISTE' TRANSI 'LIMITE_SOLUBILITE' ) 'ET' ('NON' SOLUA)) ;
  398. 'SI' SOLUL ;
  399. 'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHPOINT') ;
  400. LIMSOL = 'NOMC' 'H' TRANSI . 'LIMITE_SOLUBILITE' ;
  401. 'FINSI' ;
  402. 'SI' (('MINI' LIMSOL) < 0.D0) ;
  403. mess 'minimum de limsol' (MINI LIMSOL);
  404. 'ERREUR' 'La limite de solubilité doit etre positive.' ;
  405. 'FINSI' ;
  406. 'FINSI' ;
  407. * COEF_DISSOLUTION
  408. SOLUP = SOLUL 'ET' ('EXISTE' TRANSI 'COEF_DISSOLUTION' ) ;
  409. SOLUI = SOLUL 'ET' ('NON' SOLUP) ;
  410. 'SI' SOLUP ;
  411. codis = TRANSI . 'COEF_DISSOLUTION' ;
  412. 'SI' ('EGA' ('TYPE' CODIS) 'CHPOINT ') ;
  413. CODIS = 'NOMC' 'SCAL' TRANSI . 'COEF_DISSOLUTION' NATU DIFFUS ;
  414. 'SI' (('MINI' CODIS) < 0.D0) ;
  415. 'ERREUR' 'Le coefficient de dissolution doit etre positif.';
  416. 'FINSI' ;
  417. 'SINON' ;
  418. CODIS = TRANSI . 'COEF_DISSOLUTION' ;
  419. CODIS = MANU CHPO MCENT SCAL CODIS ;
  420. 'SI' ((MINI CODIS) < 0.D0) ;
  421. 'ERREUR' 'Le coefficient de dissolution doit etre positif.';
  422. 'FINSI' ;
  423. 'FINSI' ;
  424. 'FINSI' ;
  425.  
  426. 'SI' (SOLUP) ;
  427. * Là où la précipitation est nulle (codis =0) on choisit
  428. * une limite de solubilité élevée à des fins d'optimisation
  429. * algorithmique
  430. 'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHPOINT') ;
  431. dum = 'MASQUE' CODIS INFERIEUR 1.D-30 ;
  432. LIMSOL = ((1.D0 '-' dum) * LIMSOL) '+' ('NOMC' 'H' (1.D3*dum)) ;
  433. 'FINSI' ;
  434. 'FINSI' ;
  435. SOLUB = SOLUL 'OU' SOLUA ;
  436. *
  437. * ------------------------------------------------------------------ *
  438. * VERIFICATION DE LA BONNE STRUCTURE DES TABLES RESULTAT D'ORIGINE *
  439. * ------------------------------------------------------------------ *
  440. *
  441. * TEMPS
  442. 'SI' ( 'NON' ('EXISTE' TRANSI 'TEMPS' ) ) ;
  443. 'ERREUR'
  444. 'Indice TEMPS absent de la table de données.' ;
  445. 'QUITTER' TRANSGEN ;
  446. 'FINSI' ;
  447. * INCONNUE
  448. 'SI' ( 'NON' ('EXISTE' TRANSI NOMINC ) ) ;
  449. 'ERREUR' ('CHAINE'
  450. 'Indice ' NOMINC ' absent de la table de données.') ;
  451. 'QUITTER' TRANSGEN ;
  452. 'FINSI' ;
  453. * FLUXDIFF
  454. 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUXDIFF' ) ) ;
  455. 'ERREUR'
  456. 'Indice FLUX diffusif absent de la table de données.' ;
  457. 'QUITTER' TRANSGEN ;
  458. 'FINSI' ;
  459. * FLUXCONV
  460. 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUXCONV' ) ) ;
  461. 'ERREUR'
  462. 'Indice FLUX convectif absent de la table de données.' ;
  463. 'QUITTER' TRANSGEN ;
  464. 'FINSI' ;
  465. * DISSOLUTION, PRECIPITE
  466. 'SI' SOLUB ;
  467. * GBM mettre a jour dissolution .
  468. * 'SI' ( 'NON' ('EXISTE' TRANSI 'DISSOLUTION' ) ) ;
  469. * 'ERREUR'
  470. * 'Indice DISSOLUTION absent de la table de données.' ;
  471. * 'QUITTER' TRANSGEN ;
  472. * 'FINSI' ;
  473. 'SI' ( 'NON' ('EXISTE' TRANSI 'PRECIPITE' ) ) ;
  474. 'ERREUR'
  475. 'Indice PRECIPITE absent de la table de données.' ;
  476. 'QUITTER' TRANSGEN ;
  477. 'FINSI' ;
  478. 'FINSI' ;
  479. * TEST DES TAILLES DE TABLE
  480. IND1 = 'INDEX' ( TRANSI . 'TEMPS' ) ;
  481. IND3 = 'INDEX' ( TRANSI . NOMINC ) ;
  482. IND4 = 'INDEX' ( TRANSI . 'FLUXDIFF' ) ;
  483. LIN1 = 'DIME' IND1 ;
  484. LIN3 = 'DIME' IND3 ;
  485. LIN4 = 'DIME' IND4 ;
  486. 'SI' ( LIN1 'NEG' LIN3 ) ;
  487. 'ERREUR' ('CHAINE'
  488. 'Longueur des tables TEMPS et ' NOMINC ' différente.') ;
  489. 'QUITTER' TRANSGEN ;
  490. 'FINSI' ;
  491. 'SI' ( LIN1 'NEG' LIN4 ) ;
  492. 'ERREUR'
  493. 'Longueur des tables TEMPS et FLUXDIFF différente.' ;
  494. 'QUITTER' TRANSGEN ;
  495. 'FINSI' ;
  496. 'SI' SOLUB ;
  497. * GBM mettre a jour dissolution .
  498. * IND5 = 'INDEX' ( TRANSI . 'DISSOLUTION' ) ;
  499. IND6 = 'INDEX' ( TRANSI . 'PRECIPITE' ) ;
  500. * LIN5 = 'DIME' IND5 ;
  501. LIN6 = 'DIME' IND6 ;
  502. * 'SI' ( LIN1 'NEG' LIN5 ) ;
  503. * 'ERREUR'
  504. * 'Longueur des tables TEMPS et DISSOLUTION différente.' ;
  505. * 'QUITTER' TRANSGEN ;
  506. * 'FINSI' ;
  507. 'SI' ( LIN1 'NEG' LIN6 ) ;
  508. 'ERREUR'
  509. 'Longueur des tables TEMPS et PRECIPITE différente.' ;
  510. 'QUITTER' TRANSGEN ;
  511. 'FINSI' ;
  512. 'FINSI' ;
  513. * TEST DES INDICES DE TABLE
  514. IPO1 = 0 ;
  515. 'REPETER' BOU1 LIN1 ;
  516. IPO1 = IPO1 + 1 ;
  517. LAST1 = IND1 . IPO1 ;
  518. LAST3 = IND3 . IPO1 ;
  519. LAST4 = IND4 . IPO1 ;
  520. 'SI' SOLUB ;
  521. * GBM mettre a jour dissolution .
  522. * LAST5 = IND5 . IPO1 ;
  523. LAST6 = IND6 . IPO1 ;
  524. 'FINSI' ;
  525. 'SI' ( 'NEG' LAST1 LAST3 ) ;
  526. 'ERREUR'
  527. 'Indices des tables TEMPS et CONCENTRATION incohérents.' ;
  528. 'QUITTER' TRANSGEN ;
  529. 'FINSI' ;
  530. 'SI' ( 'NEG' LAST1 LAST4 ) ;
  531. 'ERREUR'
  532. 'Indices des tables TEMPS et FLUX incohérents.' ;
  533. 'QUITTER' TRANSGEN ;
  534. 'FINSI' ;
  535. 'SI' SOLUB ;
  536. * GBM mettre a jour dissolution .
  537. * 'SI' ( 'NEG' LAST1 LAST5 ) ;
  538. * 'ERREUR'
  539. * 'Indices des tables TEMPS et DISSOLUTION incohérents.';
  540. * 'QUITTER' TRANSGEN ;
  541. * 'FINSI' ;
  542. 'SI' ( 'NEG' LAST1 LAST6 ) ;
  543. 'ERREUR'
  544. 'Indices des tables TEMPS et PRECIPITE incohérents.' ;
  545. 'QUITTER' TRANSGEN ;
  546. 'FINSI' ;
  547. 'FINSI' ;
  548. 'FIN' BOU1 ;
  549. *
  550. * ------------------------------------------------------------------ *
  551. * RECUPERATION DES CONDITIONS INITIALES OU DU DERNIER PAS SAUVE *
  552. * ------------------------------------------------------------------ *
  553. *
  554. TPSINI = TRANSI . 'TEMPS' . LAST1 ;
  555.  
  556. *Momentanée pour la trace
  557.  
  558. * On extrait les noms de composante des espèces
  559. NOMESP = 'EXTRAIRE' ('EXTRAIRE' TRANSI . NOMINC . LAST1 'COMP') 1 ;
  560.  
  561. * nomespl est un nom d'espece locale à la procédure car
  562. * les EFMH veulent voir une composante 'H'
  563. NOMESPL = 'MOT' 'H' ;
  564.  
  565. CHRG = 'NOMC' 'H' TRANSI . NOMINC . LAST1 ;
  566. *VERRUE
  567. FLU0 = 'NOMC' 'H' TRANSI . 'FLUXDIFF' . LAST1 ;
  568. FLUCO0 = 'NOMC' 'H' TRANSI . 'FLUXCONV' . LAST1 ;
  569.  
  570.  
  571. *
  572. *--------------------------------------------------------------------*
  573. * RECUPERATION DES THETA SCHEMAS DIFFUSION-CONVECTION *
  574. *--------------------------------------------------------------------*
  575. * THETA DIFFUSION
  576. 'SI' ( 'EXISTE' TRANSI 'THETA_DIFF' ) ;
  577. TETA = TRANSI . 'THETA_DIFF' ;
  578. 'SINON' ;
  579. TETA = 1.D0 ;
  580. 'FINSI' ;
  581. * THETA CONVECTION
  582. 'SI' ( 'EXISTE' TRANSI 'THETA_CONVECTION' ) ;
  583. TETAC = TRANSI . 'THETA_CONVECTION' ;
  584. 'SINON' ;
  585. TETAC = TETA ;
  586. 'FINSI' ;
  587. *
  588. *--------------------------------------------------------------------*
  589. * RECUPERATION DES CONDITIONS AUX LIMITES ET DES CHARGEMENTS *
  590. *--------------------------------------------------------------------*
  591. *
  592. * CL de type Dirichlet : BLOCAGE
  593. * TRACE_IMPOSE
  594.  
  595. * On débranche les tests sur les longueurs de flux non nuls
  596. * en CL pour les cas ou ils ne sont pas intégrés en temps
  597. startflu = 1.D30;
  598. startflt = 1.D30;
  599. startmix = 1.D30;
  600. 'SI' ( 'EXISTE' TRANSI 'TRACE_IMPOSE') ;
  601. CHDIRI = TRANSI . 'TRACE_IMPOSE' ;
  602. 'FINSI' ;
  603. * CL de type Neumann : FLUX_IMPOSE
  604. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' ) ;
  605. FLUIMP = TRANSI . 'FLUX_IMPOSE' ;
  606. 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  607. 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  608. 'ET' ('EGA' (teta * tetac) 1.D0))) ;
  609. * On écrit les flux sous forme intégrale sauf en explicite et
  610. * kranck-Nickholson pour EFMH
  611. FLUIMP startflu = CHAMINT FLUIMP ;
  612. 'FINSI' ;
  613. 'FINSI' ;
  614. * CL de type flux total: FLUXTOT_IMP
  615. 'SI' ( 'EXISTE' TRANSI 'FLUXTOT_IMP' ) ;
  616. FLTOTIMP = TRANSI . 'FLUXTOT_IMP' ;
  617. 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  618. 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  619. 'ET' ('EGA' (teta * tetac) 1.D0))) ;
  620. * On écrit les flux sous forme intégrale sauf en explicite et
  621. * kranck-Nickholson pour EFMH
  622. FLTOTIMP startflt = CHAMINT FLTOTIMP ;
  623. 'FINSI' ;
  624. 'FINSI' ;
  625. * CL de type mixte : MIXTES
  626. 'SI' ( 'EXISTE' TRANSI 'MIXTES' ) ;
  627. FLUMIX = TRANSI . 'MIXTES' ;
  628. 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  629. 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  630. 'ET' ('EGA' (teta * tetac) 1.D0))) ;
  631. * On écrit les flux sous forme intégrale sauf en explicite et
  632. * kranck-Nickholson pour EFMH
  633. FLUMIX startmix = CHAMINT FLUMIX ;
  634. 'FINSI' ;
  635. 'FINSI' ;
  636.  
  637. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  638. TERSOU = TRANSI . 'SOURCE' ;
  639. * on crée une evolution intégrée temporellement pour assurer la
  640. * conservation lors de liste de pas de temps de chargement differente
  641. * des temps de discrétisation
  642. TERSOU startsou = CHAMINT TERSOU ;
  643. 'FINSI' ;
  644.  
  645.  
  646. *
  647. *--------------------------------------------------------------------*
  648. * RECUPERATION DES données NUMERIQUES *
  649. *--------------------------------------------------------------------*
  650. *
  651. * PARAMETRES DE SOLUBILITE
  652. 'SI' (SOLUI 'OU' SOLUP) ;
  653. 'SI' ( 'EXISTE' TRANSI 'ITMAX_LIM') ;
  654. ITMAXI = TRANSI . 'ITMAX_LIM' ;
  655. 'SINON' ;
  656. ITMAXI = 50 ;
  657. 'FINSI' ;
  658. 'FINSI' ;
  659. 'SI' SOLUI ;
  660. 'SI' ( 'EXISTE' TRANSI 'PENALISATION') ;
  661. PENAL = TRANSI . 'PENALISATION' ;
  662. MET0 = 'PENALISATION' ;
  663. 'SINON' ;
  664. 'SI' ( 'EXISTE' TRANSI 'EPSI_LIM') ;
  665. EPS1 = TRANSI . 'EPSI_LIM' ;
  666. MET0 = 'PF DISSOLUTION' ;
  667. 'SINON' ;
  668. 'MESS' 'Il manque le paramètre ' ;
  669. 'ERREUR' 'du schéma numérique de limite de solubilité';
  670. 'FINSI' ;
  671. 'FINSI' ;
  672. 'FINSI' ;
  673. * PARAMETRES DE RETARD NON LINEAIRE
  674. 'SI' RNONLIN ;
  675. 'SI' ( 'EXISTE' TRANSI 'EPSI_RET' ) ;
  676. EPSRNL = TRANSI . 'EPSI_RET' ;
  677. 'SINON' ;
  678. EPSRNL = 1.D-4 ;
  679. 'FINSI' ;
  680. 'SI' ( 'EXISTE' TRANSI 'ITMAX_RET' ) ;
  681. ITMAXRNL = TRANSI . 'ITMAX_RET' ;
  682. 'SINON' ;
  683. ITMAXRNL = 20 ;
  684. 'FINSI' ;
  685. * Petit saut de concentration minimal pour le calcul de dérivée
  686. 'SI' ( 'EXISTE' TRANSI 'EPSI_COR' ) ;
  687. EPSCORD = TRANSI 'EPSI_COR' ;
  688. 'SINON' ;
  689. 'SI' RNONLINL ;
  690. EPSCORD = ('MANU' 'CHPO' MCENT 1 'SCAL' 1.D-4)
  691. / RETAR1 * FSAT ;
  692. 'FINSI' ;
  693. 'SI' RNONLINF ;
  694. EPSCORD = 1.D-4 ;
  695. 'FINSI' ;
  696. 'FINSI' ;
  697. * Pour les isothermes linéaires la boucle de prise en compte de
  698. * retard non linéaire n'est parcourue qu'une fois :
  699. NPicard = ITMAXRNL ;
  700. 'SINON' ;
  701. NPicard = 1 ;
  702. 'FINSI' ;
  703.  
  704. * THETA_DECROISSANCE
  705. 'SI' ( 'EXISTE' TRANSI 'THETA_DEC' ) ;
  706. BETA = TRANSI . 'THETA_DEC' ;
  707. 'SINON' ;
  708. BETA = 0.5 ;
  709. 'FINSI' ;
  710. * THETA_DISSOLUTION
  711. 'SI' ( 'EXISTE' TRANSI 'THETA_DIS' ) ;
  712. GAMMA = TRANSI . 'THETA_DIS' ;
  713. 'SINON' ;
  714. GAMMA = 1. ;
  715. 'FINSI' ;
  716. * TEMPS_CALCULES
  717. 'SI' ( 'EXISTE' TRANSI 'TEMPS_CALCULES' ) ;
  718. TPCAL = 'ORDO' (TRANSI . 'TEMPS_CALCULES') ;
  719. TRANSI . 'TEMPS_CALCULES' = TPCAL ;
  720. DCAL = 'DIME' TPCAL ;
  721. TPSFIN = 'EXTR' TPCAL DCAL ;
  722. TPS0 = 'EXTR' TPCAL 1 ;
  723. 'SI' ('NEG' DCAL 1) ;
  724. TPSANT = 'EXTR' TPCAL (DCAL - 1) ;
  725. TPS1 = 'EXTR' TPCAL 2 ;
  726. DT10 = TPS1 - TPS0 ;
  727. DT21 = TPSFIN - TPSANT ;
  728. PROPS = 'PROG' DT10 DT21 ;
  729. DTEPS = 'MINIMUM' PROPS ;
  730. 'SINON' ;
  731. DTEPS = TPSFIN - TPS0 ;
  732. 'FINSI' ;
  733. *
  734. EPS0 = DTEPS * 1.D-6 ;
  735. * ICAL est le premier indice des temps à calculer supérieur à
  736. * TpsIni + eps
  737. IOK1 = 0 ;
  738. 'REPETER' BOU2 DCAL ;
  739. ICAL = &BOU2 ;
  740. TEMS = 'EXTR' TPCAL ICAL ;
  741. 'SI' ( TPSINI '<' (TEMS - EPS0) ) ;
  742. IOK1 = 1 ;
  743. 'QUITTER' BOU2 ;
  744. 'FINSI' ;
  745. 'FIN' BOU2 ;
  746. 'SI' ( IOK1 'EGA' 0) ;
  747. 'MESS' 'Listreel des temps de calcul inférieur à tini.' ;
  748. 'MESS' 'On ne fait donc rien. ' ;
  749. 'QUITTER' TRANSGEN ;
  750. 'FINSI' ;
  751. DELTAT0 = TEMS - TPSINI ;
  752. 'SINON' ;
  753. 'ERREUR'
  754. 'Indice TEMPS_CALCULES absent de la table de données.' ;
  755. 'QUITTER' TRANSGEN ;
  756. 'FINSI' ;
  757. * TEMPS_SAUVES
  758. 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ;
  759. TPSOR = 'ORDO' (TRANSI . 'TEMPS_SAUVES') ;
  760. TRANSI . 'TEMPS_SAUVES' = TPSOR ;
  761. DSOR = 'DIME' TPSOR ;
  762. ISOR = 0 ;
  763. IOK1 = 0 ;
  764. 'REPETER' BOU3 DSOR ;
  765. ISOR = ISOR + 1 ;
  766. TEMS = 'EXTR' TPSOR ISOR ;
  767. 'SI' ( ( TPSINI '<' (TEMS - EPS0) ) 'ET'
  768. ( TEMS '&lt;EG' (TPSFIN +EPS0)) ) ;
  769. IOK1 = 1 ;
  770. 'QUITTER' BOU3 ;
  771. 'FINSI' ;
  772. 'FIN' BOU3 ;
  773. 'SI' ( IOK1 'EGA' 0) ;
  774. 'ERREUR' 'Listreel des temps de sauvegarde hors tmin,tmax.';
  775. 'QUITTER' TRANSGEN ;
  776. 'FINSI' ;
  777. 'FINSI' ;
  778. *
  779. *--------------------------------------------------------------------*
  780. * RAPPEL DES DONNEES ENTREES DANS LA PROCEDURE *
  781. *--------------------------------------------------------------------*
  782. *
  783. *
  784. 'SI' (EXISTE TRANSI 'AFFICH') ;
  785. affich = TRANSI . 'AFFICH' ;
  786. SINON;
  787. affich = VRAI;
  788. 'FINSI' ;
  789.  
  790. 'SI' (affich) ;
  791. 'MESS' ' ' ;
  792. 'MESS' 'MODELISATION Darcy en transitoire.'
  793. (TRANSI . 'TYPDISCRETISATION') ;
  794. 'MESS' '--------------------------------------------' ;
  795. 'MESS' ' ' ;
  796. 'MESS' 'Donnees présentes en entrée : ' ;
  797. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  798. 'MESS' 'Calcul avec conditions aux limites sur Th (Dirichlet)' ;
  799. 'FINSI' ;
  800. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' ) ;
  801. 'MESS' 'Calcul avec conditions aux limites de flux (Neumann)' ;
  802. 'FINSI' ;
  803. 'SI' ( 'EXISTE' TRANSI 'FLUXTOT_IMP' ) ;
  804. 'MESS' 'Calcul avec conditions aux limites de flux total)' ;
  805. 'FINSI' ;
  806. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' ) ;
  807. 'MESS' 'Ce problème comporte un terme convectif' ;
  808. 'FINSI' ;
  809. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  810. 'MESS' 'Ce problème comporte un terme source' ;
  811. 'FINSI' ;
  812. 'SI' DECRO ;
  813. 'MESS' 'Ce problème comporte un terme de décroissance' ;
  814. 'FINSI' ;
  815. 'SI' RNONLIN ;
  816. 'MESS' 'Ce problème comporte la prise en compte d isothermes '
  817. 'd adsorption non-linéaires' ;
  818. 'SI' RNONLINL ;
  819. 'MESS' ' de type Langmuir (méthode de Picard)' ;
  820. 'FINSI' ;
  821. 'SI' RNONLINF ;
  822. 'MESS' ' de type Freundlich (méthode de Picard)' ;
  823. 'FINSI' ;
  824. 'FINSI' ;
  825. 'SI' SOLUL ;
  826. 'MESS' 'Ce problème comporte une limite de solubilité' ;
  827. 'SI' SOLUP ;
  828. 'MESS' ' avec dissolution progressive d ordre 1' ;
  829. 'FINSI' ;
  830. 'SI' (EGA MET0 'PENALISATION') ;
  831. 'MESS' ' avec coefficient de pénalisation = ' PENAL ;
  832. 'FINSI' ;
  833. 'SI' (EGA MET0 'PF DISSOLUTION') ;
  834. 'MESS' ' avec critère de convergence = ' EPS1 ;
  835. 'FINSI' ;
  836. 'FINSI' ;
  837. 'SI' SOLUA ;
  838. 'MESS' 'Ce probleme comporte une dissolution imposée' ;
  839. 'FINSI' ;
  840. 'MESS' ' ' ;
  841. 'MESS' 'Valeur des paramètres des schémas numériques :' ;
  842. 'MESS' ' (0:Schéma explicite, 0.5:Crank-Nicholson, 1:Implicite)' ;
  843. 'MESS' ' diffusion : ' TETA ;
  844. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' ) ;
  845. 'MESS' ' convection : ' TETAC ;
  846. 'FINSI' ;
  847. 'SI' DECRO ;
  848. 'MESS' ' décroissance : ' BETA ;
  849. 'FINSI' ;
  850. 'SI' SOLUP ;
  851. 'MESS' ' dissolution d ordre 1 : ' GAMMA ;
  852. 'FINSI' ;
  853. 'MESS' ' ' ;
  854. 'MESS' 'Valeur du temps initial : ' TPSINI ;
  855. 'MESS' 'Valeur du temps final : ' TPSFIN ;
  856. 'FINSI' ;
  857. *
  858. *
  859. *--------------------------------------------------------------------*
  860. * BOUCLE RESOLVANT LE SYSTEME POUR CHAQUE PAS DE TEMPS *
  861. *--------------------------------------------------------------------*
  862. *
  863. *==============
  864. * Préparation préliminaire de la boucle sur les pas de temps :
  865. *==============
  866. *
  867. NBIT = 0 ;
  868. NBITRNL = 0 ;
  869. PRECED = TPSINI ;
  870. IPAS = ICAL ;
  871. DELOLD = 0.D0 ;
  872. TPS1 = 'EXTR' TPCAL ICAL ;
  873. DELTAT = TPS1 - TPSINI ;
  874. * initialisation de l'indice de stockage des intégrales de concentrat
  875. LASTINT = LAST1 ;
  876. 'MESS' 'Incrément de temps initial : ' DELTAT ;
  877. 'MESS' ' ' ;
  878.  
  879. *
  880. * initialisation du terme source intégral.
  881. *
  882.  
  883. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  884. TERSC2M1 = 'NOMC' (NOMESPL) ('TIRE' TERSOU TPSINI) ;
  885. 'FINSI' ;
  886.  
  887. 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE') ;
  888. FLUIMPM1 = 'TIRE' FLUIMP TPSINI ;
  889. 'FINSI' ;
  890.  
  891. 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP') ;
  892. FLUTMPM1 = 'TIRE' FLTOTIMP TPSINI ;
  893. 'FINSI' ;
  894.  
  895. 'SI' ('EXISTE' TRANSI 'MIXTES') ;
  896. FLUMMPM1 = 'TIRE' FLUMIX TPSINI ;
  897. 'FINSI' ;
  898.  
  899. *
  900. * Initialisation du materiau - diffusivité
  901. *
  902.  
  903. 'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHARGEME') ;
  904. MAT1 = 'TIRE' TRANSI . 'CARACTERISTIQUES' TPSINI ;
  905. 'FINSI' ;
  906.  
  907.  
  908. *
  909. * Initialisation de la quantité dissoute
  910. *
  911. 'SI' SOLUB ;
  912. PRCI = 'NOMC' 'H' TRANSI . 'PRECIPITE' . LAST1 ;
  913. ANT0 = 'EXCO' 'H' ('MASQ' PRCI 'SUPERIEUR' 1.D-14) 'SCAL' ;
  914. 'FINSI' ;
  915. *
  916. * Calcul de la fraction adsorbée initiale :
  917. *
  918. CC0 = 'NOMC' CHRG 'SCAL' ;
  919. 'SI' RNONLIN ;
  920. 'SI' RNONLINL ;
  921. * Langmuir :
  922. FF0 = RETAR1M1 * CC0
  923. /(1.D0 + ( RM1SURF * CC0 )) ;
  924. 'FINSI' ;
  925. 'SI' RNONLINF ;
  926. * Freundlich :
  927. * La bidouille utilisant les masques sert à obtenir le champ point
  928. * des signes de CC0.
  929. a = 'MASQUE' CC0 'SUPERIEUR' 0. ;
  930. b = 'MASQUE' CC0 'INFERIEUR' 0. ;
  931. FF0 = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b) ;
  932. 'OUBLI' a ;
  933. 'OUBLI' b ;
  934. 'FINSI' ;
  935. 'SINON' ;
  936. * Isotherme linéaire :
  937. FF0 = RETAR1M1 * CC0 ;
  938. 'FINSI' ;
  939. FF = FF0 ;
  940.  
  941. *
  942. * Calcul du coefficient de retard initial :
  943. * Le facteur de retard s'exprime comme
  944. * R(C) = 1 + ( F'(Ct+dt + eps) - F'(Ct) ) / ( Ct+dt + eps - Ct) )
  945. * C'est la méthode de la corde employée ici pour ses vertus
  946. * stabilisantes. Le petit 'eps' (EpsCord) intervient pour évacuer les
  947. * questions de conditions initiales ou de régime permanent.
  948. 'SI' RNONLIN ;
  949. CC1 = CC0 + EpsCord ;
  950. 'SI' RNONLINL ;
  951. * Langmuir
  952. FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 )) ;
  953. 'SINON' ;
  954. * Freundlich
  955. a = 'MASQUE' CC1 'SUPERIEUR' 0. ;
  956. b = 'MASQUE' CC1 'INFERIEUR' 0. ;
  957. FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b) ;
  958. 'DETRUIT' a ;
  959. 'DETRUIT' b ;
  960. 'FINSI' ;
  961. RETARC = 1.D0 + ( (FF1 - FF0) / EpsCord ) ;
  962. 'OUBLI' CC1 ;
  963. 'OUBLI' FF1 ;
  964. 'SINON' ;
  965. RETARC = RETAR1 ;
  966. 'FINSI' ;
  967. RETAR0 = RETARC ;
  968. 'SI' SAUVRET ;
  969. TRANSI . 'RETARD' . 0 = RETARC ;
  970. 'FINSI' ;
  971. *
  972. * Fonctions de volume :
  973. PORETSU = RETARC * POROS ;
  974. DECROI0 = 'EXCO' ( (-1.) * LAMBD0 * PORETSU ) 'SCAL' 'SCAL' ;
  975. DECROI0 = VOLU1 * DECROI0 ;
  976. BDECROI0 = BETA * DECROI0 ;
  977. BMDCROI0 = (1.D0 - BETA) * DECROI0 ;
  978. POROM1V = (1.D0 - POROS) * VOLU1 ;
  979. POROV = POROS * VOLU1 ;
  980. *
  981. * Indicateur les recalculs de matrices à effectuer
  982. *
  983. TABMODI = TABLE ;
  984. TABMODI . 'POROSITE' = VRAI ;
  985. TABMODI . 'CONVECTI' = VRAI ;
  986. TABMODI . 'DELTAT' = VRAI ;
  987. TABMODI . 'COEF_LIN' = VRAI ;
  988. TABMODI . 'DIFFUSIV' = VRAI ;
  989. *
  990. * Initialisation de la table de calcul de trangeol
  991. * pour le calcul d'un itéré de l'équation de transport
  992. *
  993. CHCLIM = 'TABLE' ;
  994. GEOL1 = 'TABLE' ;
  995. * l'option abandon indique qu'en dessous de 10**-13 pour la concentr
  996. * et pour des sources ou flux n'impliquant pas des variation de C
  997. * superieurs à 10**-16 on sort 0 sans calcul
  998. 'SI' ('EXISTE' TRANSI 'SEUILCALC') ;
  999. GEOL1 . 'ABANDON' = VRAI ;
  1000. GEOL1 . 'SEUILCALC' = TRANSI . 'SEUILCALC' ;
  1001. 'SINON' ;
  1002. GEOL1 . 'ABANDON' = FAUX ;
  1003. GEOL1 . 'SEUILCALC' = 1.D-30 ;
  1004. 'FINSI' ;
  1005.  
  1006. 'SI' ('EXISTE' TRANSI 'NUM_PECLET' ) ;
  1007. GEOL1 . 'NUM_PECLET' = TRANSI . 'NUM_PECLET' ;
  1008. 'SINON' ;
  1009. GEOL1 . 'NUM_PECLET' = 2.D0 ;
  1010. 'FINSI' ;
  1011. GEOL1 . 'CONCENTRATION' = CHRG ;
  1012. GEOL1 . 'LUMP' = TRANSI . 'LUMP' ;
  1013. GEOL1 . 'TYPDISCRETISATION' = TRANSI . 'TYPDISCRETISATION' ;
  1014. GEOL1 . 'THETA_DIFFUSION' = TETA ;
  1015. GEOL1 . 'THETA_CONVECTION' = TETAC ;
  1016. GEOL1 . 'DECENTREMENT' = TRANSI . 'DECENTR' ;
  1017. GEOL1 . 'DIFFUSIVITE' = MAT1 ;
  1018. *GEOL1 . 'SOLVEUR' = TRANSI . 'METHINV' . TYPINV ;
  1019. *GEOL1 . 'PRECONDITIONNEUR' = TRANSI . 'METHINV' . PRECOND ;
  1020.  
  1021.  
  1022. GEOL1 . 'MODIFICATI' = TABMODI ;
  1023. GEOL1 . 'POROSITE' = PORETSU ;
  1024. GEOL1 . 'SOURCE' = 'NOMC' NOMESPL
  1025. ('KCHT' (TRANSI.'MODELE') SCAL 'CENTRE' 0.D0) ;
  1026.  
  1027.  
  1028.  
  1029.  
  1030.  
  1031. *=================================
  1032. 'REPETE' BOUTPS (DCAL - ICAL + 1); COMM 'Boucle sur le temps' ;
  1033. *=================================
  1034. *MESS 'TEMPS1' ; TEMPS;
  1035.  
  1036. *
  1037. IPAS = ICAL + &BOUTPS - 1 ;
  1038. *
  1039. *-- Initialisation en-tête de boucle sur le temps
  1040. *
  1041. TPS = 'EXTR' TPCAL IPAS ;
  1042. MESS 'TEMPS T =' TPS ;
  1043.  
  1044. *list ((PLACE) * 4. / 1000000.) ;
  1045.  
  1046.  
  1047. DELTAT = TPS - PRECED ;
  1048. GEOL1 . 'DELTAT' = DELTAT ;
  1049. EPSDT = DELTAT + DELOLD / 2.D0 * 1.D-6 ;
  1050. TABMODI . 'DELTAT' = TABMODI . 'DELTAT'
  1051. 'OU' ( DELOLD 'NEG' DELTAT EPSDT ) ;
  1052.  
  1053. * CONVECTION
  1054. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' ) ;
  1055. GEOL1 . 'CONVECTION' = TRANSI . 'CONVECTION' ;
  1056. 'FINSI' ;
  1057.  
  1058. * VITESSE CENTRE
  1059. 'SI' ( 'EXISTE' TRANSI 'VITELEM' ) ;
  1060. GEOL1 . 'VITELEM' = TRANSI . 'VITELEM' ;
  1061. 'FINSI' ;
  1062.  
  1063. * ALPHAL
  1064. 'SI' ( 'EXISTE' TRANSI 'ALPHAL' ) ;
  1065. GEOL1 . 'ALPHAL' = TRANSI . 'ALPHAL' ;
  1066. 'FINSI' ;
  1067.  
  1068. 'SI' ( 'EXISTE' TRANSI 'ALPHAT') ;
  1069. GEOL1 . 'ALPHAT' = TRANSI . 'ALPHAT' ;
  1070. 'FINSI' ;
  1071.  
  1072. *
  1073. *- Calcul de la contribution des termes sources
  1074. *
  1075. TERSCE = 'MANU' 'CHPO' MCENT 1 'SOUR' 0. ;
  1076. TERSCE = 'NOMC' NOMESPL TERSCE ;
  1077.  
  1078. tpsm1 = 'EXTR' TPCAL (IPAS - 1) ;
  1079. * Terme source propre :
  1080. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  1081. 'SI' (tpsm1 &lt;EG startsou) ;
  1082. * la source varie encore et est non nulle
  1083. zozo = 'TIRE' TERSOU TPS ;
  1084. 'SINON' ;
  1085. * la source est nulle on garde l'intégrale précédente
  1086. zozo = 'COPIER' TERSC2M1 ;
  1087. 'FINSI' ;
  1088. TERSC2 = 'NOMC' NOMESPL zozo ;
  1089. 'DETRUIT' zozo ;
  1090. zozo = TERSC2 '-' TERSC2M1 ;
  1091. TERSCV = zozo '/' DELTAT ;
  1092. 'DETRUIT' zozo ;
  1093. zozo = TERSCV + TERSCE ;
  1094. 'DETRUIT' TERSCE ;
  1095. TERSCE = 'COPIER' zozo ;
  1096. 'DETRUIT' zozo ;
  1097. 'DETRUIT' TERSCV ;
  1098. 'FINSI' ;
  1099.  
  1100. *
  1101. * Terme source de décroissance explicite :
  1102. * du soluté et de l'adsorbat
  1103. 'SI' DECRO ;
  1104. * modif GBM suppression de FF, revient à compter R * R
  1105. * TERSC3 = ( ('NOMC' CHRG 'SCAL') + FF ) * DECROI0 ;
  1106. zozo = 'NOMC' CHRG 'SCAL' ;
  1107. zuzu = zozo * DECROI0 ;
  1108. 'DETRUIT' zozo ;
  1109. TERSC3 = 'NOMC' NOMESPL zuzu ;
  1110. 'DETRUIT' zuzu ;
  1111. zozo = TERSC3 + TERSCE ;
  1112. 'DETRUIT' TERSCE ;
  1113. TERSCE = 'COPIER' zozo ;
  1114. 'DETRUIT' zozo ;
  1115. 'SINON' ;
  1116. * pas de filiation explicite
  1117. TERSC3 = 0.D0 * TERSCE ;
  1118. 'FINSI' ;
  1119.  
  1120.  
  1121. * Chaine de filiation source venant du pere
  1122. 'SI' ('EXISTE' TRANSI 'FILIATION') ;
  1123. * TERSC4 = VOLU1 * (TRANSI . 'FILIATION' . (IPAS '-' 1)
  1124. * '-' TRANSI . 'FILIATION' . (IPAS '-' 2)) '/' DELTAT ;
  1125. zozo = TRANSI . 'FILIATION' . (IPAS '-' 1)
  1126. '-' TRANSI . 'FILIATION' . (IPAS '-' 2) ;
  1127. * GBM attention à cette destruction
  1128. 'DETRUIT' (TRANSI . 'FILIATION' . (IPAS '-' 2)) ;
  1129. zuzu = zozo '/' DELTAT ;
  1130. 'DETRUIT' zozo ;
  1131. zozo = VOLU1 * zuzu ;
  1132. 'DETRUIT' zuzu ;
  1133. TERSC4 = 'NOMC' NOMESPL zozo ;
  1134. 'DETRUIT' zozo ;
  1135. zozo = TERSC4 '+' TERSCE ;
  1136. TERSCE = 'COPIER' zozo ;
  1137. 'DETRUIT' zozo ;
  1138. 'DETRUIT' TERSC4 ;
  1139. 'FINSI' ;
  1140.  
  1141. *
  1142. * Gestion d'un CSAT variable (t) - chargement
  1143. *
  1144. 'SI' ('EXISTE' TRANSI 'LIMITE_SOLUBILITE') ;
  1145. 'SI' ('EGA' ('TYPE' TRANSI . 'LIMITE_SOLUBILITE') 'CHARGEME') ;
  1146. * GBM ON PEUT OPTIMISER EN TESTANT SI LIMSOL CHANGE REELEMENT EN ECART
  1147. * RELATIF PAR RAPPORT AU TEMPS PRECEDENT
  1148. TABMODI . 'LIMSOL' = VRAI ;
  1149. LIMSOL = 'TIRE' TRANSI . 'LIMITE_SOLUBILITE' TPS ;
  1150. 'SI' (SOLUP) ;
  1151. * Là où la précipitation est nulle (codis =0) on choisit
  1152. * une limite de solubilité élevée à des fins d'optimisation
  1153. * algorithmique
  1154. dum = 'MASQUE' CODIS INFERIEUR 1.D-30 ;
  1155. LIMSOL = ((1.D0 '-' dum) * LIMSOL) '+' ('NOMC' 'H' (1.D3*dum)) ;
  1156. 'FINSI' ;
  1157. 'SINON' ;
  1158. TABMODI . 'LIMSOL' = FAUX ;
  1159. 'FINSI' ;
  1160. 'FINSI' ;
  1161.  
  1162. *
  1163. * Gestion de la diffusivité variable en temps - chargement. Donnée
  1164. * explicite, ie au temps n lors du traitement vers n+1
  1165. *
  1166.  
  1167. 'SI' ('EGA' ('TYPE' TRANSI . 'CARACTERISTIQUES') 'CHARGEME') ;
  1168. * On peut optimiser en testant la variation réelle de la diffusion
  1169. * elle pourrait en effet être constante par palliers
  1170. TABMODI . 'DIFFUSIV' = VRAI ;
  1171. MAT1 = 'TIRE' TRANSI . 'CARACTERISTIQUES' TPS ;
  1172. * Estelle : affectation de la diffusion dans geol1
  1173. GEOL1 . 'DIFFUSIVITE' = MAT1 ;
  1174. 'SINON' ;
  1175. TABMODI . 'DIFFUSIV' = FAUX ;
  1176. 'FINSI' ;
  1177.  
  1178. *
  1179. *- Incorporation des CLs
  1180. *
  1181. 'SI' ('EXISTE' TRANSI 'TRACE_IMPOSE') ;
  1182. CHARIMPO = 'TIRE' CHDIRI TPS ;
  1183. CHCLIM . 'DIRICHLET' = 'NOMC' NOMESPL CHARIMPO ;
  1184. 'DETRUIT' CHARIMPO ;
  1185. 'FINSI' ;
  1186.  
  1187. 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE') ;
  1188. 'SI' (tpsm1 '&lt;EG' startflu) ;
  1189. FLUIMPO = 'TIRE' FLUIMP TPS ;
  1190. 'SINON' ;
  1191. FLUIMPO = 'COPIER' FLUIMPM1 ;
  1192. 'FINSI' ;
  1193. FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET' ;
  1194. 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  1195. 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  1196. 'ET' ('EGA' (teta * tetac) 1.D0))) ;
  1197. * On écrit les flux sous forme intégrale sauf en explicite et
  1198. * kranck-Nickholson pour EFMH. On les différentie ici
  1199. zozo = FLUIMPO '-' FLUIMPM1 ;
  1200. FLUMP = zozo '/' DELTAT ;
  1201. 'DETRUIT' zozo ;
  1202. 'SINON' ;
  1203. FLUMP = 'COPIER' FLUIMPO ;
  1204. 'FINSI' ;
  1205. CHCLIM . 'NEUMANN' = 'NOMC' NOMESPL FLUMP ;
  1206. 'DETRUIT' FLUMP ;
  1207. 'FINSI' ;
  1208.  
  1209.  
  1210. 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP') ;
  1211. 'SI' (tpsm1 '&lt;EG' startflt) ;
  1212. FLUTMPO = 'TIRE' FLTOTIMP TPS ;
  1213. 'SINON' ;
  1214. FLUTMPO = 'COPIER' FLUTMPM1 ;
  1215. 'FINSI' ;
  1216. FLUTMPO = 'CHAN' 'ATTRIBUT' FLUTMPO 'NATURE' 'DISCRET' ;
  1217. 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  1218. 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  1219. 'ET' ('EGA' (teta * tetac) 1.D0))) ;
  1220. * On écrit les flux sous forme intégrale sauf en explicite et
  1221. * kranck-Nickholson pour EFMH. On les différentie ici
  1222. zozo = FLUTMPO '-' FLUTMPM1 ;
  1223. FLUTMP = zozo '/' DELTAT ;
  1224. 'DETRUIT' zozo ;
  1225. 'SINON' ;
  1226. FLUTMP = 'COPIER' FLUTMPO ;
  1227. 'FINSI' ;
  1228. CHCLIM . 'FLUTOTAL' = 'NOMC' NOMESPL FLUTMP ;
  1229. 'DETRUIT' FLUTMP ;
  1230. 'FINSI' ;
  1231.  
  1232.  
  1233. 'SI' ('EXISTE' TRANSI 'MIXTES') ;
  1234. 'SI' (tpsm1 '&lt;EG' startmix) ;
  1235. FLUMMPO = 'TIRE' FLUMIX TPS ;
  1236. 'SINON' ;
  1237. FLUMMPO = 'COPIER' FLUMMPM1 ;
  1238. 'FINSI' ;
  1239. FLUMMPO = 'CHAN' 'ATTRIBUT' FLUMMPO 'NATURE' 'DISCRET' ;
  1240. 'SI' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'VF')
  1241. 'OU' (('EGA' (TRANSI . 'TYPDISCRETISATION') 'EFMH')
  1242. 'ET' ('EGA' (teta * tetac) 1.D0))) ;
  1243. * On écrit les flux sous forme intégrale sauf en explicite et
  1244. * kranck-Nickholson pour EFMH. On les différentie ici
  1245. zozo = FLUMMPO '-' FLUMMPM1 ;
  1246. FLUMMP = zozo '/' DELTAT ;
  1247. 'DETRUIT' zozo ;
  1248. 'SINON' ;
  1249. FLUMMP = 'COPIER' FLUMMPO ;
  1250. 'FINSI' ;
  1251. CHCLIM . 'FLUMIXTE' = TABLE ;
  1252. CHCLIM . 'FLUMIXTE' . 'VAL' = 'NOMC' NOMESPL FLUMMP ;
  1253. CHCLIM . 'FLUMIXTE' . 'COEFA' = TRANSI . 'MIXCOFA' ;
  1254. CHCLIM . 'FLUMIXTE' . 'COEFB' = TRANSI . 'MIXCOFB' ;
  1255. 'DETRUIT' FLUMMP ;
  1256. 'FINSI' ;
  1257.  
  1258. GEOL1 . 'CLIMITES' = CHCLIM ;
  1259. * initialisation de la table de préconditionnement et eventuellement
  1260. * des traces de concentration en EFMH au premier temps
  1261. * A virer dès modification arguments transgeol
  1262.  
  1263. 'SI' ('EGA' &BOUTPS 1) ;
  1264. * VERRUE d'initialisation à virer
  1265. * il s'agit d'une avancée de DT = 0.
  1266. GEOL1 . 'DELTAT' = 1.D-15 ;
  1267. GEOL1 . 'METHINV' = TRANSI . 'METHINV' ;
  1268. * Dans ce cas particulier d'un dt petit, precond diagonal, bcgs
  1269. GEOL1 . 'SOLVEUR' = 3 ;
  1270. GEOL1 . 'PRECONDITIONNEUR' = 1 ;
  1271. GEOL1 GEOL2 = TRANGEOL transi.'MODELE' GEOL1 ;
  1272. * On remet les choix utilisateur pour la suite
  1273. GEOL2 . 'METHINV' . 'TYPINV' = TRANSI . 'METHINV' . 'TYPINV' ;
  1274. GEOL2 . 'METHINV' . 'PRECOND' = TRANSI . 'METHINV' . 'PRECOND' ;
  1275. SI (EGA (TRANSI . 'METHINV' . 'PRECOND') 8) ;
  1276. GEOL2 . 'METHINV' . 'ILUTDTOL' = 1.D-2 ;
  1277. FINSI ;
  1278.  
  1279. * remise à jour des paramètres d'entrée de EFMHTgen
  1280. GEOL1 . 'DELTAT' = DELTAT ;
  1281.  
  1282. TABMODI . 'POROSITE' = FAUX ;
  1283. TABMODI . 'CONVECTI' = FAUX ;
  1284. TABMODI . 'COEF_LIN' = FAUX ;
  1285. * TABMODI . DIFFUSIV est écrasé à FAUX meme si chargement car
  1286. * un premier appel à trangeol a déjà été fait pour initialiser les
  1287. * préconditionnements, seul le pas de temps change ici. C'est vrai
  1288. * uniquement pour le premier pas de temps.
  1289. TABMODI . 'DIFFUSIV' = FAUX ;
  1290. TABMODI . 'DELTAT' = VRAI ;
  1291.  
  1292. 'FINSI' ;
  1293.  
  1294. *
  1295. *
  1296. *|-------------------------------------------------------------------|
  1297. *| Boucle de Retard Non Linéaire |
  1298. *|-------------------------------------------------------------------|
  1299. *
  1300. * Boucle de Picard pour la prise en compte
  1301. * d'isothermes d'adsorption non linéaires : isothermes de Langmuir et
  1302. * de Freundlich. Ce point nécessite d'effectuer une boucle sur
  1303. * le facteur de retard qui devient fonction de la concentration
  1304. * en solution.
  1305. *
  1306. *------------------------
  1307. 'REPETE' BOURNL NPicard ; COMM 'Boucle sur le coef. de retard' ;
  1308. *------------------------
  1309. *
  1310. IRNL = &BOURNL ;
  1311. *
  1312. * Quatre grands blocs selon que l'on traite :
  1313. * -------------------------------------------
  1314. * -1- sans limite de solubilité ou avec dissolution arbitraire
  1315. * -2- avec limite de solubilité et dissolution d'ordre 1
  1316. * -3- avec limite suivant méthode de pénalisation,
  1317. * -4- avec limite suivant schéma prédicteur/correcteur itératif
  1318.  
  1319.  
  1320. *|----------------------------------------------------|
  1321. *| Pas de limite de solubilité ou Dissolution imposée |
  1322. *|----------------------------------------------------|
  1323. 'SI' (('NON' SOLUB) 'OU' SOLUA) ;
  1324. *
  1325. * On recalcule la matrice hybride si DeltaT ou le coefdt ont changé :
  1326. 'SI' (TABMODI.'DELTAT' 'OU' TABMODI . 'POROSITE') ;
  1327.  
  1328. * GBM FAIRE MEILLEUR TEST SI DELTAT BOUGE SEULEMENT TABMODI.PORO
  1329. * NE CHANGE PAS
  1330.  
  1331. * Coefficients correspondant à la décroissance expl. et implicite :
  1332. NUM0 = (BETA - 1.) * LAMBD0 * DELTAT + 1. ;
  1333. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1334. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1335. COFDC = PORETSU * DEN0 ;
  1336. * On charge le coef de DC/DT dans GEOL1 pour TRANGEOL
  1337. GEOL1 . 'POROSITE' = COFDC ;
  1338. * On indique que le coef devant DC/DT a changé
  1339. TABMODI . 'POROSITE' = VRAI ;
  1340. 'FINSI' ;
  1341. *
  1342. 'SI' SOLUA ;
  1343. * Evaluation des termes sources de dissolution :
  1344. * Dissolution arbitraire (ordre 0 par exemple), limitée au total
  1345. * précipité :
  1346. DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1 ;
  1347. DISSOLP = 'TIRE' DISIMP TPS ;
  1348. TERSC4 = 0.5 * (DISSOLT + DISSOLP - ('ABS' (DISSOLT - DISSOLP)));
  1349. TERSC4 = 'EXCO' 'H' TERSC4 'H' ;
  1350. TERSCT = TERSCE + TERSC4 ;
  1351. 'SINON' ;
  1352. TERSCT = TERSCE ;
  1353. 'FINSI' ;
  1354.  
  1355. * On charge le terme source dans GEOL1
  1356.  
  1357. GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSCT ;
  1358.  
  1359. * options de calcul
  1360. GEOL1 . 'MODIFICATI' = TABMODI ;
  1361.  
  1362. *- Résolution - APPEL TRANGEOL
  1363. * OPTIMISER STOCKAGE DOUBLE DES MATRICES PRECOND ?????
  1364. * faire transgeol mat1 mat2 et stocké tracini tracfin, conini concfin...
  1365. * on pourrait alors écraser les matrices sytématiquement et faire un
  1366. * advance comme philippe : stocke cfin dans cini
  1367. * transgeol doit pouvoir accepter matrice vide precond en option
  1368. GEOLPF1 GEOLPF2 = TRANGEOL transi.'MODELE' GEOL1 GEOL2 ;
  1369. CHA2 = GEOLPF1 . 'CONCENTRATION' ;
  1370. FLU2 = GEOLPF1 . 'FLUXDIFF' ;
  1371. FLUCO2 = GEOLPF1 . 'FLUXCONV' ;
  1372.  
  1373. * Remise à jour des coef de recalcul des matrices - dans le point
  1374. * fixe, ces propriétés ne changent pas
  1375. TABMODI . 'POROSITE' = FAUX ;
  1376. TABMODI . 'CONVECTI' = FAUX ;
  1377. TABMODI . 'DELTAT' = FAUX ;
  1378. TABMODI . 'COEF_LIN' = FAUX ;
  1379. TABMODI . 'DIFFUSIV' = FAUX ;
  1380.  
  1381. 'SI' SOLUA ;
  1382. DIS1 = 'EXCO' 'H' (TERSC4 * DELTAT) 'H' ;
  1383. PRE2 = ((NUM0 * PRCI) - (DIS1 / VOLU1)) / DEN0 ;
  1384. 'FINSI' ;
  1385.  
  1386.  
  1387. * 'Fin module dissolution arbitraire'
  1388. 'FINSI' ;
  1389.  
  1390. *list ((PLACE) * 4. / 1000000.) ;
  1391.  
  1392. *|---------------------------|
  1393. *| Dissolution d'ordre 1 |
  1394. *|---------------------------|
  1395. 'SI' SOLUP;
  1396.  
  1397. * On recalcule la matrice hybride si DeltaT ou le coefdt ont changé ou
  1398. * la limite de solubilité a changé:
  1399. 'SI' (TABMODI.'DELTAT' 'OU' TABMODI . 'POROSITE'
  1400. 'OU' TABMODI . 'LIMSOL') ;
  1401. * présence de précipitation dissolution
  1402. zozo = 'MASQ' PRCI 'SUPERIEUR' 1.D-14 ;
  1403. ANTP = 'EXCO' 'H' zozo 'SCAL' ;
  1404. 'DETRUIT' zozo ;
  1405.  
  1406. zozo = (1.D0 '+' 1.D-10) * LIMSOL ;
  1407. zuzu = 'MASQ' CHRG 'SUPERIEUR' zozo ;
  1408. 'DETRUIT' zozo ;
  1409. ANTC = 'EXCO' 'H' zuzu 'SCAL' ;
  1410. 'DETRUIT' zuzu ;
  1411.  
  1412. zozo = ANTP '+' ANTC ;
  1413. 'DETRUIT' ANTC ;
  1414. ANT0 = 'MASQ' zozo 'SUPERIEUR' 0.1 ;
  1415. 'DETRUIT' zozo ;
  1416.  
  1417. * Coefficients correspondant à la décroissance expl. et implicite :
  1418. NUM0 = (BETA - 1.) * LAMBD0 * DELTAT + 1. ;
  1419. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1420.  
  1421. * Terme implicite de dissolution, schéma Euler
  1422. zozo = DELTAT * CODIS ;
  1423. DEN1 = 'KOPS' zozo * ANT0 ;
  1424. 'DETRUIT' zozo ;
  1425.  
  1426. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1427. * solub en poros*codis Id et non Retard*poros*codis
  1428. zozo = 'KOPS' PORETSU * DEN0 ;
  1429. zeze = 'KOPS' POROS * DEN1 ;
  1430. COFDC = 'KOPS' zeze '+' zozo ;
  1431. 'DETRUIT' zozo ;
  1432. 'DETRUIT' zeze ;
  1433.  
  1434. * On charge le coef de DC/DT dans GEOL1 pour TRANGEOL
  1435. GEOL1 . 'POROSITE' = COFDC ;
  1436. * On indique que le coef devant DC/DT a changé
  1437. TABMODI . 'POROSITE' = VRAI ;
  1438.  
  1439. 'FINSI' ;
  1440.  
  1441. * ON charge dans une variable point fixe le precipité
  1442. * initial
  1443. PRCIPF = 'COPIER' PRCI ;
  1444. * On initialise le terme correspondant à une dissolution
  1445. * complete du précipité le derniere itération avant sa disparition
  1446. * Pour eviter des précipités négatifs
  1447. SRCPRCI = 0.D0 * PRCI ;
  1448.  
  1449. *---------------
  1450. 'REPETER' BOU8 ; COMM 'Boucle dissolution ordre 1' ;
  1451. *---------------
  1452.  
  1453.  
  1454. * terme source modifié en raison du coef de solubilité :
  1455. * on calcul deltat codis (cn+1 - cn) /deltat, astuce informatique
  1456. * pour créer le terme deprécipitation. Il faut retirer le - codis cn
  1457. * au second membre : - codis * posos * volume (car source integrale)
  1458. * TERSC5 = TERSCE '-' ( CODIS * ANT0 * POROS * VOLU1 * CHRG) ;
  1459. * On a également la source liée au terme codis poros limsol
  1460. * TERSC5 = TERSC5 '+' ( CODIS * ANT0 * POROS * VOLU1 * LIMSOL) ;
  1461.  
  1462. zozo = 'KOPS' CODIS * ANT0 ;
  1463. zeze = 'KOPS' zozo * POROV ;
  1464. 'DETRUIT' zozo ;
  1465. zuzu = LIMSOL '-' CHRG ;
  1466. zozo = zeze * zuzu ;
  1467. 'DETRUIT' zeze ;
  1468. 'DETRUIT' zuzu ;
  1469. TERSC5 = TERSCE '+' zozo ;
  1470. 'DETRUIT' zozo ;
  1471.  
  1472. * On a également la source liée a la dissolution
  1473. * instantanée du précipité là 'OU' il est devenu négatif
  1474. * artifice numérique permettant d'avoir toujours un
  1475. * précipité négatif. Bilan de matiere verifie, erreur en
  1476. * O(DT) su le reste de la physique
  1477. * TERSC5 = TERSC5 '+' (VOLU1 * SRCPRCI) ;
  1478.  
  1479. zozo = VOLU1 * SRCPRCI ;
  1480. zuzu = TERSC5 '+' zozo ;
  1481. 'DETRUIT' zozo ;
  1482. 'DETRUIT' TERSC5 ;
  1483. TERSC5 = 'COPIER' zuzu ;
  1484. 'DETRUIT' zuzu ;
  1485.  
  1486. * le fait que codis puisse etre gigantesque est sans influence
  1487. * sur la précision car le précipité est calulé par des bilans de flux
  1488. * et n'utilisera plus ces valeurs.
  1489.  
  1490.  
  1491. * On charge le terme source dans GEOL1.
  1492.  
  1493. GEOL1 . 'SOURCE' = 'NOMC' NOMESPL TERSC5 ;
  1494.  
  1495. * options de calcul
  1496. GEOL1 . 'MODIFICATI' = TABMODI ;
  1497.  
  1498. *- Résolution - APPEL TRANGEOL
  1499. * OPTIMISER STOCKAGE DOUBLE DES MATRICES PRECOND ?????
  1500. * faire transgeol mat1 mat2 et stocké tracini tracfin, conini concfin...
  1501. * on pourrait alors écraser les matrices sytématiquement et faire un
  1502. * advance comme philippe : stocke cfin dans cini
  1503. * transgeol doit pouvoir accepter matrice vide precond en option
  1504. *list ((PLACE) * 4. / 1000000.) ;
  1505. GEOLPF1 GEOLPF2 = TRANGEOL transi.'MODELE' GEOL1 GEOL2 ;
  1506. *list ((PLACE) * 4. / 1000000.) ;
  1507. CHA2 = GEOLPF1 . 'CONCENTRATION' ;
  1508. FLU2 = GEOLPF1 . 'FLUXDIFF' ;
  1509. FLUCO2 = GEOLPF1 . 'FLUXCONV' ;
  1510.  
  1511. * remise à jour de XINIT pour le point fixe. Pourra etre retirer
  1512. * quand transgeol stockera les état initiaux- finaux et intermédiaire
  1513. * pour un point fixe. Permet uniquement d'accélérer la convergence
  1514. * des méthodes itératives. Sans incidence sur les résultats physiques
  1515. GEOL2 . 'METHINV' . 'XINIT' = GEOLPF2 . 'METHINV' . 'XINIT' ;
  1516. 'SI' ('EXISTE' GEOL2 TRACE_CONC) ;
  1517.  
  1518. * On conserva absolument la trace en EFMH car c'est la vraie
  1519. * inconnue et on doit garder celle du début du point fixe
  1520. * (sinon on incrémente en temps dans le point fixe).
  1521. ** 'MESSAGE' 'yoyououou';
  1522. dum = GEOL2 . 'TRACE_CONC' ;
  1523. GEOL2 = GEOLPF2 ;
  1524. GEOL2 . 'TRACE_CONC' = dum ;
  1525. 'FINSI' ;
  1526.  
  1527. * Remise à jour des coef de recalcul des matrices
  1528. TABMODI . 'POROSITE' = FAUX ;
  1529. TABMODI . 'CONVECTI' = FAUX ;
  1530. TABMODI . 'DELTAT' = FAUX ;
  1531. TABMODI . 'COEF_LIN' = FAUX ;
  1532. TABMODI . 'DIFFUSIV' = FAUX ;
  1533.  
  1534. * SI ANT0 nul partout alors pas de prec dissolution donc le precipité
  1535. * est forcément nul partout. On ne fait pas de calcul.
  1536.  
  1537. 'SI' (('MAXIMUM' ANT0) 'EGA' 0) ;
  1538. * on met à 0 le précipité, on utilise le champ de concentration
  1539. * pour avoir la géométrie du support.
  1540. PRE2 = 0.D0 * CHRG ;
  1541. 'SINON' ;
  1542.  
  1543. * Calcul du précipité par bilan des variations des quantités et des
  1544. * flux diffusifs et convectifs sur la maille. On n'utilise surtout
  1545. * pas la loi de dissolution codis * (C - limsol) car si codis
  1546. * est grand (pénalisation) le calcul est très imprécis.
  1547. * le précipité est en mole/volume de milieu solide.
  1548. * On décompose le calcul en plusieurs termes :
  1549. * Faire ATTENTION rien n'est intuitif. Si des schémas sont modifiés
  1550. * en amont, que les termes sources sont modifiés, il faudra
  1551. * s'assurer que cela ne modifie pas les lignes suivantes.
  1552.  
  1553. * la variation de concentration dans le milieu (et du sorba)
  1554. * poros * retar * (delta COncentration) * volum / deltat
  1555.  
  1556. *list ((PLACE) * 4. / 1000000.) ;
  1557. PRR1 = (-1.D0 '/' DELTAT) * VOLU1 ;
  1558. PRR2 = 'KOPS' PRR1 * PORETSU ;
  1559. 'DETRUIT' PRR1;
  1560. zozo = 'KOPS' CHA2 - CHRG ;
  1561. PRR1 = 'KOPS' PRR2 * zozo ;
  1562. 'DETRUIT' zozo;
  1563. 'DETRUIT' PRR2;
  1564.  
  1565.  
  1566. * la variation de concentration disparue par filiation. Il faut
  1567. * écrire les équations de bilan pour se convaincre que tout n'est pas
  1568. * compté plusieurs fois de trop: - lambda * poros * retard * volume *
  1569. * Cn+1 multiplié par beta (valeur du theta schéma pour décroissance).
  1570. * plkus partie explicite : - lambda * poros * retard * volume *
  1571. * Cn * (1 - beta). Le signe '-' est dans DECROI0
  1572.  
  1573. zozo = 'KOPS' BDECROI0 * CHA2 ;
  1574. PRR2 = 'KOPS' PRR1 '+' zozo ;
  1575. 'DETRUIT' zozo;
  1576. 'DETRUIT' PRR1;
  1577. zozo = 'KOPS' BMDCROI0 * CHRG ;
  1578. PRR1 = 'KOPS' PRR2 '+' zozo;
  1579. 'DETRUIT' zozo;
  1580. 'DETRUIT' PRR2;
  1581.  
  1582.  
  1583. * Le terme source de C qui comprend d'ailleurs, la filiation venant
  1584. * du pere (concentration + précipité cumulés), et les termes
  1585. * source de l'utilisateur issu de l'indice SOURCE de la table
  1586. * d'entrée. La multiplication par le volume de l'élément
  1587. * est déjà effectuée. On a retirer TERSC3 (source de filiation)
  1588. * déjà comptée ligne précédente. Encore une fois faire très
  1589. * attention lors de modif de ces lignes, on peut s'y perdre.
  1590. * A la limite tout coder en Fortran
  1591.  
  1592. zaza = 'KOPS' TERSCE - TERSC3 ;
  1593. PRR2 = 'KOPS' PRR1 + zaza ;
  1594. 'DETRUIT' zaza;
  1595. 'DETRUIT' PRR1 ;
  1596.  
  1597. * le variation par filiation du précipité - partie explicite
  1598. * du theta schéma de décroissance radioactive : - (1 - poros)
  1599. * fois LAMBDA * (1 - BETA) * volume maille * precipité (instant n)
  1600.  
  1601. PRR1 = 'KOPS' POROM1V * PRCIPF ;
  1602. zozo = (LAMBD0 * (1.D0 '-' BETA)) * PRR1 ;
  1603. 'DETRUIT' PRR1 ;
  1604. PRR1 = 'KOPS' PRR2 - zozo ;
  1605. 'DETRUIT' zozo;
  1606. 'DETRUIT' PRR2 ;
  1607.  
  1608. * Les flux diffusifs et convectifs intégrés au travers des faces de
  1609. * l'élément. flux SORTANT donc diminuant le précipité
  1610. * fluttot = (TETA * FLU2) + (1.D0 '-' TETA) * FLU0) '+'
  1611. * ((TETAC * FLUCO2) '+' ((1.D0 '-' TETAC) * FLUCO0)) ;
  1612.  
  1613. zozo = TETA '*' FLU2 ;
  1614. zizi = ((1.D0 '-' TETA) * FLU0) ;
  1615. zaza = 'KOPS' zizi + zozo ;
  1616. 'DETRUIT' zizi;
  1617. 'DETRUIT' zozo;
  1618.  
  1619. zozo = TETAC * FLUCO2 ;
  1620. zizi = ((1.D0 '-' TETAC) * FLUCO0) ;
  1621. zeze = 'KOPS' zizi + zozo;
  1622. 'DETRUIT' zizi;
  1623. 'DETRUIT' zozo;
  1624.  
  1625. fluttot1 = 'KOPS' zeze + zaza;
  1626. 'DETRUIT' zeze;
  1627. 'DETRUIT' zaza;
  1628.  
  1629. * DIVU veut une composante FLUX
  1630.  
  1631. fluttot = 'NOMC' 'FLUX' fluttot1 ;
  1632. 'DETRUIT' fluttot1;
  1633.  
  1634. divflt = 'DIVU' MODDARCY fluttot MCHYB ;
  1635. 'DETRUIT' fluttot;
  1636.  
  1637. divfll = 'NOMC' 'SCAL' divflt;
  1638. 'DETRUIT' divflt;
  1639.  
  1640. PRR2 = 'KOPS' PRR1 '-' divfll ;
  1641. 'DETRUIT' divfll;
  1642. 'DETRUIT' PRR1;
  1643.  
  1644.  
  1645.  
  1646. * La quantité de précipité à l'instant précédent :
  1647. * (1 '-' POROS) PRCI * VOLU1 '/' DELTAT
  1648. * AFIN d'éviter les problemes lorsque DELTAT = 0
  1649. * on multiplie tout par DELTAT, donc ne pas le faire 'PLUS' loin
  1650. * PRE2 = (deltat * PRE2) '+' ((1.D0 '-' POROS) * VOLU1 * PRCIPF) ;
  1651.  
  1652. PRR1 = deltat * PRR2 ;
  1653. 'DETRUIT' PRR2 ;
  1654. zozo = 'KOPS' POROM1V * PRCIPF ;
  1655. PRR2 = 'KOPS' PRR1 '+' zozo ;
  1656. 'DETRUIT' zozo;
  1657. 'DETRUIT' PRR1;
  1658.  
  1659.  
  1660. * On résout alors (vol * (1 - poros) /deltat) + (1 - poros) vol lambda
  1661. * beta le tout fois precipité instant N+1 = PRE2 calculé au dessus
  1662. *
  1663. * coff = ((DELTAT * LAMBD0 * BETA) '+' 1.D0) * (1.D0 '-' POROS)
  1664. * * VOLU1 ;
  1665. * PRE2 = PRE2 '/' coff ;
  1666.  
  1667.  
  1668. coff = ((DELTAT * LAMBD0 * BETA) '+' 1.D0) * POROM1V;
  1669. PRR1 = 'KOPS' PRR2 '/' coff;
  1670. 'DETRUIT' coff;
  1671. 'DETRUIT' PRR2;
  1672.  
  1673.  
  1674. * ATTENTION SI IL Y A DU PR2CIPITE, POROSITE = 1 INTERDIT CAR
  1675. * DIT QU'il NY A PAS DE SOLIDE
  1676.  
  1677. * On multiplie par la fonction caractéristique de présence du précipité
  1678. * inutile en pratique si les schémas sont conservatifs sur l'élément
  1679. * c'est mis pour ne pas propager d'erreurs machines.
  1680.  
  1681. PRE2 = 'KOPS' ANT0 * PRR1 ;
  1682. 'DETRUIT' PRR1 ;
  1683.  
  1684.  
  1685. zozo = 'NOMC' NOMESPL PRE2 ;
  1686. PRE2 = 'COPIER' zozo;
  1687. 'DETRUIT' zozo ;
  1688.  
  1689. * Fin du calcul du précipité par bilan
  1690. 'FINSI' ;
  1691.  
  1692. * mess 'maximum precipite' (maxi pre2);
  1693. * mess 'maximum concentration' (maxi cha2);
  1694.  
  1695. * Critère de sortie et mise à jour de la nouvelle distribution :
  1696. * On sort la distribution de précipité n'a pas changé (ie le précpité
  1697. * est toujours positif aux memes endoits) et que la concentration ne
  1698. * dépasse pas limsol dans la partie ou il n'y a pas de précipité
  1699.  
  1700. * LTI2 = 'CHAINE' 'Front 1D-h temps ' ;
  1701. * 'TITR' LTI2 ;
  1702. * drmil = 'QUELCONQUE' 'SEG2' ('EXTRAIRE' CHA2 maillage) ;
  1703. * drmil = 'INVERSE' drmil ;
  1704. * AV2 = 'EVOL' 'ROUG' 'CHPO' CHA2 'H' drmil ;
  1705. * 'DESS' AV2 'MIMA' ;
  1706. * AV2 = 'EVOL' 'ROUG' 'CHPO' PRE2 'H' drmil ;
  1707. * 'DESS' AV2 'MIMA' 'NCLK' ;
  1708.  
  1709.  
  1710. * ANT2 caracterise la présence de dissolution quant vaut 1
  1711. zozo = 'MASQ' PRE2 'SUPERIEUR' 1.D-14 ;
  1712. ANT2 = 'EXCO' 'H' zozo 'SCAL' ;
  1713. 'DETRUIT' zozo ;
  1714. * ANT3 caracterise la précipitation quand vaut 1
  1715. * On autorise une seuil supérieur legerement car on relache du
  1716. * precipité en solution dans le point fixe
  1717. zozo = 'MASQ' CHA2 'SUPERIEUR' ((1.D00 '+' 1.D-7) * LIMSOL) ;
  1718. ANT3 = 'EXCO' 'H' zozo 'SCAL' ;
  1719. 'DETRUIT' zozo ;
  1720.  
  1721.  
  1722. * 'MESSAGE' 'maxi cha2' ('MAXIMUM' cha2);
  1723. * 'MESSAGE' 'maxi ant2' ('MAXIMUM' ant2);
  1724. * 'MESSAGE' 'maxi ant3' ('MAXIMUM' ant3);
  1725.  
  1726. * 'SI' (TPS 'EGA' 8.0D5 1.) ;
  1727. * DEPLA (DOMA moddarcy 'MAILLAGE') AFFI (1./15.) (0. 0.) (1. 0.) ;
  1728. * toto = kcha MODDARCY pre2 'CHAM';
  1729. * trac MODDARCY toto;
  1730. * toto = kcha MODDARCY ANT2 'CHAM';
  1731. ** trac MODDARCY toto;
  1732. * tutu = kcha MODDARCY ANT3 'CHAM';
  1733. * trac MODDARCY (toto '-' tutu);
  1734. * DEPLA (DOMA moddarcy 'MAILLAGE') AFFI (15.) (0. 0.) (1. 0.) ;
  1735. * 'FINSI' ;
  1736.  
  1737.  
  1738. * le nouveau ANT2 caracterise la précipitation dissolution
  1739.  
  1740. zozo = ANT2 '+' ANT3 ;
  1741. 'DETRUIT' ANT3 ;
  1742. ANT4 = 'MASQ' zozo 'SUPERIEUR' 0.1 ;
  1743. 'DETRUIT' zozo ;
  1744.  
  1745. zozo = ANT4 '-' ANT0 ;
  1746. CRIT = 'ABS' zozo ;
  1747. 'DETRUIT' zozo ;
  1748.  
  1749. 'DETRUIT' ANT0 ;
  1750. ANT0 = 'COPIER' ANT4 ;
  1751. 'DETRUIT' ANT4 ;
  1752.  
  1753. MAXCRIT = 'MAXIMUM' CRIT ;
  1754. 'DETRUIT' CRIT ;
  1755.  
  1756. 'SI' (MAXCRIT < 0.1) ;
  1757. * On a convergé dans le point fixe.
  1758. * On stocke la nouvelle distribution du précipité, ANT0
  1759. * est déjà à jour (prec et dissolution)
  1760. 'DETRUIT' ANTP ;
  1761. ANTP = 'COPIER' ANT2 ;
  1762. 'DETRUIT' ANT2 ;
  1763. 'DETRUIT' PRCIPF ;
  1764. 'DETRUIT' SRCPRCI ;
  1765. 'QUITTER' BOU8 ;
  1766. 'SINON' ;
  1767.  
  1768. * On résout par point fixe
  1769.  
  1770. * La où precipitation apparait, on ne fait rien de
  1771. * special. C'est pris en compte dans le nouveau ant0
  1772.  
  1773. * La 'OU' le precipité a diparu, le fait qu'il devienne
  1774. * négatif est génant. ON dissout en fait le précipité
  1775. * entierement au début du nouvel itéré de point fixe.
  1776. * Dans la mesure où il est censé disparaître, on espere
  1777. * que la concentration finale sera en dessous de la
  1778. * saturation. Sinon plusieurs itérés auront lieu avant
  1779. * convergence. Par contre on ne veut SURTOUT pas faire
  1780. * précipiter cet exces temporaire de concentration. DOnc
  1781. * on met valeur de ant0 indiquant qu'il n'a 'PLUS'
  1782. * de précipité ('ET' pas de précipitation 'NON' 'PLUS')
  1783.  
  1784. * 1 là 'OU' le précipité vient de diparaitre par
  1785. * rapport à l'état initial du point fixe !!!!
  1786. ANTDMM = ANTP '-' ANT2 ;
  1787. * on teste les valeurs positives correspondant à une disparition
  1788. * effective du precipité (le cas écheant c'est une apparition)
  1789. ANTDUM = 'MASQUE' ANTDMM SUPERIEUR 0.1 ;
  1790. 'DETRUIT' ANTDMM ;
  1791. * 0 la 'OU' le precipité vient de disparaitre, 1 ailleurs
  1792. ANTDM1 = 1.D0 '-' ANTDUM ;
  1793. * On s'assure d'interdire toute précipitation là 'OU' le
  1794. * précipité à disparu 'ET' l'on relache son état initial
  1795. * sous forme d'exces de concentration dans le point fixe.
  1796. zozo = ANTDM1 * ANT0 ;
  1797. 'DETRUIT' ANT0 ;
  1798. ANT0 = 'COPIER' zozo ;
  1799. 'DETRUIT' zozo ;
  1800. * relâchement de l'état initial du précipité par volume
  1801. * quantité (1 '-' poros) prec * volum '/' DELTAT
  1802. * on retire la patie explicite de décroissance du précipité
  1803. * SRCPRCI = coff
  1804. * * ((1.D0 '-' POROS)) * ANTDUM * PRCI '/' DELTAT ;
  1805. coff = (1.D0 '-' (DELTAT * LAMBD0 * (1.D0 '-' BETA))) ;
  1806. zozo = 1.D0 '-' POROS ;
  1807. zizi = zozo * ANTDUM ;
  1808. 'DETRUIT' zozo ;
  1809. zozo = zizi * PRCI ;
  1810. 'DETRUIT' zizi ;
  1811. SRCPRCI = (coff '/' DELTAT) * zozo ;
  1812. 'DETRUIT' zozo ;
  1813. * On retire cela du precipité initial
  1814. 'DETRUIT' PRCIPF ;
  1815. PRCIPF = ANTDM1 * PRCI ;
  1816.  
  1817. 'DETRUIT' ANTDUM ;
  1818. 'DETRUIT' ANTDM1 ;
  1819. 'DETRUIT' ANT2 ;
  1820.  
  1821. * Coefficient correspondant à la décroissance implicite :
  1822. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1823. * Terme implicite de dissolution, schéma Euler
  1824. DEN1 = CODIS * ANT0 * DELTAT ;
  1825. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1826. * solub en poros*codis Id et non Retard*poros*codis
  1827. zozo = DEN0 * PORETSU ;
  1828. zeze = DEN1 * POROS ;
  1829. 'DETRUIT' COFDC ;
  1830. COFDC = zozo '+' zeze ;
  1831. 'DETRUIT' zozo ;
  1832. 'DETRUIT' zeze ;
  1833. * On charge le coef de DC/DT dans GEOL1 pour TRANGEOL
  1834. GEOL1 . 'POROSITE' = COFDC ;
  1835. * On indque que le coef devant DC/DT a changé
  1836. TABMODI . 'POROSITE' = VRAI ;
  1837. NBIT = NBIT + 1 ;
  1838. *
  1839. 'FINSI' ;
  1840. 'SI' (&BOU8 EGA ITMAXI) ;
  1841. 'MESS' 'Sortie sans convergence' ;
  1842. 'QUITTER' BOU8 ;
  1843. 'FINSI' ;
  1844.  
  1845.  
  1846.  
  1847. *-----------
  1848. 'FIN' BOU8 ;
  1849. *-----------
  1850.  
  1851.  
  1852.  
  1853. NB0 = &BOU8 - 2 ;
  1854. 'SI' (NB0 >EG 1) ;
  1855. 'MESS' ' ' NB0 ' itérations supp de solubilité' ;
  1856. 'FINSI' ;
  1857.  
  1858. * 'Fin module Dissolution d ordre 1'
  1859. 'FINSI' ;
  1860.  
  1861. *|-------------------------------------|
  1862. *| Retard Non Linéaire |
  1863. *|-------------------------------------|
  1864. *- Calcul du nouveau coefficient de retard et de l'adsorbât
  1865. 'SI' RNONLIN ;
  1866. * 1. F(C)
  1867. CC0 = 'NOMC' CHRG 'SCAL' ;
  1868. 'SI' RNONLINL ;
  1869. * Langmuir
  1870. FF = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 )) ;
  1871. 'SINON' ;
  1872. * Freundlich
  1873. * La bidouille utilisant les masques sert à obtenir le champ
  1874. * des signes de CC0.
  1875. a = 'MASQUE' CC0 'SUPERIEUR' 0. ;
  1876. b = 'MASQUE' CC0 'INFERIEUR' 0. ;
  1877. FF = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b) ;
  1878. 'DETRUIT' a ;
  1879. 'DETRUIT' b ;
  1880. 'FINSI' ;
  1881. FF0 = FF ;
  1882. * 2. F(C + dC)
  1883. CC1 = ('NOMC' CHA2 'SCAL') + EpsCord ;
  1884. 'SI' RNONLINL ;
  1885. * Langmuir
  1886. FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 )) ;
  1887. 'SINON' ;
  1888. * Freundlich
  1889. a = 'MASQUE' CC1 'SUPERIEUR' 0. ;
  1890. b = 'MASQUE' CC1 'INFERIEUR' 0. ;
  1891. FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b) ;
  1892. 'DETRUIT' a ;
  1893. 'DETRUIT' b ;
  1894. 'FINSI' ;
  1895. * 3. R = 1 + dF/dC
  1896. RETAR2 = 1.D0 + ( (FF1 - FF0) / (CC1 - CC0) ) ;
  1897. 'DETRUIT' CC1 ;
  1898. 'DETRUIT' FF1 ;
  1899. *
  1900. *- Critère de sortie de la boucle de retard non linéaire :
  1901. CRIT = ( 'ABS' (RETAR2 - RETARC) ) / RETARC ;
  1902. CRIT = 'MAXIMUM' CRIT ;
  1903. 'SI' (CRIT < EPSRNL) ;
  1904. 'QUITTER' BOURNL ;
  1905. 'SINON' ;
  1906. * Le coef. de retard ne change que si le critère est mauvais.
  1907. * Réévaluation des variables affectées devant DC/DT:
  1908. RETARC = RETAR2 ;
  1909. PORETSU = RETARC * POROS ;
  1910. * coef devant dc/dt modifié
  1911. TABMODI . 'POROSITE' = VRAI ;
  1912.  
  1913. * Fonctions de volume réactualisées
  1914. PORETSU = RETARC * POROS ;
  1915. DECROI0 = 'EXCO' ((-1.) * LAMBD0 * PORETSU ) 'SCAL' 'SCAL';
  1916. DECROI0 = VOLU1 * DECROI0 ;
  1917.  
  1918. 'FINSI' ;
  1919. 'SINON' ;
  1920. * NE sert a rien GBM
  1921. * Retard linéaire F = (R - 1) C
  1922. * zozo = 'NOMC' CHRG 'SCAL' ;
  1923. * FF = RETAR1M1 * zozo ;
  1924. * 'DETRUIT' zozo ;
  1925. 'FINSI' ;
  1926.  
  1927. * Si on atteint le nombre maximum d'itérations - en particulier si le
  1928. * retard est linéaire - on quitte la boucle :
  1929. 'SI' (IRNL 'EGA' ITMAXRNL) ;
  1930. ' sortie boucle de retard non linéaire sans convergence' ;
  1931. 'QUITTER' BOURNL ;
  1932. 'FINSI' ;
  1933.  
  1934. *-------------
  1935. 'FIN' BOURNL ;
  1936. *-------------
  1937.  
  1938. 'SI' (IRNL > 1 ) ;
  1939. 'MESS' ' ' IRNL ' itérations de retard non linéaire.'
  1940. ' Critère = ' crit ;
  1941. 'FINSI' ;
  1942. NBITRNL = NBITRNL + (IRNL - 1) ;
  1943.  
  1944.  
  1945. *|-----------------------------------|
  1946. *| Fin des modules de résolution |
  1947. *|-----------------------------------|
  1948. *
  1949. *- Calcul de la dissolution par unité de volume et par unité de temps :
  1950. *
  1951. 'SI' SOLUA ;
  1952. DIS2 = DIS1 / VOLU1 / DELTAT ;
  1953. 'FINSI' ;
  1954. *
  1955. *- Archivage des resultats
  1956. *
  1957. * si le prochain temps à sauver tombe entre les deux
  1958. * incréments de temps, on obtient les valeurs à sauver
  1959. * par interpolation linéaire puis on les stocke.
  1960. *
  1961.  
  1962. *list ((PLACE) * 4. / 1000000.) ;
  1963.  
  1964. 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ;
  1965.  
  1966. * Sauvegarde de tous les temps intermédiaires, s'il y en a :
  1967. 'REPETER' BOU9 ;
  1968. TEMS = 'EXTR' TPSOR ISOR ;
  1969.  
  1970. 'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '&lt;EG' TPS ) ) ;
  1971. * on interpole les variables linéairement entre
  1972. * les pas de temps
  1973. LAST1 = LAST1 + 1 ;
  1974. DTEM = ( TPS - TEMS ) / DELTAT ;
  1975. UNMO = 1.D0 - DTEM ;
  1976. CHBIS = 'COLI' CHRG DTEM CHA2 UNMO ;
  1977. FLUBIS = 'COLI' FLU0 DTEM FLU2 UNMO ;
  1978. FLUBISCO= 'COLI' FLUCO0 DTEM FLUCO2 UNMO ;
  1979. 'SI' ('EGA' ('TYPE' RETAR0) 'CHPOINT ') ;
  1980. RETBIS = 'COLI' RETAR0 DTEM RETARC UNMO ;
  1981. 'SINON' ;
  1982. RETBIS = (DTEM * RETAR0) + (UNMO * RETARC) ;
  1983. 'FINSI' ;
  1984. 'SI' SOLUB ;
  1985. PRBIS = 'COLI' PRCI DTEM PRE2 UNMO ;
  1986. 'FINSI' ;
  1987. TRANSI . 'TEMPS' . LAST1 = TEMS ;
  1988. TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP
  1989. CHBIS ;
  1990. 'DETRUIT' CHBIS ;
  1991. TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP
  1992. FLUBIS ;
  1993. 'DETRUIT' FLUBIS ;
  1994. TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP
  1995. FLUBISCO ;
  1996. 'DETRUIT' FLUBISCO ;
  1997. 'SI' SOLUB ;
  1998. TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP
  1999. PRBIS ;
  2000. 'DETRUIT' PRBIS ;
  2001. 'FINSI' ;
  2002. 'SI' SAUVRET ;
  2003. TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP
  2004. RETBIS ;
  2005. 'DETRUIT' RETBIS ;
  2006. 'FINSI' ;
  2007. ISOR = ISOR + 1 ;
  2008. 'SINON' ;
  2009. 'QUITTER' BOU9 ;
  2010. 'FINSI' ;
  2011. 'SI' ( ISOR '>' DSOR ) ;
  2012. ISOR = DSOR ;
  2013. 'QUITTER' BOU9 ;
  2014. 'FINSI' ;
  2015. 'FIN' BOU9 ;
  2016. 'SINON' ;
  2017. * sinon, on sauvegarde tous les temps :
  2018. LAST1 = LAST1 + 1 ;
  2019. TRANSI . 'TEMPS' . LAST1 = TPS ;
  2020. TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP CHA2 ;
  2021. TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP FLU2 ;
  2022. TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP FLUCO2 ;
  2023. 'SI' SOLUB ;
  2024. TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP PRE2 ;
  2025. 'FINSI' ;
  2026. 'SI' SAUVRET ;
  2027. TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP RETARC ;
  2028. 'FINSI' ;
  2029. 'FINSI' ;
  2030.  
  2031.  
  2032. *On sauve l'intégrale de decroissance*concentration et
  2033. * decroissance*précipité au cours du temps
  2034. * stockée pour l'instant à tous les temps. Peut s'améliorer.
  2035.  
  2036. * indice pour le stockage des intégrales de concentration
  2037. LASTINT = LASTINT '+' 1;
  2038.  
  2039. BM1 = 1.D0 '-' BETA ;
  2040. * si l'espece fait l'objet d'une filiation on stocke les intégrales
  2041. * du précipité et de la concentration.
  2042. SI ((EXISTE TRANSI DECROISSANCE) 'ET' (TRANSI . 'PERE')) ;
  2043. *TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1)
  2044. * '+' (DELTAT * PORETSU * TRANSI . 'DECROISSANCE' *
  2045. * ((BM1 * ('NOMC' NOMESP CHRG))
  2046. * '+' ((BETA) * ('NOMC' NOMESP CHA2)))) ;
  2047. zozo = (DELTAT * TRANSI . 'DECROISSANCE') * PORETSU ;
  2048. zaza = BM1 * CHRG ;
  2049. zuzu = BETA * CHA2 ;
  2050. zeze = 'KOPS' zaza '+' zuzu ;
  2051. 'DETRUIT' zaza ;
  2052. 'DETRUIT' zuzu ;
  2053. zaza = 'KOPS' zozo * zeze ;
  2054. 'DETRUIT' zozo ;
  2055. 'DETRUIT' zeze ;
  2056. zozo = 'NOMC' NOMESP zaza ;
  2057. 'DETRUIT' zaza ;
  2058. TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1)
  2059. '+' zozo ;
  2060. 'DETRUIT' zozo ;
  2061.  
  2062. 'SI' SOLUB ;
  2063.  
  2064. * TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . LASTINT
  2065. * '+' (DELTAT * (1.D0 '-' POROS) * TRANSI . 'DECROISSANCE' *
  2066. * ((BM1 * ('NOMC' NOMESP PRCI))
  2067. * '+' ((BETA) * ('NOMC' NOMESP PRE2)))) ;
  2068. zizi = (1.D0 '-' POROS) ;
  2069. zozo = (DELTAT * TRANSI . 'DECROISSANCE') * zizi ;
  2070. 'DETRUIT' zizi ;
  2071. zaza = BM1 * PRCI ;
  2072. zuzu = BETA * PRE2 ;
  2073. zeze = 'KOPS' zaza '+' zuzu ;
  2074. 'DETRUIT' zaza ;
  2075. 'DETRUIT' zuzu ;
  2076. zaza = 'KOPS' zozo * zeze ;
  2077. 'DETRUIT' zozo ;
  2078. 'DETRUIT' zeze ;
  2079. zozo = 'NOMC' NOMESP zaza ;
  2080. 'DETRUIT' zaza ;
  2081. zaza = TRANSI . 'INTFILI' . LASTINT '+' zozo ;
  2082. 'DETRUIT' TRANSI . 'INTFILI' . LASTINT ;
  2083. 'DETRUIT' zozo ;
  2084. TRANSI . 'INTFILI' . LASTINT = zaza ;
  2085. 'OUBLIER' zaza ;
  2086.  
  2087. 'FINSI' ;
  2088. FINSI ;
  2089.  
  2090.  
  2091. * Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait.
  2092. * Le dernier pas de temps ne peut pas faire l'objet d'interpolation
  2093. * linéaire si le temps sauvegardé ne tombe pas pile dessus.
  2094. 'SI' (IPAS 'EGA' DCAL) ;
  2095. TPSDER = TRANSI. 'TEMPS' . LAST1 ;
  2096. 'SI' ( TPSDER 'NEG' TPSFIN EPSDT ) ;
  2097. LAST1 = LAST1 + 1 ;
  2098. TRANSI . 'TEMPS' . LAST1 = TPS ;
  2099. TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP CHA2 ;
  2100. TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP FLU2 ;
  2101. TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP FLUCO2 ;
  2102. 'SI' SOLUB ;
  2103. TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP PRE2 ;
  2104. 'FINSI' ;
  2105. 'SI' SAUVRET ;
  2106. TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP RETARC ;
  2107. 'FINSI' ;
  2108. 'FINSI' ;
  2109. 'FINSI' ;
  2110. *
  2111. *- Initialisations pour le pas suivant
  2112. *
  2113. PRECED = TPS ;
  2114. 'DETRUIT' CHRG ;
  2115. CHRG = CHA2 ;
  2116. 'DETRUIT' FLU0 ;
  2117. FLU0 = FLU2 ;
  2118. 'DETRUIT' FLUCO0 ;
  2119. FLUCO0 = FLUCO2 ;
  2120. 'OUBLIER' GEOL1 ;
  2121. 'OUBLIER' GEOL2 ;
  2122. GEOL1 = 'TABLE' GEOLPF1 ;
  2123. GEOL2 = 'TABLE' GEOLPF2 ;
  2124. 'DETRUIT' RETAR0 ;
  2125. RETAR0 = RETARC ;
  2126. DELOLD = DELTAT ;
  2127. 'SI' SOLUB ;
  2128. 'DETRUIT' PRCI ;
  2129. PRCI = PRE2 ;
  2130. 'FINSI' ;
  2131. 'SI' ('EXISTE' TRANSI 'SOURCE') ;
  2132. 'DETRUIT' TERSC2M1 ;
  2133. TERSC2M1 = TERSC2 ;
  2134. 'FINSI' ;
  2135. 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE') ;
  2136. 'DETRUIT' FLUIMPM1 ;
  2137. FLUIMPM1 = FLUIMPO ;
  2138. 'FINSI' ;
  2139. 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP') ;
  2140. 'DETRUIT' FLUTMPM1 ;
  2141. FLUTMPM1 = FLUTMPO ;
  2142. 'FINSI' ;
  2143. 'SI' ('EXISTE' TRANSI 'MIXTES') ;
  2144. 'DETRUIT' FLUMMPM1 ;
  2145. FLUMMPM1 = FLUMMPO ;
  2146. 'FINSI' ;
  2147.  
  2148. *MESS 'TEMPS2' ; TEMPS;
  2149. *'MENAGE' ;
  2150.  
  2151. * destruction d'objets temporaires
  2152.  
  2153. 'DETRUIT' TERSC3 ;
  2154.  
  2155. *list ((PLACE) * 4. / 1000000.) ;
  2156.  
  2157. *temps impr;
  2158.  
  2159. *
  2160. *=============
  2161. 'FIN' BOUTPS ; COMM 'Boucle sur le temps' ;
  2162. *=============
  2163. *
  2164. *- Indication du nombre de passages :
  2165. *
  2166. 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION')) ;
  2167. 'MESS' 'Nombre total de matrices calculées : ' NBIT ;
  2168. 'FINSI' ;
  2169. 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION')) ;
  2170. 'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
  2171. NBIT ;
  2172. 'FINSI' ;
  2173. 'SI' SOLUP ;
  2174. 'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
  2175. NBIT;
  2176. 'FINSI' ;
  2177. 'SI' RNONLIN ;
  2178. 'MESS' 'Nombre total d itérations supplémentaires de '
  2179. 'retard non linéaire : ' NBITRNL;
  2180. 'FINSI' ;
  2181.  
  2182. 'FINP' ;
  2183.  
  2184.  
  2185.  
  2186.  

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