Télécharger darcytra.procedur

Retour à la liste

Numérotation des lignes :

  1. * DARCYTRA PROCEDUR JC220346 14/02/19 21:15:02 7941
  2. 'DEBP' DARCYTRA TRANSI*'TABLE';
  3.  
  4. *---------------------------------------------------------------------
  5. * Résolution de l'équation de darcy pour un problème d'écoulement ou
  6. * de transport par une méthode d'éléments
  7. * finis mixtes hybrides. Les inconnues du problème sont
  8. * - pour l'écoulement, la charge (H) et trace de charge (TH) et le
  9. * débit hydraulique (FLUX)
  10. * - pour le transport, la concentration (H), la trace de
  11. * concentration (TH) et le débit diffusif (FLUX).
  12. *---------------------------------------------------------------------
  13. *
  14. *------------------------------
  15. * Phrase d'appel (en GIBIANE) :
  16. *------------------------------
  17. *
  18. * DARCYTRA TABLE ;
  19. *
  20. *----------------------------------
  21. * Opérandes (à mettre dans TABLE) :
  22. *----------------------------------
  23. *
  24. * ___________________________________________________________________
  25. * | |
  26. * | Indice Contenu |
  27. * | |
  28. * -------------------------------------------------------------------
  29. * | |
  30. * |------------------------------------------------ |
  31. * |Données physiques, géométriques et materielles : |
  32. * |------------------------------------------------ |
  33. * | |
  34. * | ------ Indices communs à l'écoulement et au transport ------ |
  35. * | ------------------------------------------------------------ |
  36. * | |
  37. * |'SOUSTYPE' 'DARCY' (type MOT) |
  38. * | |
  39. * |'MODELE' Objet modèle (MMODEL créé par MODE) |
  40. * | |
  41. * | |
  42. * | ------ 1ère possibilité : Résolution de l'écoulement ------ |
  43. * | ----------------------------------------------------------- |
  44. * | |
  45. * |'CARACTERISTIQUES' Données physiques et materielles : |
  46. * | conductivité hydraulique (CHAMELEM créé par MATE) |
  47. * | |
  48. * |'EMMAGASINEMENT' Valeur du coefficient d'emmagasinement |
  49. * | (Type CHPO Centre, Comp 'CK', ou FLOTTANT) |
  50. * | - Défaut 1. |
  51. * | |
  52. * | ------ 2ème possibilité : Résolution du transport ------ |
  53. * | ------------------------------------------------------------ |
  54. * | |
  55. * |'CARACTERISTIQUES' Données physiques et materielles : |
  56. * | diffusivité effective (CHAMELEM créé par MATE) |
  57. * | |
  58. * |'POROSITE' Valeur de la porosité (Type CHAMELEM, Comp |
  59. * | 'CK', ou FLOTTANT) - Défaut 1. |
  60. * | |
  61. * |'DECROISSANCE' Valeur du terme de décroissance (Type FLOTTANT) |
  62. * | Tel que dC/dt = - Lambda * C - Défaut 0. |
  63. * | |
  64. * |'COEF_RETARD' Coefficient de retard linéaire dans le cas simple, |
  65. * | ou Pente à l'origine de la fonction F(C) dans le |
  66. * | cas d'isotherme non linéaire de Langmuir |
  67. * | ou Coefficient K de l'isotherme de Freundlich |
  68. * | (Type CHPO Centre 'SCAL', ou FLOTTANT) |
  69. * | |
  70. * |'LANGMUIR' Quantité maximale 'Fsat' adsorbée sur le solide |
  71. * | rapportée à l'unité de volume du fluide et exprimée|
  72. * | dans la meme unité que le soluté. |
  73. * | (Type CHPO Centre 'SCAL', ou FLOTTANT). |
  74. * | F = (R-1) C / [1 + ((R-1) C / Fsat)] |
  75. * | Si cet indice et le suivant sont absents, |
  76. * | l'équilibre d'adsorption est linéaire. Cet indice a|
  77. * | priorité sur l'indice FREUNDLICH. |
  78. * | |
  79. * |'FREUNDLICH' Exposant de la loi de Freundlich F = K (C ^ 1/n) |
  80. * | (Type FLOTTANT). |
  81. * | Dans ce cas (et si l'indice LANGMUIR n'existe pas),|
  82. * | l'indice 'COEF_RETARD' contient le coefficient |
  83. * | K ramené à une unité de volume de fluide. |
  84. * | - Non disponible pour l'instant - |
  85. * | |
  86. * |'LIMITE_SOLUBILITE' Limite de solubilité (Type MCHAML, Comp 'H') |
  87. * | |
  88. * |'COEF_DISSOLUTION' Coef. de dissolution (Type CHPO Centre, Comp |
  89. * | 'SCAL'). Tel que dC/dt = Coef * (Csat - C) - Par |
  90. * | défaut, la dissolution est instantanée |
  91. * | |
  92. * |'CONVECTION' Débit intégré de la vitesse convective à travers |
  93. * | chaque face des éléments (Type CHPO Face, comp. |
  94. * | 'FLUX') |
  95. * | |
  96. * |---------------------- |
  97. * |Conditions initiales : |
  98. * |---------------------- |
  99. * | |
  100. * | ------ Indices communs à l'écoulement et au transport ------ |
  101. * | ------------------------------------------------------------ |
  102. * | |
  103. * |'TEMPS' TABLE contenant à l'indice 0 la valeur du temps |
  104. * | initial (FLOTTANT) |
  105. * | |
  106. * | |
  107. * | ------ 1ère possibilité : Résolution de l'écoulement ------ |
  108. * | ----------------------------------------------------------- |
  109. * | |
  110. * |'CHARGE' TABLE contenant à l'indice 0 la charge hydraulique |
  111. * | (quantité d'élément par unité de volume d'eau) |
  112. * | (Type CHPO Centre, Comp 'H') |
  113. * | |
  114. * |'TRACE_CHARGE' TABLE contenant à l'indice 0 la trace de |
  115. * | charge initiale (CHPO, 'TH') |
  116. * | |
  117. * |'FLUX' TABLE contenant à l'indice 0 le flux hydraulique |
  118. * | initial intégré sur chaque face (Type CHPO Face, |
  119. * | comp. 'FLUX') |
  120. * | |
  121. * | ------ 2ème possibilité : Résolution du transport ------ |
  122. * | ------------------------------------------------------------ |
  123. * | |
  124. * |'CONCENTRATION' TABLE contenant à l'indice 0 la concentration |
  125. * | (quantité d'élément par unité de volume d'eau) |
  126. * | (Type CHPO Centre, Comp 'H') |
  127. * | |
  128. * |'TRACE_CONC' TABLE contenant à l'indice 0 la trace de |
  129. * | concentration initiale (Type CHPO Face, Comp 'TH') |
  130. * | |
  131. * |'FLUX' TABLE contenant à l'indice 0 le flux diffusif |
  132. * | initial intégré sur chaque face (Type CHPO Face, |
  133. * | comp. 'FLUX') |
  134. * | |
  135. * |'PRECIPITE' TABLE contenant à l'indice 0 la quantité initiale |
  136. * | de précipité par unité de volume de milieu |
  137. * | (Type CHPO Centre, Comp 'H') |
  138. * | |
  139. * |'DISSOLUTION' TABLE contenant à l'indice 0 la quantité initiale |
  140. * | pour estimer la dissolution au premier pas de temps|
  141. * | (Type CHPO, Comp 'H'), voir plus loin. |
  142. * | |
  143. * |-------------------------------------- |
  144. * |Conditions aux limites / chargements : |
  145. * |-------------------------------------- |
  146. * | |
  147. * | ------ Indices communs à l'écoulement et au transport ------ |
  148. * | ------------------------------------------------------------ |
  149. * | |
  150. * |'BLOCAGE' Contient les matrices de blocage (RIGIDITE) |
  151. * | |
  152. * |'TRACE_IMPOSE' Valeurs des traces imposées (charge ou concentra- |
  153. * | -tion) (CHARGEMENT 'TH' - Obligatoire si BLOCAGE) |
  154. * | |
  155. * |'FLUX_IMPOSE' Valeurs des flux imposés intégrés par face |
  156. * | (Type CHARGEMENT de CHPO Face, comp. 'FLUX'- |
  157. * | défaut 0.) |
  158. * | |
  159. * |'SOURCE' Valeurs du terme source par maille et par unité de |
  160. * | temps (ex : puits, filiation) |
  161. * | Les valeurs à l'indice i sont les valeurs entre |
  162. * | les temps i-1 et i. |
  163. * | (CHARGEMENT de CHPO Centre, comp 'SOUR'- défaut 0.)|
  164. * | |
  165. * | ------ 2ème possibilité : Résolution du transport ------ |
  166. * | ------------------------------------------------------------ |
  167. * | |
  168. * |'DISSOLUTION_IMPOSEE' Valeurs des dissolutions imposées par unité|
  169. * | de temps et par maille. (Type CHARGEMENT de CHPO, |
  170. * | Comp 'H'). Les valeurs à l'indice i sont les |
  171. * | valeurs moyennes de dissolution par unité de temps |
  172. * | entre les temps i-1 et i. |
  173. * | Priorité de la dissolution imposée sur les |
  174. * | cinétiques. |
  175. * | |
  176. * |-------------------- |
  177. * |Données numériques : |
  178. * |-------------------- |
  179. * | |
  180. * | ------ Indices communs à l'écoulement et au transport ------ |
  181. * | ------------------------------------------------------------ |
  182. * | |
  183. * |'TEMPS_CALCULES' Valeur des temps calculés (LISTREEL) |
  184. * | Contient obligatoirement le temps final. |
  185. * | |
  186. * |'TEMPS_SAUVES' Valeur des temps sauvegardés (LISTREEL - défaut : |
  187. * | on sauve tous les pas de temps) |
  188. * | |
  189. * | ------ 1ère possibilité : Résolution de l'écoulement ------ |
  190. * | ----------------------------------------------------------- |
  191. * | |
  192. * |'THETA' Coefficient de relaxation compris entre 0. et 1. |
  193. * | (theta-méthode diffusion) (FLOTTANT - défaut 1.) |
  194. * | Possibilité de non-convergence lorsque theta<1/2 |
  195. * | Valeurs de theta généralement utilisées : |
  196. * | Schéma de Euler explicite : 0. |
  197. * | Schéma de Crank-Nicholson : 1/2 |
  198. * | Schéma de Galerkin : 2/3 |
  199. * | Schéma de Euler implicite : 1. |
  200. * | |
  201. * | ------ 2ème possibilité : Résolution du transport ------ |
  202. * | ------------------------------------------------------------ |
  203. * | |
  204. * |'THETA_DIFF' Coefficient de relaxation compris entre 0. et 1. |
  205. * | (theta-méthode diffusion) ('FLOTTANT' - défaut 1.) |
  206. * | |
  207. * |'THETA_CONVECTION' Idem pour la convection |
  208. * | ('FLOTTANT', Défaut = THETA_DIFF) |
  209. * | |
  210. * |'THETA_DEC' Idem mais pour la décroissance |
  211. * | ('FLOTTANT' - défaut 1/2) |
  212. * | |
  213. * |'THETA_DISS' Idem mais pour la dissolution |
  214. * | ('FLOTTANT' - défaut 1.) |
  215. * | |
  216. * |'PENALISATION' Coefficient de pénalisation pour la prise en |
  217. * | compte de la limite de solubilité. La présence de |
  218. * | cet indice ou du suivant indique quel schéma a été |
  219. * | choisi. |
  220. * | (Type 'FLOTTANT') - Valeur conseillée 1.D7 |
  221. * | |
  222. * |'EPSI_LIM' Précision relative d'arrêt pour le shéma limite de |
  223. * | solubilité prédicteur-correcteur itératif |
  224. * | (Type FLOTTANT) - Valeur conseillée 5.D-3 |
  225. * | |
  226. * |'ITMAX_LIM' Nombre maxi d'itérations correspondant aux modules |
  227. * | de dissolution avant d'abandonner |
  228. * | (Type 'ENTIER') - Défaut 50 |
  229. * | |
  230. * |'EPSI_RET' Précision relative d'arrêt pour la résolution |
  231. * | itérative (Picard) de l'adsorption non linéaire |
  232. * | (Type FLOTTANT) - Défaut 1.D-4 |
  233. * | |
  234. * |'EPSI_COR' Petit saut de concentration pour calculer le coef. |
  235. * | de retard par la méthode de la corde lorsque le |
  236. * | retard est non-linéaire. |
  237. * | (Type FLOTTANT) - Défaut 1.D-4 |
  238. * | |
  239. * |'ITMAX_RET' Nombre maxi d'itérations correspondant au retard |
  240. * | non linéaire avant d'abandonner. |
  241. * | (Type 'ENTIER') - Défaut 20 |
  242. * |_________________________________________________________________|
  243. *
  244. *
  245. *
  246. *---------------------------------
  247. * Résultats (stockés dans TABLE) :
  248. *---------------------------------
  249. *
  250. * ___________________________________________________________________
  251. * | |
  252. * | Indice Contenu |
  253. * | |
  254. * -------------------------------------------------------------------
  255. * | |
  256. * | ------ Indices communs à l'écoulement et au transport ------ |
  257. * | ------------------------------------------------------------ |
  258. * | |
  259. * |'TEMPS' TABLE contenant les temps sauvegardés (FLOTTANT) |
  260. * | |
  261. * | ------ 1ère possibilité : Résolution de l'écoulement ------ |
  262. * | ------------------------------------------------------------ |
  263. * | |
  264. * |'CHARGE' TABLE contenant les charges |
  265. * | (Type CHPO Centre, Comp 'H') |
  266. * | |
  267. * |'TRACE_CHARGE' TABLE contenant les traces de charge |
  268. * | (Type CHPO Face, Comp 'TH') |
  269. * | |
  270. * |'FLUX' TABLE contenant les débits hydrauliques intégrés |
  271. * | par face (Type CHPO Face, comp. 'FLUX') |
  272. * | |
  273. * | ------ 2ème possibilité : Résolution du transport ------ |
  274. * | ------------------------------------------------------------ |
  275. * | |
  276. * |'CONCENTRATION' TABLE contenant les concentrations |
  277. * | (Type CHPO Centre, Comp 'H') |
  278. * | |
  279. * |'TRACE_CONC' TABLE contenant les traces de concentration |
  280. * | (Type CHPO Face, Comp 'TH') |
  281. * | |
  282. * |'FLUX' TABLE contenant les débits diffusifs intégrés |
  283. * | par face (Type CHPO Face, comp. 'FLUX') |
  284. * | |
  285. * |'PRECIPITE' TABLE contenant la quantité de précipité par maille|
  286. * | (Type CHPO Centre, Comp 'H') |
  287. * | |
  288. * |'DISSOLUTION' TABLE contenant la quantité de précipité dissoute |
  289. * | entre deux pas de temps par unité de volume et par |
  290. * | unité de temps. La valeur stockée à l'indice i, |
  291. * | est valable entre les temps i-1 et i |
  292. * | (Type CHPO, Comp 'H'). |
  293. * | ATTENTION, les valeurs de cette table résultat |
  294. * | n'ont aucun sens lorsque les temps sauvegardés ne |
  295. * | sont pas les memes que les temps calculés. Toute |
  296. * | tentative d'exploitation donnera alors des |
  297. * | résultats incohérents (erreurs de bilan) |
  298. * | |
  299. * |'RETARD' Si cet indice a été préalablement défini comme une |
  300. * | TABLE, alors il contient les valeurs du coefficient|
  301. * | de retard (Type 'CHPO' centre, Comp 'SCAL'). Sinon,|
  302. * | les valeurs du coefficient de retard ne sont pas |
  303. * | sauvegardées. |
  304. * |_________________________________________________________________|
  305. *
  306. *
  307. * ___________________________________________________________________
  308. * | |
  309. * | Les tables résultats sont indicées par des entiers variant de 0 |
  310. * | à N . |
  311. * | A l'indice 0 on stocke les valeurs initiales, aux indices |
  312. * | suivants les champs correspondant au temps de sortie TEMPS.I . |
  313. * | Les champs servant en cas de reprise sont ceux correpondant au |
  314. * | dernier indice. |
  315. * |_________________________________________________________________|
  316. *
  317. *
  318. * SOUSTYPE
  319. 'SI' ( 'NON' ('EXISTE' TRANSI 'SOUSTYPE' ) ) ;
  320. 'ERREUR'
  321. 'Indice SOUSTYPE absent de la table de données.' ;
  322. 'QUITTER' DARCYTRA ;
  323. 'FINSI' ;
  324. 'SI' ( 'NEG' ( TRANSI.'SOUSTYPE' ) 'DARCY_TRANSITOIRE' ) ;
  325. 'ERREUR' 'Le sous-type de la table d entree est incorrect' ;
  326. 'QUITTER' DARCYTRA ;
  327. 'FINSI' ;
  328. * TYPE DE RESOLUTION
  329. * Le nom des variables suit la logique du transport.
  330. * Par la suite, si la résolution est hydro, la porosité tiendra lieu
  331. * de coef. d'emmagasinement, la concentration de charge (idem pour les
  332. * traces), et theta-diff de theta.
  333. * Pour chaque variable spécifique, on vérifie que sa présence
  334. * correspond bien au type de résolution envisagé.
  335. HYDRO = 'EXISTE' TRANSI 'CHARGE' ;
  336. TRANSP ='NON' HYDRO ;
  337. 'SI' HYDRO ;
  338. NOMINC = 'MOT' 'CHARGE' ;
  339. NOMTINC = 'MOT' 'TRACE_CHARGE' ;
  340. NOMTETA = 'MOT' 'THETA' ;
  341. NOMDDT = 'MOT' 'EMMAGASINEMENT' ;
  342. 'SI' ('EXISTE' TRANSI 'POROSITE') ;
  343. 'ERREUR' ('CHAINE' 'L indice POROSITE ne correspond '
  344. 'pas à la résolution hydraulique') ;
  345. 'FINSI' ;
  346. 'SI' ('EXISTE' TRANSI 'DECROISSANCE') ;
  347. 'ERREUR' ('CHAINE' 'L indice DECROISSANCE ne correspond '
  348. 'pas à la résolution hydraulique') ;
  349. 'FINSI' ;
  350. 'SI' ('EXISTE' TRANSI 'COEF_RETARD') ;
  351. 'ERREUR' ('CHAINE' 'L indice COEF_RETARD ne correspond '
  352. 'pas à la résolution hydraulique') ;
  353. 'FINSI' ;
  354. 'SI' ('EXISTE' TRANSI 'LANGMUIR') ;
  355. 'ERREUR' ('CHAINE' 'L indice LANGMUIR ne correspond '
  356. 'pas à la résolution hydraulique') ;
  357. 'FINSI' ;
  358. 'SI' ('EXISTE' TRANSI 'FREUNDLICH') ;
  359. 'ERREUR' ('CHAINE' 'L indice FREUNDLICH ne correspond '
  360. 'pas à la résolution hydraulique') ;
  361. 'FINSI' ;
  362. 'SI' ('EXISTE' TRANSI 'LIMITE_SOLUBILITE') ;
  363. 'ERREUR' ('CHAINE' 'L indice LIMITE_SOLUBILITE ne correspond '
  364. 'pas à la résolution hydraulique') ;
  365. 'FINSI' ;
  366. 'SI' ('EXISTE' TRANSI 'COEF_DISSOLUTION') ;
  367. 'ERREUR' ('CHAINE' 'L indice COEF_DISSOLUTION ne correspond '
  368. 'pas à la résolution hydraulique') ;
  369. 'FINSI' ;
  370. 'SI' ('EXISTE' TRANSI 'CONVECTION') ;
  371. 'ERREUR' ('CHAINE' 'L indice CONVECTION ne correspond '
  372. 'pas à la résolution hydraulique') ;
  373. 'FINSI' ;
  374. 'SI' ('EXISTE' TRANSI 'CONCENTRATION') ;
  375. 'ERREUR' ('CHAINE' 'L indice CONCENTRATION ne correspond '
  376. 'pas à la résolution hydraulique') ;
  377. 'FINSI' ;
  378. 'SI' ('EXISTE' TRANSI 'TRACE_CONC') ;
  379. 'ERREUR' ('CHAINE' 'L indice TRACE_CONC ne correspond '
  380. 'pas à la résolution hydraulique') ;
  381. 'FINSI' ;
  382. 'SI' ('EXISTE' TRANSI 'PRECIPITE') ;
  383. 'ERREUR' ('CHAINE' 'L indice PRECIPITE ne correspond '
  384. 'pas à la résolution hydraulique') ;
  385. 'FINSI' ;
  386. 'SI' ('EXISTE' TRANSI 'DISSOLUTION') ;
  387. 'ERREUR' ('CHAINE' 'L indice DISSOLUTION ne correspond '
  388. 'pas à la résolution hydraulique') ;
  389. 'FINSI' ;
  390. 'SI' ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE') ;
  391. 'ERREUR' ('CHAINE' 'L indice DISSOLUTION_IMPOSEE ne correspond '
  392. 'pas à la résolution hydraulique') ;
  393. 'FINSI' ;
  394. 'SI' ('EXISTE' TRANSI 'THETA_DIFF') ;
  395. 'ERREUR' ('CHAINE' 'L indice THETA_DIFF ne correspond '
  396. 'pas à la résolution hydraulique') ;
  397. 'FINSI' ;
  398. 'SI' ('EXISTE' TRANSI 'THETA_CONVECTION') ;
  399. 'ERREUR' ('CHAINE' 'L indice THETA_CONVECTION ne correspond '
  400. 'pas à la résolution hydraulique') ;
  401. 'FINSI' ;
  402. 'SI' ('EXISTE' TRANSI 'THETA_DEC') ;
  403. 'ERREUR' ('CHAINE' 'L indice THETA_DEC ne correspond '
  404. 'pas à la résolution hydraulique') ;
  405. 'FINSI' ;
  406. 'SI' ('EXISTE' TRANSI 'THETA_DISS') ;
  407. 'ERREUR' ('CHAINE' 'L indice THETA_DISS ne correspond '
  408. 'pas à la résolution hydraulique') ;
  409. 'FINSI' ;
  410. 'SI' ('EXISTE' TRANSI 'PENALISATION') ;
  411. 'ERREUR' ('CHAINE' 'L indice PENALISATION ne correspond '
  412. 'pas à la résolution hydraulique') ;
  413. 'FINSI' ;
  414. 'SI' ('EXISTE' TRANSI 'EPSI_LIM') ;
  415. 'ERREUR' ('CHAINE' 'L indice EPSI_LIM ne correspond '
  416. 'pas à la résolution hydraulique') ;
  417. 'FINSI' ;
  418. 'SI' ('EXISTE' TRANSI 'EPSI_RET') ;
  419. 'ERREUR' ('CHAINE' 'L indice EPSI_RET ne correspond '
  420. 'pas à la résolution hydraulique') ;
  421. 'FINSI' ;
  422. 'SI' ('EXISTE' TRANSI 'EPSI_COR') ;
  423. 'ERREUR' ('CHAINE' 'L indice EPSI_COR ne correspond '
  424. 'pas à la résolution hydraulique') ;
  425. 'FINSI' ;
  426. 'SI' ('EXISTE' TRANSI 'ITMAX_RET') ;
  427. 'ERREUR' ('CHAINE' 'L indice ITMAX_RET ne correspond '
  428. 'pas à la résolution hydraulique') ;
  429. 'FINSI' ;
  430. 'SI' ('EXISTE' TRANSI 'RETARD') ;
  431. 'ERREUR' ('CHAINE' 'L indice RETARD ne correspond '
  432. 'pas à la résolution hydraulique') ;
  433. 'FINSI' ;
  434. 'FINSI' ;
  435. 'SI' TRANSP ;
  436. NOMINC = 'MOT' 'CONCENTRATION' ;
  437. NOMTINC = 'MOT' 'TRACE_CONC' ;
  438. NOMTETA = 'MOT' 'THETA_DIFF' ;
  439. NOMDDT = 'MOT' 'POROSITE' ;
  440. 'SI' ('EXISTE' TRANSI 'THETA') ;
  441. 'ERREUR' ('CHAINE' 'L indice THETA ne correspond '
  442. 'pas à la résolution du transport') ;
  443. 'FINSI' ;
  444. 'SI' ('EXISTE' TRANSI 'TRACE_CHARGE') ;
  445. 'ERREUR' ('CHAINE' 'L indice TRACE_CHARGE ne correspond '
  446. 'pas à la résolution du transport') ;
  447. 'FINSI' ;
  448. 'SI' ('EXISTE' TRANSI 'EMMAGASINEMENT') ;
  449. 'ERREUR' ('CHAINE' 'L indice EMMAGASINEMENT ne correspond '
  450. 'pas à la résolution du transport') ;
  451. 'FINSI' ;
  452. 'FINSI' ;
  453. *
  454. *--------------------------------------------------------------------*
  455. * RECUPERATION DES DONNEES PHYSIQUES, GEOMETRIQUES ET MATERIELLES *
  456. *--------------------------------------------------------------------*
  457. *
  458. * MODELE
  459. 'SI' ( 'EXISTE' TRANSI 'MODELE' ) ;
  460. MODHYB = TRANSI . 'MODELE' ;
  461. HYTOT = DOMA MODHYB 'TABLE' ;
  462. 'SINON' ;
  463. 'ERREUR' 'Il manque le modele.' ;
  464. 'QUITTER' DARCYTRA ;
  465. 'FINSI' ;
  466. * VOLUME
  467. VOLU1 = 'DOMA' MODHYB 'VOLUME' ;
  468. * ORIENTATION
  469. MCHYB = 'DOMA' MODHYB 'ORIENTAT' ;
  470. * CARACTERISTIQUES
  471. 'SI' ( 'EXISTE' TRANSI 'CARACTERISTIQUES' ) ;
  472. MAT1 = TRANSI . 'CARACTERISTIQUES' ;
  473. 'SINON' ;
  474. 'ERREUR'
  475. 'Indice CARACTERISTIQUES absent de la table de données.' ;
  476. 'QUITTER' DARCYTRA ;
  477. 'FINSI' ;
  478. * EMMAGASINEMENT OU POROSITE
  479. 'SI' ( 'EXISTE' TRANSI NOMDDT ) ;
  480. POROS = TRANSI . NOMDDT ;
  481. 'SINON' ;
  482. POROS = 1.D0 ;
  483. 'FINSI' ;
  484. 'SI' ('EGA' ('TYPE' POROS) 'FLOTTANT') ;
  485. POROS = 'MANU' 'CHML' HYTOT.'MAILLAGE' 'CK' POROS ;
  486. 'FINSI' ;
  487. 'SI' ('NEG' ('TYPE' POROS) 'MCHAML ') ;
  488. 'ERREUR' ('CHAINE' 'L indice ' NOMDDT ' doit etre de type '
  489. 'FLOTTANT ou CHAMELEM') ;
  490. 'FINSI' ;
  491. 'SI' HYDRO ;
  492. 'SI' (('MINI' POROS) < 0.D0 ) ;
  493. 'ERREUR' 'Le coef. d emmagasinement doit etre positif' ;
  494. 'FINSI' ;
  495. 'FINSI' ;
  496. 'SI' TRANSP ;
  497. 'SI' ((('MINI' POROS) < 0.D0) 'OU' (('MAXI' POROS) > 1.D0)) ;
  498. 'ERREUR' 'La porosité doit etre comprise entre 0 et 1' ;
  499. 'FINSI' ;
  500. 'FINSI' ;
  501. POROS1 = 'KCHA' MODHYB POROS 'CHPO' ;
  502. * CONVECTION
  503. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' ) ;
  504. CONV1 = TRANSI . 'CONVECTION' ;
  505. 'FINSI' ;
  506. * RETARD
  507. 'SI' ( 'EXISTE' TRANSI 'COEF_RETARD' ) ;
  508. RETAR1 = TRANSI . 'COEF_RETARD' ;
  509. 'SI' ('EGA' ('TYPE' RETAR1) 'CHPOINT ') ;
  510. MINRET = ('MINI' RETAR1) ;
  511. 'SINON' ;
  512. MINRET = RETAR1 ;
  513. 'FINSI' ;
  514. 'SI' (MINRET < 1.D0) ;
  515. 'ERREUR' 'Le coefficient de retard doit etre supérieur à 1';
  516. 'FINSI' ;
  517. 'SINON' ;
  518. RETAR1 = 1. ;
  519. 'FINSI' ;
  520. RETAR1M1 = RETAR1 - 1. ;
  521. SAUVRET = 'EXISTE' TRANSI 'RETARD' ;
  522. * RETARD NON LINEAIRE
  523. RNONLINL = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
  524. ('EXISTE' TRANSI 'LANGMUIR') ;
  525. RNONLINF = ('EXISTE' TRANSI 'COEF_RETARD') 'ET'
  526. ('EXISTE' TRANSI 'FREUNDLICH') 'ET' ('NON' RNONLINL) ;
  527. RNONLIN = RNONLINL 'OU' RNONLINF ;
  528. * Cas Langmuir :
  529. 'SI' RNONLINL ;
  530. FSAT = TRANSI . 'LANGMUIR' ;
  531. 'SI' ('EGA' ('TYPE' FSAT) 'CHPOINT ') ;
  532. 'SI' (('MINI' FSAT) < 0.D0) ;
  533. 'ERREUR'
  534. 'La masse adsorbee a saturation doit etre positive' ;
  535. 'QUITTER' DARCYTRA ;
  536. 'FINSI' ;
  537. 'SINON' ;
  538. 'SI' (FSAT < 0.D0) ;
  539. 'ERREUR'
  540. 'La masse adsorbee a saturation doit etre positive' ;
  541. 'QUITTER' DARCYTRA ;
  542. 'FINSI' ;
  543. 'FINSI' ;
  544. RM1SURF = RETAR1M1 / FSAT ;
  545. 'FINSI' ;
  546. * Cas Freundlich :
  547. 'SI' RNONLINF ;
  548. UNSURN = TRANSI . 'FREUNDLICH' ;
  549. 'SI' ((UNSURN &lt;EG 0.D0) 'OU' (UNSURN > 1.D0)) ;
  550. 'ERREUR'
  551. 'L exposant de Freundlich doit etre compris entre 0'
  552. ' (exclu) et 1' ;
  553. 'QUITTER' DARCYTRA ;
  554. 'FINSI' ;
  555. 'FINSI' ;
  556. * DECROISSANCE
  557. DECRO = ( 'EXISTE' TRANSI 'DECROISSANCE' ) ;
  558. 'SI' DECRO ;
  559. LAMBD0 = TRANSI . 'DECROISSANCE' ;
  560. 'SI' (LAMBD0 < 0.D0) ;
  561. 'ERREUR' 'Le coefficient de décroissance doit etre positif.';
  562. 'FINSI' ;
  563. 'SINON' ;
  564. LAMBD0 = 0. ;
  565. 'FINSI' ;
  566. * DISSOLUTION_IMPOSEE
  567. SOLUA = ('EXISTE' TRANSI 'DISSOLUTION_IMPOSEE' ) ;
  568. 'SI' SOLUA ;
  569. DISIMP = TRANSI . 'DISSOLUTION_IMPOSEE' ;
  570. 'FINSI' ;
  571. * LIMITE_SOLUBILITE
  572. * (Priorité de la dissolution imposée sur les autres processus)
  573. SOLUL = (('EXISTE' TRANSI 'LIMITE_SOLUBILITE' ) 'ET' ('NON' SOLUA)) ;
  574. 'SI' SOLUL ;
  575. LIMSOL = TRANSI . 'LIMITE_SOLUBILITE' ;
  576. 'SI' (('MINI' LIMSOL) < 0.D0) ;
  577. 'ERREUR' 'La limite de solubilité doit etre positive.' ;
  578. 'FINSI' ;
  579. 'FINSI' ;
  580. * COEF_DISSOLUTION
  581. SOLUP = SOLUL 'ET' ('EXISTE' TRANSI 'COEF_DISSOLUTION' ) ;
  582. SOLUI = SOLUL 'ET' ('NON' SOLUP) ;
  583. 'SI' SOLUP ;
  584. CODIS = TRANSI . 'COEF_DISSOLUTION' ;
  585. 'SI' ('EGA' ('TYPE' CODIS) 'CHPOINT ') ;
  586. 'SI' (('MINI' CODIS) < 0.D0) ;
  587. 'ERREUR' 'Le coefficient de dissolution doit etre positif.';
  588. 'FINSI' ;
  589. 'SINON' ;
  590. 'SI' (CODIS < 0.D0) ;
  591. 'ERREUR' 'Le coefficient de dissolution doit etre positif.';
  592. 'FINSI' ;
  593. 'FINSI' ;
  594. 'FINSI' ;
  595. SOLUB = SOLUL 'OU' SOLUA ;
  596. *
  597. * ------------------------------------------------------------------ *
  598. * VERIFICATION DE LA BONNE STRUCTURE DES TABLES RESULTAT D'ORIGINE *
  599. * ------------------------------------------------------------------ *
  600. *
  601. * TEMPS
  602. 'SI' ( 'NON' ('EXISTE' TRANSI 'TEMPS' ) ) ;
  603. 'ERREUR'
  604. 'Indice TEMPS absent de la table de données.' ;
  605. 'QUITTER' DARCYTRA ;
  606. 'FINSI' ;
  607. * TRACE_INCONNUE
  608. 'SI' ( 'NON' ('EXISTE' TRANSI NOMTINC ) ) ;
  609. 'ERREUR' ('CHAINE'
  610. 'Indice ' NOMTINC ' absent de la table de données.') ;
  611. 'QUITTER' DARCYTRA ;
  612. 'FINSI' ;
  613. * INCONNUE
  614. 'SI' ( 'NON' ('EXISTE' TRANSI NOMINC ) ) ;
  615. 'ERREUR' ('CHAINE'
  616. 'Indice ' NOMINC ' absent de la table de données.') ;
  617. 'QUITTER' DARCYTRA ;
  618. 'FINSI' ;
  619. * FLUX
  620. 'SI' ( 'NON' ('EXISTE' TRANSI 'FLUX' ) ) ;
  621. 'ERREUR'
  622. 'Indice FLUX absent de la table de données.' ;
  623. 'QUITTER' DARCYTRA ;
  624. 'FINSI' ;
  625. * DISSOLUTION, PRECIPITE
  626. 'SI' SOLUB ;
  627. 'SI' ( 'NON' ('EXISTE' TRANSI 'DISSOLUTION' ) ) ;
  628. 'ERREUR'
  629. 'Indice DISSOLUTION absent de la table de données.' ;
  630. 'QUITTER' DARCYTRA ;
  631. 'FINSI' ;
  632. 'SI' ( 'NON' ('EXISTE' TRANSI 'PRECIPITE' ) ) ;
  633. 'ERREUR'
  634. 'Indice PRECIPITE absent de la table de données.' ;
  635. 'QUITTER' DARCYTRA ;
  636. 'FINSI' ;
  637. 'FINSI' ;
  638. * TEST DES TAILLES DE TABLE
  639. IND1 = 'INDEX' ( TRANSI . 'TEMPS' ) ;
  640. IND2 = 'INDEX' ( TRANSI . NOMTINC ) ;
  641. IND3 = 'INDEX' ( TRANSI . NOMINC ) ;
  642. IND4 = 'INDEX' ( TRANSI . 'FLUX' ) ;
  643. LIN1 = 'DIME' IND1 ;
  644. LIN2 = 'DIME' IND2 ;
  645. LIN3 = 'DIME' IND3 ;
  646. LIN4 = 'DIME' IND4 ;
  647. 'SI' ( LIN1 'NEG' LIN2 ) ;
  648. 'ERREUR' ('CHAINE'
  649. 'Longueur des tables TEMPS et ' NOMTINC ' différente.') ;
  650. 'QUITTER' DARCYTRA ;
  651. 'FINSI' ;
  652. 'SI' ( LIN1 'NEG' LIN3 ) ;
  653. 'ERREUR' ('CHAINE'
  654. 'Longueur des tables TEMPS et ' NOMINC ' différente.') ;
  655. 'QUITTER' DARCYTRA ;
  656. 'FINSI' ;
  657. 'SI' ( LIN1 'NEG' LIN4 ) ;
  658. 'ERREUR'
  659. 'Longueur des tables TEMPS et FLUX différente.' ;
  660. 'QUITTER' DARCYTRA ;
  661. 'FINSI' ;
  662. 'SI' SOLUB ;
  663. IND5 = 'INDEX' ( TRANSI . 'DISSOLUTION' ) ;
  664. IND6 = 'INDEX' ( TRANSI . 'PRECIPITE' ) ;
  665. LIN5 = 'DIME' IND5 ;
  666. LIN6 = 'DIME' IND6 ;
  667. 'SI' ( LIN1 'NEG' LIN5 ) ;
  668. 'ERREUR'
  669. 'Longueur des tables TEMPS et DISSOLUTION différente.' ;
  670. 'QUITTER' DARCYTRA ;
  671. 'FINSI' ;
  672. 'SI' ( LIN1 'NEG' LIN6 ) ;
  673. 'ERREUR'
  674. 'Longueur des tables TEMPS et PRECIPITE différente.' ;
  675. 'QUITTER' DARCYTRA ;
  676. 'FINSI' ;
  677. 'FINSI' ;
  678. * TEST DES INDICES DE TABLE
  679. IPO1 = 0 ;
  680. 'REPETER' BOU1 LIN1 ;
  681. IPO1 = IPO1 + 1 ;
  682. LAST1 = IND1 . IPO1 ;
  683. LAST2 = IND2 . IPO1 ;
  684. LAST3 = IND3 . IPO1 ;
  685. LAST4 = IND4 . IPO1 ;
  686. 'SI' SOLUB ;
  687. LAST5 = IND5 . IPO1 ;
  688. LAST6 = IND6 . IPO1 ;
  689. 'FINSI' ;
  690. 'SI' ( 'NEG' LAST1 LAST2 ) ;
  691. 'ERREUR'
  692. 'Indices des tables TEMPS et TRACE_CONC incohérents.' ;
  693. 'QUITTER' DARCYTRA ;
  694. 'FINSI' ;
  695. 'SI' ( 'NEG' LAST1 LAST3 ) ;
  696. 'ERREUR'
  697. 'Indices des tables TEMPS et CONCENTRATION incohérents.' ;
  698. 'QUITTER' DARCYTRA ;
  699. 'FINSI' ;
  700. 'SI' ( 'NEG' LAST1 LAST4 ) ;
  701. 'ERREUR'
  702. 'Indices des tables TEMPS et FLUX incohérents.' ;
  703. 'QUITTER' DARCYTRA ;
  704. 'FINSI' ;
  705. 'SI' SOLUB ;
  706. 'SI' ( 'NEG' LAST1 LAST5 ) ;
  707. 'ERREUR'
  708. 'Indices des tables TEMPS et DISSOLUTION incohérents.';
  709. 'QUITTER' DARCYTRA ;
  710. 'FINSI' ;
  711. 'SI' ( 'NEG' LAST1 LAST6 ) ;
  712. 'ERREUR'
  713. 'Indices des tables TEMPS et PRECIPITE incohérents.' ;
  714. 'QUITTER' DARCYTRA ;
  715. 'FINSI' ;
  716. 'FINSI' ;
  717. 'FIN' BOU1 ;
  718. *
  719. * ------------------------------------------------------------------ *
  720. * RECUPERATION DES CONDITIONS INITIALES OU DU DERNIER PAS SAUVE *
  721. * ------------------------------------------------------------------ *
  722. *
  723. TPSINI = TRANSI . 'TEMPS' . LAST1 ;
  724. TRP = TRANSI . NOMTINC . LAST1 ;
  725. CHRG = TRANSI . NOMINC . LAST1 ;
  726. FLU0 = TRANSI . 'FLUX' . LAST1 ;
  727. * On enlève les multiplicateurs de Lagrange éventuels
  728. LIST1 = 'EXTR' TRP 'COMP' ;
  729. LLIS1 = 'DIME' LIST1 ;
  730. 'SI' ( LLIS1 'NEG' 1 ) ;
  731. TRP = 'EXCO' TRP 'TH' 'TH' ;
  732. 'FINSI' ;
  733. *
  734. *--------------------------------------------------------------------*
  735. * RECUPERATION DES CONDITIONS AUX LIMITES ET DES CHARGEMENTS *
  736. *--------------------------------------------------------------------*
  737. *
  738. * CL de type Dirichlet : BLOCAGE
  739. * TRACE_IMPOSE
  740. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  741. MATBLOC = TRANSI . 'BLOCAGE' ;
  742. 'SI' ( 'EXISTE' TRANSI 'TRACE_IMPOSE') ;
  743. CHAIMP = TRANSI . 'TRACE_IMPOSE' ;
  744. 'SINON' ;
  745. 'ERREUR'
  746. 'Valeurs de traces de concentration imposées absentes' ;
  747. 'QUITTER' DARCYTRA ;
  748. 'FINSI' ;
  749. 'FINSI' ;
  750. * CL de type Neumann : FLUX_IMPOSE
  751. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' ) ;
  752. FLUIMP = TRANSI . 'FLUX_IMPOSE' ;
  753. 'FINSI' ;
  754. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  755. TERSOU = TRANSI . 'SOURCE' ;
  756. 'FINSI' ;
  757. *
  758. *--------------------------------------------------------------------*
  759. * RECUPERATION DES DONNEES NUMERIQUES *
  760. *--------------------------------------------------------------------*
  761. * THETA DIFFUSION
  762. 'SI' ( 'EXISTE' TRANSI NOMTETA ) ;
  763. TETA = TRANSI . NOMTETA ;
  764. 'SINON' ;
  765. TETA = 1.D0 ;
  766. 'FINSI' ;
  767. * THETA CONVECTION
  768. 'SI' ( 'EXISTE' TRANSI 'THETA_CONVECTION' ) ;
  769. TETAC = TRANSI . 'THETA_CONVECTION' ;
  770. 'SINON' ;
  771. TETAC = TETA ;
  772. 'FINSI' ;
  773. * PARAMETRES DE SOLUBILITE
  774. 'SI' (SOLUI 'OU' SOLUP) ;
  775. 'SI' ( 'EXISTE' TRANSI 'ITMAX_LIM') ;
  776. ITMAXI = TRANSI . 'ITMAX_LIM' ;
  777. 'SINON' ;
  778. ITMAXI = 50 ;
  779. 'FINSI' ;
  780. 'FINSI' ;
  781. 'SI' SOLUI ;
  782. 'SI' ( 'EXISTE' TRANSI 'PENALISATION') ;
  783. PENAL = TRANSI . 'PENALISATION' ;
  784. MET0 = 'PENALISATION' ;
  785. 'SINON' ;
  786. 'SI' ( 'EXISTE' TRANSI 'EPSI_LIM') ;
  787. EPS1 = TRANSI . 'EPSI_LIM' ;
  788. MET0 = 'PF DISSOLUTION' ;
  789. 'SINON' ;
  790. 'MESS' 'Il manque le paramètre ' ;
  791. 'ERREUR' 'du schéma numérique de limite de solubilité';
  792. 'FINSI' ;
  793. 'FINSI' ;
  794. 'FINSI' ;
  795. * PARAMETRES DE RETARD NON LINEAIRE
  796. 'SI' RNONLIN ;
  797. 'SI' ( 'EXISTE' TRANSI 'EPSI_RET' ) ;
  798. EPSRNL = TRANSI . 'EPSI_RET' ;
  799. 'SINON' ;
  800. EPSRNL = 1.D-4 ;
  801. 'FINSI' ;
  802. 'SI' ( 'EXISTE' TRANSI 'ITMAX_RET' ) ;
  803. ITMAXRNL = TRANSI . 'ITMAX_RET' ;
  804. 'SINON' ;
  805. ITMAXRNL = 20 ;
  806. 'FINSI' ;
  807. * Petit saut de concentration minimal pour le calcul de dérivée
  808. 'SI' ( 'EXISTE' TRANSI 'EPSI_COR' ) ;
  809. EPSCORD = TRANSI 'EPSI_COR' ;
  810. 'SINON' ;
  811. 'SI' RNONLINL ;
  812. EPSCORD = ('MANU' 'CHPO' HYTOT.CENTRE 1 'SCAL' 1.D-4)
  813. / RETAR1 * FSAT ;
  814. 'FINSI' ;
  815. 'SI' RNONLINF ;
  816. EPSCORD = 1.D-4 ;
  817. 'FINSI' ;
  818. 'FINSI' ;
  819. * Pour les isothermes linéaires la boucle de prise en compte de
  820. * retard non linéaire n'est parcourue qu'une fois :
  821. NPicard = ITMAXRNL ;
  822. 'SINON' ;
  823. NPicard = 1 ;
  824. 'FINSI' ;
  825. * THETA_DECROISSANCE
  826. 'SI' ( 'EXISTE' TRANSI 'THETA_DEC' ) ;
  827. BETA = TRANSI . 'THETA_DEC' ;
  828. 'SINON' ;
  829. BETA = 0.5 ;
  830. 'FINSI' ;
  831. * THETA_DISSOLUTION
  832. 'SI' ( 'EXISTE' TRANSI 'THETA_DIS' ) ;
  833. GAMMA = TRANSI . 'THETA_DIS' ;
  834. 'SINON' ;
  835. GAMMA = 1. ;
  836. 'FINSI' ;
  837. * TEMPS_CALCULES
  838. 'SI' ( 'EXISTE' TRANSI 'TEMPS_CALCULES' ) ;
  839. TPCAL = ORDO (TRANSI.'TEMPS_CALCULES') ;
  840. TRANSI . 'TEMPS_CALCULES' = TPCAL ;
  841. DCAL = 'DIME' TPCAL ;
  842. TPSFIN = 'EXTR' TPCAL DCAL ;
  843. TPS0 = 'EXTR' TPCAL 1 ;
  844. 'SI' ('NEG' DCAL 1) ;
  845. TPSANT = 'EXTR' TPCAL (DCAL - 1) ;
  846. TPS1 = 'EXTR' TPCAL 2 ;
  847. DT10 = TPS1 - TPS0 ;
  848. DT21 = TPSFIN - TPSANT ;
  849. PROPS = 'PROG' DT10 DT21 ;
  850. DTEPS = 'MINIMUM' PROPS ;
  851. 'SINON' ;
  852. DTEPS = TPSFIN - TPS0 ;
  853. 'FINSI' ;
  854. EPS0 = DTEPS * 1.D-6 ;
  855. *
  856. * ICAL est le premier indice des temps à calculer supérieur à
  857. * TpsIni + eps
  858. IOK1 = 0 ;
  859. 'REPETER' BOU2 DCAL ;
  860. ICAL = &BOU2 ;
  861. TEMS = 'EXTR' TPCAL ICAL ;
  862. 'SI' ( TPSINI '<' (TEMS - EPS0) ) ;
  863. IOK1 = 1 ;
  864. 'QUITTER' BOU2 ;
  865. 'FINSI' ;
  866. 'FIN' BOU2 ;
  867. 'SI' ( IOK1 'EGA' 0) ;
  868. 'MESS' 'Listreel des temps de calcul inférieur à tini.' ;
  869. 'MESS' 'On ne fait donc rien. ' ;
  870. 'QUITTER' DARCYTRA ;
  871. 'FINSI' ;
  872. DELTAT0 = TEMS - TPSINI ;
  873. 'SINON' ;
  874. 'ERREUR'
  875. 'Indice TEMPS_CALCULES absent de la table de données.' ;
  876. 'QUITTER' DARCYTRA ;
  877. 'FINSI' ;
  878. * TEMPS_SAUVES
  879. 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ;
  880. TPSOR = ORDO (TRANSI . 'TEMPS_SAUVES') ;
  881. TRANSI . 'TEMPS_SAUVES' = TPSOR ;
  882. DSOR = 'DIME' TPSOR ;
  883. ISOR = 0 ;
  884. IOK1 = 0 ;
  885. 'REPETER' BOU3 DSOR ;
  886. ISOR = ISOR + 1 ;
  887. TEMS = 'EXTR' TPSOR ISOR ;
  888. 'SI' ( ( TPSINI '<' (TEMS - EPS0) ) 'ET'
  889. ( TEMS '&lt;EG' (TPSFIN +EPS0)) ) ;
  890. IOK1 = 1 ;
  891. 'QUITTER' BOU3 ;
  892. 'FINSI' ;
  893. 'FIN' BOU3 ;
  894. 'SI' ( IOK1 'EGA' 0) ;
  895. 'ERREUR' 'Listreel des temps de sauvegarde hors tmin,tmax.';
  896. 'QUITTER' DARCYTRA ;
  897. 'FINSI' ;
  898. 'FINSI' ;
  899. *
  900. *--------------------------------------------------------------------*
  901. * RAPPEL DES DONNEES ENTREES DANS LA PROCEDURE *
  902. *--------------------------------------------------------------------*
  903. *
  904. 'MESS' ' ' ;
  905. 'MESS' 'MODELISATION Darcy EFMH en transitoire.' ;
  906. 'MESS' '---------------------------------------' ;
  907. 'MESS' ' ' ;
  908. 'MESS' 'Donnees présentes en entrée : ' ;
  909. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  910. 'MESS' 'Calcul avec conditions aux limites sur Th (Dirichlet)' ;
  911. 'FINSI' ;
  912. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' ) ;
  913. 'MESS' 'Calcul avec conditions aux limites de flux (Neumann)' ;
  914. 'FINSI' ;
  915. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' ) ;
  916. 'MESS' 'Ce problème comporte un terme convectif' ;
  917. 'FINSI' ;
  918. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  919. 'MESS' 'Ce problème comporte un terme source' ;
  920. 'FINSI' ;
  921. 'SI' DECRO ;
  922. 'MESS' 'Ce problème comporte un terme de décroissance' ;
  923. 'FINSI' ;
  924. 'SI' RNONLIN ;
  925. 'MESS' 'Ce problème comporte la prise en compte d isothermes '
  926. 'd adsorption non-linéaires' ;
  927. 'SI' RNONLINL ;
  928. 'MESS' ' de type Langmuir (méthode de Picard)' ;
  929. 'FINSI' ;
  930. 'SI' RNONLINF ;
  931. 'MESS' ' de type Freundlich (méthode de Picard)' ;
  932. 'FINSI' ;
  933. 'FINSI' ;
  934. 'SI' SOLUL ;
  935. 'MESS' 'Ce problème comporte une limite de solubilité' ;
  936. 'SI' SOLUP ;
  937. 'MESS' ' avec dissolution progressive d ordre 1' ;
  938. 'FINSI' ;
  939. 'SI' (EGA MET0 'PENALISATION') ;
  940. 'MESS' ' avec coefficient de pénalisation = ' PENAL ;
  941. 'FINSI' ;
  942. 'SI' (EGA MET0 'PF DISSOLUTION') ;
  943. 'MESS' ' avec critère de convergence = ' EPS1 ;
  944. 'FINSI' ;
  945. 'FINSI' ;
  946. 'SI' SOLUA ;
  947. 'MESS' 'Ce probleme comporte une dissolution imposée' ;
  948. 'FINSI' ;
  949. 'MESS' ' ' ;
  950. 'MESS' 'Valeur des paramètres des schémas numériques :' ;
  951. 'MESS' ' (0:Schéma explicite, 0.5:Crank-Nicholson, 1:Implicite)' ;
  952. 'MESS' ' diffusion : ' TETA ;
  953. 'SI' ( 'EXISTE' TRANSI 'CONVECTION' ) ;
  954. 'MESS' ' convection : ' TETAC ;
  955. 'FINSI' ;
  956. 'SI' DECRO ;
  957. 'MESS' ' décroissance : ' BETA ;
  958. 'FINSI' ;
  959. 'SI' SOLUP ;
  960. 'MESS' ' dissolution d ordre 1 : ' GAMMA ;
  961. 'FINSI' ;
  962. 'MESS' ' ' ;
  963. 'MESS' 'Valeur du temps initial : ' TPSINI ;
  964. 'MESS' 'Valeur du temps final : ' TPSFIN ;
  965. *
  966. *
  967. *--------------------------------------------------------------------*
  968. * BOUCLE RESOLVANT LE SYSTEME POUR CHAQUE PAS DE TEMPS *
  969. *--------------------------------------------------------------------*
  970. *
  971. *==============
  972. * Préparation préliminaire de la boucle sur les pas de temps :
  973. *==============
  974. *
  975. NBIT = 0 ;
  976. NBITRNL = 0 ;
  977. MASHYB = 'MHYB' MODHYB MAT1 ;
  978. PRECED = TPSINI ;
  979. IPAS = ICAL ;
  980. DELOLD = 0.D0 ;
  981. TPS1 = 'EXTR' TPCAL ICAL ;
  982. DELTAT = TPS1 - TPSINI ;
  983. 'MESS' 'Incrément de temps initial : ' DELTAT ;
  984. 'MESS' ' ' ;
  985. *
  986. * Initialisation de la quantité dissoute
  987. *
  988. 'SI' SOLUB ;
  989. * DIS et DIS1 sont des valeurs de dissolution valables pour
  990. * l'ensemble d'une maille, et intégré sur tout le pas de temps
  991. DIS = (TRANSI . 'DISSOLUTION' . LAST1) * VOLU1 * DELTAT0 ;
  992. PRCI = TRANSI . 'PRECIPITE' . LAST1 ;
  993. ANT0 = 'EXCO' 'H' ('MASQ' PRCI 'SUPERIEUR' 1.D-14) 'SCAL' ;
  994. 'FINSI' ;
  995. *
  996. * Calcul de la fraction adsorbée initiale :
  997. *
  998. CC0 = 'NOMC' CHRG 'SCAL' ;
  999. 'SI' RNONLIN ;
  1000. 'SI' RNONLINL ;
  1001. * Langmuir :
  1002. FF0 = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 )) ;
  1003. 'FINSI' ;
  1004. 'SI' RNONLINF ;
  1005. * Freundlich :
  1006. * La bidouille utilisant les masques sert à obtenir le champ point
  1007. * des signes de CC0.
  1008. a = 'MASQUE' CC0 'SUPERIEUR' 0. ;
  1009. b = 'MASQUE' CC0 'INFERIEUR' 0. ;
  1010. FF0 = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b) ;
  1011. 'DETRUIT' a ;
  1012. 'DETRUIT' b ;
  1013. 'FINSI' ;
  1014. 'SINON' ;
  1015. * Isotherme linéaire :
  1016. FF0 = RETAR1M1 * CC0 ;
  1017. 'FINSI' ;
  1018. FF = FF0 ;
  1019.  
  1020. *
  1021. * Calcul du coefficient de retard initial :
  1022. * Le facteur de retard s'exprime comme
  1023. * R(C) = 1 + ( F'(Ct+dt + eps) - F'(Ct) ) / ( Ct+dt + eps - Ct) )
  1024. * C'est la méthode de la corde employée ici pour ses vertus
  1025. * stabilisantes. Le petit 'eps' (EpsCord) intervient pour évacuer les
  1026. * questions de conditions initiales ou de régime permanent.
  1027. 'SI' RNONLIN ;
  1028. CC1 = CC0 + EpsCord ;
  1029. 'SI' RNONLINL ;
  1030. * Langmuir
  1031. FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 )) ;
  1032. 'SINON' ;
  1033. * Freundlich
  1034. a = 'MASQUE' CC1 'SUPERIEUR' 0. ;
  1035. b = 'MASQUE' CC1 'INFERIEUR' 0. ;
  1036. FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b) ;
  1037. 'DETRUIT' a ;
  1038. 'DETRUIT' b ;
  1039. 'FINSI' ;
  1040. RETARC = 1.D0 + ( (FF1 - FF0) / EpsCord ) ;
  1041. 'DETRUIT' CC1 ;
  1042. 'DETRUIT' FF1 ;
  1043. 'SINON' ;
  1044. RETARC = RETAR1 ;
  1045. 'FINSI' ;
  1046. RETAR0 = RETARC ;
  1047. 'SI' SAUVRET ;
  1048. TRANSI . 'RETARD' . 0 = RETARC ;
  1049. 'FINSI' ;
  1050. *
  1051. * Fonctions de volume :
  1052. PORSUF1 = POROS1 * VOLU1 ;
  1053. PORETSU1 = PORSUF1 * RETARC ;
  1054. DECROI0 = 'EXCO' ( (-1.) * LAMBD0 * PORSUF1 ) 'CK' 'SCAL' ;
  1055. *
  1056. * Indicateur stipulant si l'on doit recalculer la matrice hybride :
  1057. * Au départ, vrai car il faut la calculer au moins la première fois,
  1058. * puis faux par la suite jusqu-à notification contraire.
  1059. Recalhyb = VRAI ;
  1060. *
  1061. * Table utilisée par les opérateurs MATP, SMTP, HYBP et HDEB
  1062. * TAB.'SURF' est le coefficient devant dC/dt (défini dans chaque module).
  1063. * On y intègre la contribution de la décroissance implicite
  1064. * ainsi que celle du terme de pénalisation
  1065. TAB = 'TABLE' ;
  1066. TAB.'SOUSTYPE' = 'DARCY_TRANSITOIRE' ;
  1067. TAB.'THETA' = TETA ;
  1068. TAB.'THETA_CONVECTION' = TETAC ;
  1069. *
  1070. *
  1071. *=================================
  1072. 'REPETE' BOUTPS (DCAL - ICAL + 1); COMM 'Boucle sur le temps' ;
  1073. *=================================
  1074. *
  1075. IPAS = ICAL + &BOUTPS - 1 ;
  1076. *
  1077. *-- Initialisation en-tête de boucle sur le temps
  1078. *
  1079. TPS = 'EXTR' TPCAL IPAS ;
  1080. MESS 'TEMPS T=' TPS ;
  1081. DELTAT = TPS - PRECED ;
  1082. EPSDT = DELTAT + DELOLD / 2.D0 * 1.D-6 ;
  1083. TAB.'PAS' = DELTAT ;
  1084. TAB.'TRACE' = TRP ;
  1085. TAB.'CHARGE' = CHRG ;
  1086. Recalhyb = Recalhyb 'OU' ( DELOLD 'NEG' DELTAT EPSDT ) ;
  1087.  
  1088. *
  1089. *- Calcul de la contribution des termes sources et de la convection
  1090. *
  1091. TERSCE = 'MANU' 'CHPO' HYTOT.CENTRE 1 'SOUR' 0. ;
  1092. *
  1093. * terme source de convection :
  1094. 'SI' ((TETAC < 1.D-4) 'ET' ( 'EXISTE' TRANSI 'CONVECTION' )) ;
  1095. CONC1 = 'NOMC' 'SCAL' TRP ;
  1096. CONC2 = (-1.) * CONC1 * CONV1 ;
  1097. TERSC1 = 'DIVU' MODHYB CONC2 MCHYB ;
  1098. TERSC1 = 'NOMC' 'SOUR' TERSC1 ;
  1099. TERSCE = TERSC1 + TERSCE ;
  1100. 'FINSI' ;
  1101. *
  1102. * Terme source propre :
  1103. 'SI' ( 'EXISTE' TRANSI 'SOURCE' ) ;
  1104. TERSC2 = 'EXCO' ('TIRE' TERSOU TPS) 'SOUR' 'SOUR' ;
  1105. TERSCE = TERSC2 + TERSCE ;
  1106. 'FINSI' ;
  1107. *
  1108. * Terme source de décroissance explicite :
  1109. * du soluté et de l'adsorbat
  1110. 'SI' DECRO ;
  1111. TERSC3 = ( ('NOMC' CHRG 'SCAL') + FF ) * DECROI0 ;
  1112. TERSC3 = 'NOMC' 'SOUR' TERSC3 ;
  1113. TERSCE = TERSC3 + TERSCE ;
  1114. 'FINSI' ;
  1115. *
  1116. *- Incorporation des CLs
  1117. *
  1118. SMTR0 = 'MANU' 'CHPO' HYTOT.FACE 1 'FLUX' 0. 'NATURE' 'DISCRET' ;
  1119. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1120. CHARIMPO = 'TIRE' CHAIMP TPS ;
  1121. CHARIMPO = 'CHAN' 'ATTRIBUT' CHARIMPO 'NATURE' 'DISCRET' ;
  1122. SMTR0 = SMTR0 'ET' CHARIMPO ;
  1123. 'FINSI' ;
  1124. *
  1125. 'SI' ( 'EXISTE' TRANSI 'FLUX_IMPOSE' ) ;
  1126. FLUIMPO = 'TIRE' FLUIMP TPS ;
  1127. FLUIMPO = 'CHAN' 'ATTRIBUT' FLUIMPO 'NATURE' 'DISCRET' ;
  1128. SMTR0 = SMTR0 'ET' FLUIMPO ;
  1129. 'FINSI' ;
  1130.  
  1131. *
  1132. *
  1133. *|-------------------------------------------------------------------|
  1134. *| Boucle de Retard Non Linéaire |
  1135. *|-------------------------------------------------------------------|
  1136. *
  1137. * Boucle de Picard pour la prise en compte
  1138. * d'isothermes d'adsorption non linéaires : isothermes de Langmuir et
  1139. * de Freundlich. Ce point nécessite d'effectuer une boucle sur
  1140. * le facteur de retard qui devient fonction de la concentration
  1141. * en solution.
  1142. *
  1143. *------------------------
  1144. 'REPETE' BOURNL NPicard ; COMM 'Boucle sur le coef. de retard' ;
  1145. *------------------------
  1146. *
  1147. IRNL = &BOURNL ;
  1148. *
  1149. * Quatre grands blocs selon que l'on traite :
  1150. * -------------------------------------------
  1151. * -1- sans limite de solubilité ou avec dissolution arbitraire
  1152. * -2- avec limite de solubilité et dissolution d'ordre 1
  1153. * -3- avec limite suivant méthode de pénalisation,
  1154. * -4- avec limite suivant schéma prédicteur/correcteur itératif
  1155.  
  1156.  
  1157. *|----------------------------------------------------|
  1158. *| Pas de limite de solubilité ou Dissolution imposée |
  1159. *|----------------------------------------------------|
  1160. 'SI' (('NON' SOLUB) 'OU' SOLUA) ;
  1161. *
  1162. * On recalcule la matrice hybride si DeltaT a changé :
  1163. 'SI' Recalhyb ;
  1164. * Coefficients correspondant à la décroissance expl. et implicite :
  1165. NUM0 = (BETA - 1.) * LAMBD0 * DELTAT + 1. ;
  1166. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1167. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1168. COFDC = PORETSU1 * DEN0 ;
  1169. TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM' ;
  1170. * Matrice du problème
  1171. HNDTR = 'MATP' MODHYB MASHYB TAB ;
  1172. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1173. HNDT1 = HNDTR 'ET' MATBLOC ;
  1174. 'SINON' ;
  1175. HNDT1 = HNDTR ;
  1176. 'FINSI' ;
  1177. * On ajoute la rigidité associée à la convection implicite
  1178. * Elle ne dépend ni du terme source, ni de la masse hybride,
  1179. * ni de la valeur de la trace de concentration. On doit néanmoins
  1180. * fournir des valeurs de ces grandeurs :
  1181. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1182. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCE CONV1 ;
  1183. HNDT1 = HNDT1 'ET' RIGC ;
  1184. 'FINSI' ;
  1185. Recalhyb = FAUX ;
  1186. 'FINSI' ;
  1187. *
  1188. 'SI' SOLUA ;
  1189. * Evaluation des termes sources de dissolution :
  1190. * Dissolution arbitraire (ordre 0 par exemple), limitée au total précipité :
  1191. DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1 ;
  1192. DISSOLP = 'TIRE' DISIMP TPS ;
  1193. TERSC4 = 0.5 * (DISSOLT + DISSOLP - ('ABS' (DISSOLT - DISSOLP)));
  1194. TERSC4 = 'EXCO' 'H' TERSC4 'SOUR' ;
  1195. TERSCT = TERSCE + TERSC4 ;
  1196. 'SINON' ;
  1197. TERSCT = TERSCE ;
  1198. 'FINSI' ;
  1199.  
  1200. * Rigidité et flux convectifs :
  1201. * SMTR0 contient la contribution en termes de flux des conditions aux
  1202. * limites
  1203. * Dans tous les cas, SMTR1 contient les termes sources et la convection,
  1204. * mais TERSCT ne contient la convection que si TETAC > 10-4
  1205. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1206. * TERSCT contient tous les T.S. sauf la convection,
  1207. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT CONV1 ;
  1208. 'SINON' ;
  1209. * TERSCT contient tous les T.S. convection comprise
  1210. SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT ;
  1211. 'FINSI' ;
  1212. SMTR = SMTR0 + SMTR1 ;
  1213.  
  1214. *- Résolution en trace de concentration :
  1215. TP1 = 'RESO' HNDT1 SMTR ;
  1216. TP2 = 'EXCO' TP1 'TH' 'TH' ;
  1217.  
  1218. *- Terme source de convection :
  1219. * TERSCVT contient tous les T.S. y compris la convection
  1220. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1221. * On l'ajoute au terme source total maintenant seulement car on
  1222. * a besoin de C(t+dt).
  1223. CONC1 = ((1.D0 - TETAC) * TRP) + (TETAC * TP2) ;
  1224. CONC2 = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1 ;
  1225. CONC3 = 'DIVU' MODHYB CONC2 MCHYB ;
  1226. TERSCV = 'NOMC' 'SOUR' CONC3 ;
  1227. TERSCVT = TERSCT + TERSCV ;
  1228. 'SINON' ;
  1229. * Le terme source de convection est déjà dans TERSCT
  1230. TERSCVT = TERSCT ;
  1231. 'FINSI' ;
  1232.  
  1233. * Calcul des nouveaux concentration, flux et précipité :
  1234. CHA2 = 'HYBP' MODHYB MASHYB TP2 TAB TERSCVT ;
  1235. FLU2 = 'HDEB' MODHYB MASHYB CHA2 TP2 ;
  1236. 'SI' SOLUA ;
  1237. DIS1 = 'EXCO' 'SOUR' (TERSC4 * DELTAT) 'H' ;
  1238. PRE2 = ((NUM0 * PRCI) - (DIS1 / VOLU1)) / DEN0 ;
  1239. 'FINSI' ;
  1240.  
  1241. * 'Fin module dissolution arbitraire'
  1242. 'FINSI' ;
  1243.  
  1244. *|---------------------------|
  1245. *| Dissolution d'ordre 1 |
  1246. *|---------------------------|
  1247. 'SI' SOLUP;
  1248.  
  1249. * On recalcule la matrice hybride si DeltaT a changé :
  1250. 'SI' Recalhyb ;
  1251. * Coefficients correspondant à la décroissance expl. et implicite :
  1252. NUM0 = (BETA - 1.) * LAMBD0 * DELTAT + 1. ;
  1253. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1254. * Terme implicite de dissolution :
  1255. 'SI' (GAMMA > 1.D-4) ;
  1256. DEN1 = ANT0 * GAMMA * CODIS * DELTAT ;
  1257. 'SINON' ;
  1258. DEN1 = 0. ;
  1259. 'FINSI' ;
  1260. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1261. COFDC = PORETSU1 * (DEN0 + DEN1) ;
  1262. TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM' ;
  1263. * Matrice du problème
  1264. HNDTR = 'MATP' MODHYB MASHYB TAB ;
  1265. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1266. HNDT1 = HNDTR 'ET' MATBLOC ;
  1267. 'SINON' ;
  1268. HNDT1 = HNDTR ;
  1269. 'FINSI' ;
  1270. * On ajoute la rigidité associée à la convection implicite
  1271. * Elle ne dépend ni du terme source, ni de la masse hybride,
  1272. * ni de la valeur de la trace de concentration. On doit néanmoins
  1273. * fournir des valeurs de ces grandeurs :
  1274. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1275. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCE CONV1 ;
  1276. HNDT1 = HNDT1 'ET' RIGC ;
  1277. 'FINSI' ;
  1278. NBIT = NBIT + 1 ;
  1279. Recalhyb = FAUX ;
  1280. 'FINSI' ;
  1281. *
  1282. * Evaluation des termes sources de dissolution (ici partie explicite) :
  1283. * terme de dissolution totale puis terme de dissolution d'ordre 1
  1284. * ramenés à une unité de temps :
  1285. DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1 ;
  1286. DIFCHG = LIMSOL - CHRG ;
  1287. DISSOLP = ('EXCO' 'CK' PORETSU1 'SCAL') * CODIS * (LIMSOL - CHRG) ;
  1288.  
  1289. *---------------
  1290. 'REPETER' BOU8 ; COMM 'Boucle dissolution ordre 1' ;
  1291. *---------------
  1292. * On prend celui des deux termes correspondant à la présence de
  1293. * précipité de l'itération d'avant.
  1294. TERSC4 = 'EXCO' 'H' (((1.D0 - ANT0) * DISSOLT)
  1295. + (ANT0 * DISSOLP)) 'SOUR' ;
  1296. TERSCT = TERSCE + TERSC4 ;
  1297.  
  1298. * Rigidité et flux convectifs :
  1299. * SMTR0 contient la contribution en termes de flux des conditions aux
  1300. * limites
  1301. * Dans tous les cas, SMTR1 contient les termes sources et la
  1302. * convection, mais TERSCT ne contient la convection que si TETAC>10-4
  1303. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1304. * TERSCT contient tous les T.S. sauf la convection,
  1305. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT CONV1 ;
  1306. 'SINON' ;
  1307. * TERSCT contient tous les T.S. convection comprise
  1308. SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT ;
  1309. 'FINSI' ;
  1310. SMTR = SMTR0 + SMTR1 ;
  1311.  
  1312. *- Résolution en trace de concentration :
  1313. TP1 = 'RESO' HNDT1 SMTR ;
  1314. TP2 = 'EXCO' TP1 'TH' 'TH' ;
  1315.  
  1316. *- Terme source de convection :
  1317. * TERSCVT contient tous les T.S. y compris la convection
  1318. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4));
  1319. * On l'ajoute au terme source total maintenant seulement car on
  1320. * a besoin de C(t+dt).
  1321. CONC1 = ((1.D0 - TETAC) * TRP) + (TETAC * TP2) ;
  1322. CONC2 = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1 ;
  1323. CONC3 = 'DIVU' MODHYB CONC2 MCHYB ;
  1324. TERSCV = 'NOMC' 'SOUR' CONC3 ;
  1325. TERSCVT = TERSCT + TERSCV ;
  1326. 'SINON' ;
  1327. * Le terme source de convection est déjà dans TERSCT
  1328. TERSCVT = TERSCT ;
  1329. 'FINSI' ;
  1330.  
  1331. * Calcul des nouveaux concentration et précipité :
  1332. CHA2 = 'HYBP' MODHYB MASHYB TP2 TAB TERSCVT ;
  1333. CHA3 = (1.D0 - GAMMA) * CHRG + (GAMMA * CHA2) ;
  1334. DIS3 = ('EXCO' 'CK' PORETSU1 'SCAL') / VOLU1
  1335. * CODIS * (LIMSOL - CHA3) * DELTAT;
  1336. PRE2 = ANT0 / DEN0 * ( (NUM0 * PRCI) - DIS3 ) ;
  1337. * On calcule le flux ici pour etre sur qu'il correspond aux données
  1338. * en cours de la résolution. En effet, en cas de sortie sans
  1339. * convergence, si on avait calculé le flux après les itérations,
  1340. * MCHYB et TAB.SURF auraient tenu compte de la nouvelle distribution
  1341. * ANT2 alors que les calculs ont été effectués avec la courante ANT0.
  1342. FLU2 = 'HDEB' MODHYB MASHYB CHA2 TP2 ;
  1343.  
  1344. * Critère de sortie et mise à jour de la nouvelle distribution :
  1345. DIFOK = (NUM0 * PRCI) - DIS3 ;
  1346. ANT2 = 'EXCO' 'H' ('MASQ' DIFOK 'SUPERIEUR' 0.) 'SCAL' ;
  1347. CRIT = 'ABS' (ANT2 - ANT0) ;
  1348. ANT0 = ANT2 ;
  1349. 'SI' ( ('MAXI' CRIT) < 1.D-14) ;
  1350. 'QUITTER' BOU8 ;
  1351. 'SINON' ;
  1352. * on reconstruit la matrice hybride
  1353. 'SI' (GAMMA > 1.D-4) ;
  1354. DEN1 = ANT0 * GAMMA * CODIS * DELTAT ;
  1355. 'SINON' ;
  1356. DEN1 = 0. ;
  1357. 'FINSI' ;
  1358. COFDC = PORETSU1 * (DEN0 + DEN1) ;
  1359. TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM' ;
  1360. HNDTR = 'MATP' MODHYB MASHYB TAB ;
  1361. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1362. HNDT1 = HNDTR 'ET' MATBLOC ;
  1363. 'SINON' ;
  1364. HNDT1 = HNDTR ;
  1365. 'FINSI' ;
  1366. * N'oublions pas la rigidité associée à la convection implicite
  1367. * calculée avec le nouveau terme de pénalisation
  1368. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4));
  1369. HNDT1 = HNDT1 'ET' RIGC ;
  1370. 'FINSI' ;
  1371. NBIT = NBIT + 1 ;
  1372. 'FINSI' ;
  1373. 'SI' (&BOU8 EGA ITMAXI) ;
  1374. 'MESS' 'Sortie sans convergence' ;
  1375. 'QUITTER' BOU8 ;
  1376. 'FINSI' ;
  1377.  
  1378. *-----------
  1379. 'FIN' BOU8 ;
  1380. *-----------
  1381.  
  1382. * Calcul de la dissolution :
  1383. DIS1 = ((NUM0 * PRCI) - (PRE2 * DEN0)) * VOLU1 ;
  1384.  
  1385. NB0 = &BOU8 - 1 ;
  1386. 'SI' (NB0 > 1);
  1387. 'MESS' ' ' NB0 ' itérations de solubilité' ;
  1388. 'FINSI';
  1389.  
  1390. * 'Fin module Dissolution d ordre 1'
  1391. 'FINSI' ;
  1392.  
  1393. *|-------------------------------------------------------------|
  1394. *| Dissolution instantanée - Point fixe Présence précipité |
  1395. *|-------------------------------------------------------------|
  1396. 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION')) ;
  1397.  
  1398. * On recalcule la matrice hybride si DeltaT a changé :
  1399. 'SI' Recalhyb ;
  1400. * Coefficients correspondant à la décroissance expl. et implicite :
  1401. NUM0 = (BETA - 1.) * LAMBD0 * DELTAT + 1. ;
  1402. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1403. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1404. PEN0 = ANT0 * PENAL * DELTAT ;
  1405. COFDC = PORETSU1 + (('EXCO' 'SCAL' PEN0 'CK') * VOLU1) * DEN0 ;
  1406. TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM' ;
  1407. * Matrice du problème
  1408. HNDTR = 'MATP' MODHYB MASHYB TAB ;
  1409. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1410. HNDT1 = HNDTR 'ET' MATBLOC ;
  1411. 'SINON' ;
  1412. HNDT1 = HNDTR ;
  1413. 'FINSI' ;
  1414. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1415. * On ajoute la rigidité associée à la convection implicite
  1416. * Elle ne dépend ni du terme source, ni de la masse hybride,
  1417. * ni de la valeur de la trace de concentration. On doit néanmoins
  1418. * fournir des valeurs de ces grandeurs :
  1419. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCE CONV1 ;
  1420. HNDT1 = HNDT1 'ET' RIGC ;
  1421. 'FINSI' ;
  1422. NBIT = NBIT + 1 ;
  1423. Recalhyb = FAUX ;
  1424. 'FINSI' ;
  1425. *
  1426. *- Résolution en trace de concentration et calcul des inconnues secondaires
  1427. * Méthode du point fixe sur la présence de précipité
  1428. * par pénalisation (dissolution instantanée)
  1429.  
  1430. * Décroissance et dissolution du précipité avant les itérations
  1431. * La filiation du précipité est en principe comprise dans
  1432. * celle du soluté (à la charge de l'utilisateur)
  1433. * On prendra au début la répartition de précipité du pas précédent
  1434. DISSOLT = PRCI * (NUM0 / DELTAT) * VOLU1 ;
  1435. TERSC4 = 'EXCO' 'H' DISSOLT 'SOUR' ;
  1436. TERSCE = TERSCE + TERSC4 ;
  1437. TERSCT = TERSCE ;
  1438.  
  1439. *--------------
  1440. 'REPETE' BOU6 ; 'COMM Boucle limite de sol., Pénalisation' ;
  1441. *--------------
  1442. * Terme de pénalisation :
  1443. TERSC5 = PEN0 * (DEN0 / DELTAT) * (LIMSOL - CHRG) * VOLU1 ;
  1444. TERSC5 = 'EXCO' 'H' TERSC5 'SOUR' ;
  1445. TERSCT = TERSCE + TERSC5 ;
  1446.  
  1447. * Calcul du second membre :
  1448. * SMTR0 contient la contribution en termes de flux des conditions aux
  1449. * limites. Dans tous les cas, SMTR1 contient le second membre du
  1450. * à la diffusion, les termes sources et le second membre du à
  1451. * la convection, mais TERSCT ne contient la convection que si
  1452. * TETAC est supérieur à 10-4
  1453. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1454. * TERSCT contient tous les T.S. sauf la convection, on
  1455. * regénére la rigidité au passage mais on ne s'en sert pas.
  1456. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT CONV1 ;
  1457. 'SINON' ;
  1458. * TERSCT contient tous les T.S. convection comprise
  1459. SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT ;
  1460. 'FINSI' ;
  1461. SMTR = SMTR0 + SMTR1 ;
  1462.  
  1463. *- Résolution :
  1464. TP1 = 'RESO' HNDT1 SMTR ;
  1465. TP2 = 'EXCO' TP1 'TH' 'TH' ;
  1466.  
  1467. *- Terme source de convection pour calculer la nouvelle concentration
  1468. * TERSCVT contient tous les T.S. y compris la convection
  1469. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1470. * On l'ajoute au terme source total maintenant seulement car on
  1471. * a besoin de C(t+dt).
  1472. CONC1 = ((1.D0 - TETAC) * TRP) + (TETAC * TP2) ;
  1473. CONC2 = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1 ;
  1474. CONC3 = 'DIVU' MODHYB CONC2 MCHYB ;
  1475. TERSCV = 'NOMC' 'SOUR' CONC3 ;
  1476. TERSCVT = TERSCT + TERSCV ;
  1477. 'SINON' ;
  1478. * Le terme source de convection est déjà dans TERSCT
  1479. TERSCVT = TERSCT ;
  1480. 'FINSI' ;
  1481. CHA2 = 'HYBP' MODHYB MASHYB TP2 TAB TERSCVT ;
  1482. DIFCHG = CHA2 - LIMSOL ;
  1483. PRE2 = PEN0 * DIFCHG ;
  1484. * On calcule le flux ici pour etre sur qu'il correspond aux données
  1485. * en cours de la résolution. En effet, en cas de sortie sans
  1486. * convergence, si on avait calculé le flux après les itérations,
  1487. * MCHYB et TAB.SURF auraient tenu compte de la nouvelle distribution
  1488. * ANT2 alors que les calculs ont été effectués avec la courante ANT0.
  1489. FLU2 = 'HDEB' MODHYB MASHYB CHA2 TP2 ;
  1490.  
  1491. * Critère de sortie de la boucle d'itération :
  1492. MASK1 = 'MASQ' DIFCHG 'SUPERIEUR' 0. ;
  1493. MASK2 = 'MASQ' PRE2 'SUPERIEUR' 0. ;
  1494. ANT2 = 'EXCO' 'H'
  1495. ('MASQ' (MASK1 + MASK2) 'EGSUPE' .5) 'SCAL' ;
  1496. CRIT = 'ABS' (ANT2 - ANT0) ;
  1497. ANT0 = ANT2 ;
  1498. 'SI' ( ('MAXI' CRIT) < 1.D-14) ;
  1499. 'QUITTER' BOU6 ;
  1500. 'SINON' ;
  1501. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1502. PEN0 = ANT0 * PENAL * DELTAT ;
  1503. COFDC = PORETSU1 + (('EXCO' 'SCAL' PEN0 'CK') * VOLU1) * DEN0;
  1504. TAB.'SURF' = 'KCHA' MODHYB COFDC 'CHAM' ;
  1505. * Matrice du problème
  1506. HNDTR = 'MATP' MODHYB MASHYB TAB ;
  1507. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1508. HNDT1 = HNDTR 'ET' MATBLOC ;
  1509. 'SINON' ;
  1510. HNDT1 = HNDTR ;
  1511. 'FINSI' ;
  1512. * N'oublions pas la rigidité associée à la convection implicite
  1513. * calculée avec le nouveau terme de pénalisation
  1514. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4));
  1515. HNDT1 = HNDT1 'ET' RIGC ;
  1516. 'FINSI' ;
  1517. NBIT = NBIT + 1 ;
  1518. 'FINSI' ;
  1519. 'SI' (&BOU6 EGA ITMAXI) ;
  1520. 'MESS' 'Sortie sans convergence' ;
  1521. 'QUITTER' BOU6 ;
  1522. 'FINSI' ;
  1523. *-----------
  1524. 'FIN' BOU6 ;
  1525. *-----------
  1526.  
  1527. * Calcul de la dissolution :
  1528. DIS1 = ((NUM0 * PRCI) - (PRE2 * DEN0)) * VOLU1 ;
  1529.  
  1530. NB0 = &BOU6 - 1 ;
  1531. 'SI' (NB0 > 1);
  1532. 'MESS' ' ' NB0 ' itérations de solubilité' ;
  1533. 'FINSI';
  1534.  
  1535. 'MENA' ;
  1536. * Fin module PENALISATION
  1537. 'FINSI' ;
  1538.  
  1539. *|------------------------------------------------------|
  1540. *| Dissolution instantanée - Point fixe Dissolution |
  1541. *|------------------------------------------------------|
  1542. 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION')) ;
  1543.  
  1544. * On recalcule la matrice hybride si DeltaT a changé :
  1545. 'SI' Recalhyb ;
  1546. * Coefficients correspondant à la décroissance expl. et implicite :
  1547. NUM0 = (BETA - 1.) * LAMBD0 * DELTAT + 1. ;
  1548. DEN0 = LAMBD0 * BETA * DELTAT + 1. ;
  1549. * Calcul du coefficient devant (Ct+dt - Ct)/dt
  1550. TAB.'SURF' = ('KCHA' MODHYB PORETSU1 'CHAM') * DEN0 ;
  1551. * Matrice du problème
  1552. HNDTR = 'MATP' MODHYB MASHYB TAB ;
  1553. 'SI' ( 'EXISTE' TRANSI 'BLOCAGE' ) ;
  1554. HNDT1 = HNDTR 'ET' MATBLOC ;
  1555. 'SINON' ;
  1556. HNDT1 = HNDTR ;
  1557. 'FINSI' ;
  1558. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1559. * On ajoute la rigidité associée à la convection implicite
  1560. * Elle ne dépend ni du terme source, ni de la masse hybride,
  1561. * ni de la valeur de la trace de concentration. On doit néanmoins
  1562. * fournir des valeurs de ces grandeurs :
  1563. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCE CONV1 ;
  1564. HNDT1 = HNDT1 'ET' RIGC ;
  1565. 'FINSI' ;
  1566. NBIT = NBIT + 1 ;
  1567. Recalhyb = FAUX ;
  1568. 'FINSI' ;
  1569. *
  1570. *- Résolution en trace de concentration et calcul des inconnues secondaires
  1571. * Méthode du point fixe sur la dissolution (dissolution instantanée)
  1572.  
  1573. * Décroissance et dissolution du précipité avant les itérations
  1574. * La filiation du précipité est en principe comprise dans
  1575. * celle du soluté (à la charge de l'utilisateur)
  1576. PRE1 = PRCI * (NUM0 / DEN0) ;
  1577. * On prend pour commencer le T.S. de solubilité du pas précédent,
  1578. * limité au maximum précipité :
  1579. DIS0 = PRCI * VOLU1 * NUM0 ;
  1580. DIS1 = 0.5 * (DIS0 + DIS - ('ABS' (DIS0 - DIS)));
  1581.  
  1582. *--------------
  1583. 'REPETE' BOU7 ; 'COMM Boucle limite de sol., Picard dissolution' ;
  1584. *--------------
  1585. TERSC4 = DIS1 / DELTAT ;
  1586. TERSC4 = 'EXCO' TERSC4 'H' 'SOUR' ;
  1587. TERSCT = TERSCE + TERSC4 ;
  1588. PRE2 = PRE1 - (DIS1 / VOLU1 / DEN0) ;
  1589.  
  1590. * Calcul du second membre :
  1591. * SMTR0 contient la contribution en termes de flux des conditions aux
  1592. * limites. Dans tous les cas, SMTR1 contient le second membre du
  1593. * à la diffusion, les termes sources et le second membre du à
  1594. * la convection, mais TERSCT ne contient la convection que si
  1595. * TETAC est supérieur à 10-4
  1596. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1597. * TERSCT contient tous les T.S. sauf la convection, on
  1598. * regénére la rigidité au passage mais on ne s'en sert pas.
  1599. RIGC SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT CONV1 ;
  1600. 'SINON' ;
  1601. * TERSCT contient tous les T.S. convection comprise
  1602. SMTR1 = 'SMTP' MODHYB MASHYB TAB TERSCT ;
  1603. 'FINSI' ;
  1604. SMTR = SMTR0 + SMTR1 ;
  1605.  
  1606. *- Résolution
  1607. TP1 = 'RESO' HNDT1 SMTR ;
  1608. TP2 = 'EXCO' TP1 'TH' 'TH' ;
  1609.  
  1610. *- Terme source de convection pour calculer la nouvelle concentration
  1611. * TERSCVT contient tous les T.S. y compris la convection
  1612. 'SI' (( 'EXISTE' TRANSI 'CONVECTION' ) 'ET' (TETAC > 1.D-4)) ;
  1613. * On l'ajoute au terme source total maintenant seulement car on
  1614. * a besoin de C(t+dt).
  1615. CONC1 = ((1.D0 - TETAC) * TRP) + (TETAC * TP2) ;
  1616. CONC2 = (-1.D0) * ('NOMC' 'SCAL' CONC1) * CONV1 ;
  1617. CONC3 = 'DIVU' MODHYB CONC2 MCHYB ;
  1618. TERSCV = 'NOMC' 'SOUR' CONC3 ;
  1619. TERSCVT = TERSCT + TERSCV ;
  1620. 'SINON' ;
  1621. * Le terme source de convection est déjà dans TERSCT
  1622. TERSCVT = TERSCT ;
  1623. 'FINSI' ;
  1624. CHA2 = 'HYBP' MODHYB MASHYB TP2 TAB TERSCVT ;
  1625. DIS3 = DEN0 * (LIMSOL - CHA2) * ('EXCO' PORETSU1 'CK' 'SCAL') ;
  1626. * On limite la dissolution au maximum précipité :
  1627. DIS0 = PRE2 * VOLU1 * DEN0 ;
  1628. DIS4 = 0.5 * (DIS0 + DIS3 - ('ABS' (DIS0 - DIS3)));
  1629.  
  1630. * Critère de sortie de la boucle d'itération :
  1631. CRIT = (ABS DIS4) -
  1632. (DEN0 * EPS1 * ('EXCO' PORETSU1 'CK' 'SCAL') * LIMSOL) ;
  1633. 'SI' ( (MAXI CRIT) < 0.) ;
  1634. 'QUITTER' BOU7 ;
  1635. 'FINSI' ;
  1636. 'SI' (&BOU7 EGA ITMAXI) ;
  1637. 'MESS' 'Sortie sans convergence' ;
  1638. 'QUITTER' BOU7 ;
  1639. 'FINSI' ;
  1640.  
  1641. DIS1 = DIS1 + DIS4 ;
  1642.  
  1643. *-----------
  1644. 'FIN' BOU7 ;
  1645. *-----------
  1646.  
  1647. NB0 = &BOU7 - 1 ;
  1648. 'SI' (NB0 > 1);
  1649. 'MESS' ' ' NB0 ' itérations de solubilité' ;
  1650. 'FINSI';
  1651. NBIT = NBIT + NB0 ;
  1652.  
  1653. * Calcul des nouveaux flux :
  1654. FLU2 = 'HDEB' MODHYB MASHYB CHA2 TP2 ;
  1655.  
  1656. 'MENA' ;
  1657. * Fin module POINT FIXE DISSOLUTION
  1658. 'FINSI' ;
  1659.  
  1660. *|-------------------------------------|
  1661. *| Fin boucle de Retard Non Linéaire |
  1662. *|-------------------------------------|
  1663. *- Calcul du nouveau coefficient de retard et de l'adsorbât
  1664. 'SI' RNONLIN ;
  1665. * 1. F(C)
  1666. CC0 = 'NOMC' CHRG 'SCAL' ;
  1667. 'SI' RNONLINL ;
  1668. * Langmuir
  1669. FF = RETAR1M1 * CC0 / (1.D0 + ( RM1SURF * CC0 )) ;
  1670. 'SINON' ;
  1671. * Freundlich
  1672. * La bidouille utilisant les masques sert à obtenir le champ point
  1673. * des signes de CC0.
  1674. a = 'MASQUE' CC0 'SUPERIEUR' 0. ;
  1675. b = 'MASQUE' CC0 'INFERIEUR' 0. ;
  1676. FF = RETAR1 * (('ABS' CC0) ** UNSURN) * (a - b) ;
  1677. 'DETRUIT' a ;
  1678. 'DETRUIT' b ;
  1679. 'FINSI' ;
  1680. FF0 = FF ;
  1681. * 2. F(C + dC)
  1682. CC1 = ('NOMC' CHA2 'SCAL') + EpsCord ;
  1683. 'SI' RNONLINL ;
  1684. * Langmuir
  1685. FF1 = RETAR1M1 * CC1 / (1.D0 + ( RM1SURF * CC1 )) ;
  1686. 'SINON' ;
  1687. * Freundlich
  1688. a = 'MASQUE' CC1 'SUPERIEUR' 0. ;
  1689. b = 'MASQUE' CC1 'INFERIEUR' 0. ;
  1690. FF1 = RETAR1 * (('ABS' CC1) ** UNSURN) * (a - b) ;
  1691. 'DETRUIT' a ;
  1692. 'DETRUIT' b ;
  1693. 'FINSI' ;
  1694. * 3. R = 1 + dF/dC
  1695. RETAR2 = 1.D0 + ( (FF1 - FF0) / (CC1 - CC0) ) ;
  1696. 'DETRUIT' CC1 ;
  1697. 'DETRUIT' FF1 ;
  1698. *
  1699. *- Critère de sortie de la boucle de retard non linéaire :
  1700. CRIT = ( 'ABS' (RETAR2 - RETARC) ) / RETARC ;
  1701. CRIT = 'MAXIMUM' CRIT ;
  1702. 'SI' (CRIT < EPSRNL) ;
  1703. 'QUITTER' BOURNL ;
  1704. 'SINON' ;
  1705. * Le coef. de retard ne change que si le critère est mauvais.
  1706. * Réévaluation des variables affectées :
  1707. RETARC = RETAR2 ;
  1708. PORETSU1 = PORSUF1 * RETARC ;
  1709. Recalhyb = VRAI ;
  1710. 'FINSI' ;
  1711. 'SINON' ;
  1712. * Retard linéaire F = (R - 1) C
  1713. CC0 = 'NOMC' CHRG 'SCAL' ;
  1714. FF = RETAR1M1 * CC0 ;
  1715. 'FINSI' ;
  1716.  
  1717. * Si on atteint le nombre maximum d'itérations - en particulier si le
  1718. * retard est linéaire - on quitte la boucle :
  1719. 'SI' (IRNL 'EGA' ITMAXRNL) ;
  1720. ' sortie boucle de retard non linéaire sans convergence' ;
  1721. 'QUITTER' BOURNL ;
  1722. 'FINSI' ;
  1723.  
  1724. *-------------
  1725. 'FIN' BOURNL ;
  1726. *-------------
  1727.  
  1728. 'SI' (IRNL > 1 ) ;
  1729. 'MESS' ' ' IRNL ' itérations de retard non linéaire.'
  1730. ' Critère = ' crit ;
  1731. 'FINSI' ;
  1732. NBITRNL = NBITRNL + (IRNL - 1) ;
  1733.  
  1734.  
  1735. *|-----------------------------------|
  1736. *| Fin des modules de résolution |
  1737. *|-----------------------------------|
  1738. *
  1739. *- Calcul de la dissolution par unité de volume et par unité de temps :
  1740. *
  1741. 'SI' SOLUB ;
  1742. DIS2 = DIS1 / VOLU1 / DELTAT ;
  1743. 'FINSI' ;
  1744. *
  1745. *- Archivage des resultats
  1746. *
  1747. * si le prochain temps à sauver tombe entre les deux
  1748. * incréments de temps, on obtient les valeurs à sauver
  1749. * par interpolation linéaire puis on les stocke.
  1750. *
  1751. 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ;
  1752. * Sauvegarde de tous les temps intermédiaires, s'il y en a :
  1753. 'REPETER' BOU9 ;
  1754. TEMS = 'EXTR' TPSOR ISOR ;
  1755. 'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '&lt;EG' TPS ) ) ;
  1756. LAST1 = LAST1 + 1 ;
  1757. DTEM = ( TPS - TEMS ) / DELTAT ;
  1758. UNMO = 1. - DTEM ;
  1759. TPBIS = 'COLI' TRP DTEM TP2 UNMO ;
  1760. CHBIS = 'COLI' CHRG DTEM CHA2 UNMO ;
  1761. FLUBIS = 'COLI' FLU0 DTEM FLU2 UNMO ;
  1762. 'SI' ('EGA' ('TYPE' RETAR0) 'CHPOINT ') ;
  1763. RETBIS = 'COLI' RETAR0 DTEM RETARC UNMO ;
  1764. 'SINON' ;
  1765. RETBIS = (DTEM * RETAR0) + (UNMO * RETARC) ;
  1766. 'FINSI' ;
  1767. 'SI' SOLUB ;
  1768. PRBIS = 'COLI' PRCI DTEM PRE2 UNMO ;
  1769. DISBIS = 'COLI' (DIS / VOLU1 / DELTAT) DTEM
  1770. DIS2 UNMO ;
  1771. 'FINSI' ;
  1772. TRANSI . 'TEMPS' . LAST1 = TEMS ;
  1773. TRANSI . NOMTINC . LAST1 = TPBIS ;
  1774. TRANSI . NOMINC . LAST1 = CHBIS ;
  1775. TRANSI . 'FLUX' . LAST1 = FLUBIS ;
  1776. 'SI' SOLUB ;
  1777. TRANSI . 'PRECIPITE' . LAST1 = PRBIS ;
  1778. TRANSI . 'DISSOLUTION' . LAST1 = DISBIS ;
  1779. 'FINSI' ;
  1780. 'SI' SAUVRET ;
  1781. TRANSI . 'RETARD' . LAST1 = RETBIS ;
  1782. 'FINSI' ;
  1783. ISOR = ISOR + 1 ;
  1784. 'SINON' ;
  1785. 'QUITTER' BOU9 ;
  1786. 'FINSI' ;
  1787. 'SI' ( ISOR '>' DSOR ) ;
  1788. ISOR = DSOR ;
  1789. 'QUITTER' BOU9 ;
  1790. 'FINSI' ;
  1791. 'FIN' BOU9 ;
  1792. 'SINON' ;
  1793. * sinon, on sauvegarde tous les temps :
  1794. LAST1 = LAST1 + 1 ;
  1795. TRANSI . 'TEMPS' . LAST1 = TPS ;
  1796. TRANSI . NOMTINC . LAST1 = TP2 ;
  1797. TRANSI . NOMINC . LAST1 = CHA2 ;
  1798. TRANSI . 'FLUX' . LAST1 = FLU2 ;
  1799. 'SI' SOLUB ;
  1800. TRANSI . 'PRECIPITE' . LAST1 = PRE2 ;
  1801. TRANSI . 'DISSOLUTION' . LAST1 = DIS2 ;
  1802. 'FINSI' ;
  1803. 'SI' SAUVRET ;
  1804. TRANSI . 'RETARD' . LAST1 = RETARC ;
  1805. 'FINSI' ;
  1806. 'FINSI' ;
  1807. * Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait.
  1808. * Le dernier pas de temps ne peut pas faire l'objet d'interpolation
  1809. * linéaire si le temps sauvegardé ne tombe pas pile dessus.
  1810. 'SI' (IPAS 'EGA' DCAL) ;
  1811. TPSDER = TRANSI. 'TEMPS' . LAST1 ;
  1812. 'SI' ( TPSDER 'NEG' TPSFIN EPSDT ) ;
  1813. LAST1 = LAST1 + 1 ;
  1814. TRANSI . 'TEMPS' . LAST1 = TPS ;
  1815. TRANSI . NOMTINC . LAST1 = TP2 ;
  1816. TRANSI . NOMINC . LAST1 = CHA2 ;
  1817. TRANSI . 'FLUX' . LAST1 = FLU2 ;
  1818. 'SI' SOLUB ;
  1819. TRANSI . 'PRECIPITE' . LAST1 = PRE2 ;
  1820. TRANSI . 'DISSOLUTION' . LAST1 = DIS2 ;
  1821. 'FINSI' ;
  1822. 'SI' SAUVRET ;
  1823. TRANSI . 'RETARD' . LAST1 = RETARC ;
  1824. 'FINSI' ;
  1825. 'FINSI' ;
  1826. 'FINSI' ;
  1827. *
  1828. *- Initialisations pour le pas suivant
  1829. *
  1830. PRECED = TPS ;
  1831. TRP = TP2 ;
  1832. CHRG = CHA2 ;
  1833. FLU0 = FLU2 ;
  1834. RETAR0 = RETARC ;
  1835. DELOLD = DELTAT ;
  1836. 'SI' SOLUB ;
  1837. DIS = DIS1 ;
  1838. PRCI = PRE2 ;
  1839. 'FINSI' ;
  1840.  
  1841. 'MENAGE' ;
  1842. *
  1843. *=============
  1844. 'FIN' BOUTPS ; COMM 'Boucle sur le temps' ;
  1845. *=============
  1846. *
  1847. *- Indication du nombre de passages :
  1848. *
  1849. 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION')) ;
  1850. 'MESS' 'Nombre total de matrices calculées : ' NBIT ;
  1851. 'FINSI' ;
  1852. 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION')) ;
  1853. 'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
  1854. NBIT ;
  1855. 'FINSI' ;
  1856. 'SI' SOLUP ;
  1857. 'MESS' 'Nombre total d itérations supplémentaires de solubilité : '
  1858. NBIT;
  1859. 'FINSI' ;
  1860. 'SI' RNONLIN ;
  1861. 'MESS' 'Nombre total d itérations supplémentaires de '
  1862. 'retard non linéaire : ' NBITRNL;
  1863. 'FINSI' ;
  1864.  
  1865. 'FINP' ;
  1866.  

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