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 <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 '<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 <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 '<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 '<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 '<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 : MESS 1931 : ' sortie boucle de retard non linéaire sans convergence' ; 1932 : 'QUITTER' BOURNL ; 1933 : 'FINSI' ; 1934 : 1935 : *------------- 1936 : 'FIN' BOURNL ; 1937 : *------------- 1938 : 1939 : 'SI' (IRNL > 1 ) ; 1940 : 'MESS' ' ' IRNL ' itérations de retard non linéaire.' 1941 : ' Critère = ' crit ; 1942 : 'FINSI' ; 1943 : NBITRNL = NBITRNL + (IRNL - 1) ; 1944 : 1945 : 1946 : *|-----------------------------------| 1947 : *| Fin des modules de résolution | 1948 : *|-----------------------------------| 1949 : * 1950 : *- Calcul de la dissolution par unité de volume et par unité de temps : 1951 : * 1952 : 'SI' SOLUA ; 1953 : DIS2 = DIS1 / VOLU1 / DELTAT ; 1954 : 'FINSI' ; 1955 : * 1956 : *- Archivage des resultats 1957 : * 1958 : * si le prochain temps à sauver tombe entre les deux 1959 : * incréments de temps, on obtient les valeurs à sauver 1960 : * par interpolation linéaire puis on les stocke. 1961 : * 1962 : 1963 : *list ((PLACE) * 4. / 1000000.) ; 1964 : 1965 : 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ; 1966 : 1967 : * Sauvegarde de tous les temps intermédiaires, s'il y en a : 1968 : 'REPETER' BOU9 ; 1969 : TEMS = 'EXTR' TPSOR ISOR ; 1970 : 1971 : 'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '<EG' TPS ) ) ; 1972 : * on interpole les variables linéairement entre 1973 : * les pas de temps 1974 : LAST1 = LAST1 + 1 ; 1975 : DTEM = ( TPS - TEMS ) / DELTAT ; 1976 : UNMO = 1.D0 - DTEM ; 1977 : CHBIS = 'COLI' CHRG DTEM CHA2 UNMO ; 1978 : FLUBIS = 'COLI' FLU0 DTEM FLU2 UNMO ; 1979 : FLUBISCO= 'COLI' FLUCO0 DTEM FLUCO2 UNMO ; 1980 : 'SI' ('EGA' ('TYPE' RETAR0) 'CHPOINT ') ; 1981 : RETBIS = 'COLI' RETAR0 DTEM RETARC UNMO ; 1982 : 'SINON' ; 1983 : RETBIS = (DTEM * RETAR0) + (UNMO * RETARC) ; 1984 : 'FINSI' ; 1985 : 'SI' SOLUB ; 1986 : PRBIS = 'COLI' PRCI DTEM PRE2 UNMO ; 1987 : 'FINSI' ; 1988 : TRANSI . 'TEMPS' . LAST1 = TEMS ; 1989 : TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP 1990 : CHBIS ; 1991 : 'DETRUIT' CHBIS ; 1992 : TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP 1993 : FLUBIS ; 1994 : 'DETRUIT' FLUBIS ; 1995 : TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP 1996 : FLUBISCO ; 1997 : 'DETRUIT' FLUBISCO ; 1998 : 'SI' SOLUB ; 1999 : TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP 2000 : PRBIS ; 2001 : 'DETRUIT' PRBIS ; 2002 : 'FINSI' ; 2003 : 'SI' SAUVRET ; 2004 : TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP 2005 : RETBIS ; 2006 : 'DETRUIT' RETBIS ; 2007 : 'FINSI' ; 2008 : ISOR = ISOR + 1 ; 2009 : 'SINON' ; 2010 : 'QUITTER' BOU9 ; 2011 : 'FINSI' ; 2012 : 'SI' ( ISOR '>' DSOR ) ; 2013 : ISOR = DSOR ; 2014 : 'QUITTER' BOU9 ; 2015 : 'FINSI' ; 2016 : 'FIN' BOU9 ; 2017 : 'SINON' ; 2018 : * sinon, on sauvegarde tous les temps : 2019 : LAST1 = LAST1 + 1 ; 2020 : TRANSI . 'TEMPS' . LAST1 = TPS ; 2021 : TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP CHA2 ; 2022 : TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP FLU2 ; 2023 : TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP FLUCO2 ; 2024 : 'SI' SOLUB ; 2025 : TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP PRE2 ; 2026 : 'FINSI' ; 2027 : 'SI' SAUVRET ; 2028 : TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP RETARC ; 2029 : 'FINSI' ; 2030 : 'FINSI' ; 2031 : 2032 : 2033 : *On sauve l'intégrale de decroissance*concentration et 2034 : * decroissance*précipité au cours du temps 2035 : * stockée pour l'instant à tous les temps. Peut s'améliorer. 2036 : 2037 : * indice pour le stockage des intégrales de concentration 2038 : LASTINT = LASTINT '+' 1; 2039 : 2040 : BM1 = 1.D0 '-' BETA ; 2041 : * si l'espece fait l'objet d'une filiation on stocke les intégrales 2042 : * du précipité et de la concentration. 2043 : SI ((EXISTE TRANSI DECROISSANCE) 'ET' (TRANSI . 'PERE')) ; 2044 : *TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1) 2045 : * '+' (DELTAT * PORETSU * TRANSI . 'DECROISSANCE' * 2046 : * ((BM1 * ('NOMC' NOMESP CHRG)) 2047 : * '+' ((BETA) * ('NOMC' NOMESP CHA2)))) ; 2048 : zozo = (DELTAT * TRANSI . 'DECROISSANCE') * PORETSU ; 2049 : zaza = BM1 * CHRG ; 2050 : zuzu = BETA * CHA2 ; 2051 : zeze = 'KOPS' zaza '+' zuzu ; 2052 : 'DETRUIT' zaza ; 2053 : 'DETRUIT' zuzu ; 2054 : zaza = 'KOPS' zozo * zeze ; 2055 : 'DETRUIT' zozo ; 2056 : 'DETRUIT' zeze ; 2057 : zozo = 'NOMC' NOMESP zaza ; 2058 : 'DETRUIT' zaza ; 2059 : TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . (LASTINT '-' 1) 2060 : '+' zozo ; 2061 : 'DETRUIT' zozo ; 2062 : 2063 : 'SI' SOLUB ; 2064 : 2065 : * TRANSI . 'INTFILI' . LASTINT = TRANSI . 'INTFILI' . LASTINT 2066 : * '+' (DELTAT * (1.D0 '-' POROS) * TRANSI . 'DECROISSANCE' * 2067 : * ((BM1 * ('NOMC' NOMESP PRCI)) 2068 : * '+' ((BETA) * ('NOMC' NOMESP PRE2)))) ; 2069 : zizi = (1.D0 '-' POROS) ; 2070 : zozo = (DELTAT * TRANSI . 'DECROISSANCE') * zizi ; 2071 : 'DETRUIT' zizi ; 2072 : zaza = BM1 * PRCI ; 2073 : zuzu = BETA * PRE2 ; 2074 : zeze = 'KOPS' zaza '+' zuzu ; 2075 : 'DETRUIT' zaza ; 2076 : 'DETRUIT' zuzu ; 2077 : zaza = 'KOPS' zozo * zeze ; 2078 : 'DETRUIT' zozo ; 2079 : 'DETRUIT' zeze ; 2080 : zozo = 'NOMC' NOMESP zaza ; 2081 : 'DETRUIT' zaza ; 2082 : zaza = TRANSI . 'INTFILI' . LASTINT '+' zozo ; 2083 : 'DETRUIT' TRANSI . 'INTFILI' . LASTINT ; 2084 : 'DETRUIT' zozo ; 2085 : TRANSI . 'INTFILI' . LASTINT = zaza ; 2086 : 'OUBLIER' zaza ; 2087 : 2088 : 'FINSI' ; 2089 : FINSI ; 2090 : 2091 : 2092 : * Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait. 2093 : * Le dernier pas de temps ne peut pas faire l'objet d'interpolation 2094 : * linéaire si le temps sauvegardé ne tombe pas pile dessus. 2095 : 'SI' (IPAS 'EGA' DCAL) ; 2096 : TPSDER = TRANSI. 'TEMPS' . LAST1 ; 2097 : 'SI' ( TPSDER 'NEG' TPSFIN EPSDT ) ; 2098 : LAST1 = LAST1 + 1 ; 2099 : TRANSI . 'TEMPS' . LAST1 = TPS ; 2100 : TRANSI . NOMINC . LAST1 = 'NOMC' NOMESP CHA2 ; 2101 : TRANSI . 'FLUXDIFF' . LAST1 = 'NOMC' NOMESP FLU2 ; 2102 : TRANSI . 'FLUXCONV' . LAST1 = 'NOMC' NOMESP FLUCO2 ; 2103 : 'SI' SOLUB ; 2104 : TRANSI . 'PRECIPITE' . LAST1 = 'NOMC' NOMESP PRE2 ; 2105 : 'FINSI' ; 2106 : 'SI' SAUVRET ; 2107 : TRANSI . 'RETARD' . LAST1 = 'NOMC' NOMESP RETARC ; 2108 : 'FINSI' ; 2109 : 'FINSI' ; 2110 : 'FINSI' ; 2111 : * 2112 : *- Initialisations pour le pas suivant 2113 : * 2114 : PRECED = TPS ; 2115 : 'DETRUIT' CHRG ; 2116 : CHRG = CHA2 ; 2117 : 'DETRUIT' FLU0 ; 2118 : FLU0 = FLU2 ; 2119 : 'DETRUIT' FLUCO0 ; 2120 : FLUCO0 = FLUCO2 ; 2121 : 'OUBLIER' GEOL1 ; 2122 : 'OUBLIER' GEOL2 ; 2123 : GEOL1 = 'TABLE' GEOLPF1 ; 2124 : GEOL2 = 'TABLE' GEOLPF2 ; 2125 : 'DETRUIT' RETAR0 ; 2126 : RETAR0 = RETARC ; 2127 : DELOLD = DELTAT ; 2128 : 'SI' SOLUB ; 2129 : 'DETRUIT' PRCI ; 2130 : PRCI = PRE2 ; 2131 : 'FINSI' ; 2132 : 'SI' ('EXISTE' TRANSI 'SOURCE') ; 2133 : 'DETRUIT' TERSC2M1 ; 2134 : TERSC2M1 = TERSC2 ; 2135 : 'FINSI' ; 2136 : 'SI' ('EXISTE' TRANSI 'FLUX_IMPOSE') ; 2137 : 'DETRUIT' FLUIMPM1 ; 2138 : FLUIMPM1 = FLUIMPO ; 2139 : 'FINSI' ; 2140 : 'SI' ('EXISTE' TRANSI 'FLUXTOT_IMP') ; 2141 : 'DETRUIT' FLUTMPM1 ; 2142 : FLUTMPM1 = FLUTMPO ; 2143 : 'FINSI' ; 2144 : 'SI' ('EXISTE' TRANSI 'MIXTES') ; 2145 : 'DETRUIT' FLUMMPM1 ; 2146 : FLUMMPM1 = FLUMMPO ; 2147 : 'FINSI' ; 2148 : 2149 : *MESS 'TEMPS2' ; TEMPS; 2150 : *'MENAGE' ; 2151 : 2152 : * destruction d'objets temporaires 2153 : 2154 : 'DETRUIT' TERSC3 ; 2155 : 2156 : *list ((PLACE) * 4. / 1000000.) ; 2157 : 2158 : *temps impr; 2159 : 2160 : * 2161 : *============= 2162 : 'FIN' BOUTPS ; COMM 'Boucle sur le temps' ; 2163 : *============= 2164 : * 2165 : *- Indication du nombre de passages : 2166 : * 2167 : 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION')) ; 2168 : 'MESS' 'Nombre total de matrices calculées : ' NBIT ; 2169 : 'FINSI' ; 2170 : 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION')) ; 2171 : 'MESS' 'Nombre total d itérations supplémentaires de solubilité : ' 2172 : NBIT ; 2173 : 'FINSI' ; 2174 : 'SI' SOLUP ; 2175 : 'MESS' 'Nombre total d itérations supplémentaires de solubilité : ' 2176 : NBIT; 2177 : 'FINSI' ; 2178 : 'SI' RNONLIN ; 2179 : 'MESS' 'Nombre total d itérations supplémentaires de ' 2180 : 'retard non linéaire : ' NBITRNL; 2181 : 'FINSI' ; 2182 : 2183 : 'FINP' ; 2184 : 2185 : 2186 :
© Cast3M 2003 - All rights reserved.
Disclaimer