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 <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 '<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 : MESS 1721 : ' sortie boucle de retard non linéaire sans convergence' ; 1722 : 'QUITTER' BOURNL ; 1723 : 'FINSI' ; 1724 : 1725 : *------------- 1726 : 'FIN' BOURNL ; 1727 : *------------- 1728 : 1729 : 'SI' (IRNL > 1 ) ; 1730 : 'MESS' ' ' IRNL ' itérations de retard non linéaire.' 1731 : ' Critère = ' crit ; 1732 : 'FINSI' ; 1733 : NBITRNL = NBITRNL + (IRNL - 1) ; 1734 : 1735 : 1736 : *|-----------------------------------| 1737 : *| Fin des modules de résolution | 1738 : *|-----------------------------------| 1739 : * 1740 : *- Calcul de la dissolution par unité de volume et par unité de temps : 1741 : * 1742 : 'SI' SOLUB ; 1743 : DIS2 = DIS1 / VOLU1 / DELTAT ; 1744 : 'FINSI' ; 1745 : * 1746 : *- Archivage des resultats 1747 : * 1748 : * si le prochain temps à sauver tombe entre les deux 1749 : * incréments de temps, on obtient les valeurs à sauver 1750 : * par interpolation linéaire puis on les stocke. 1751 : * 1752 : 'SI' ( 'EXISTE' TRANSI 'TEMPS_SAUVES' ) ; 1753 : * Sauvegarde de tous les temps intermédiaires, s'il y en a : 1754 : 'REPETER' BOU9 ; 1755 : TEMS = 'EXTR' TPSOR ISOR ; 1756 : 'SI' ( ( PRECED '<' TEMS ) 'ET' ( TEMS '<EG' TPS ) ) ; 1757 : LAST1 = LAST1 + 1 ; 1758 : DTEM = ( TPS - TEMS ) / DELTAT ; 1759 : UNMO = 1. - DTEM ; 1760 : TPBIS = 'COLI' TRP DTEM TP2 UNMO ; 1761 : CHBIS = 'COLI' CHRG DTEM CHA2 UNMO ; 1762 : FLUBIS = 'COLI' FLU0 DTEM FLU2 UNMO ; 1763 : 'SI' ('EGA' ('TYPE' RETAR0) 'CHPOINT ') ; 1764 : RETBIS = 'COLI' RETAR0 DTEM RETARC UNMO ; 1765 : 'SINON' ; 1766 : RETBIS = (DTEM * RETAR0) + (UNMO * RETARC) ; 1767 : 'FINSI' ; 1768 : 'SI' SOLUB ; 1769 : PRBIS = 'COLI' PRCI DTEM PRE2 UNMO ; 1770 : DISBIS = 'COLI' (DIS / VOLU1 / DELTAT) DTEM 1771 : DIS2 UNMO ; 1772 : 'FINSI' ; 1773 : TRANSI . 'TEMPS' . LAST1 = TEMS ; 1774 : TRANSI . NOMTINC . LAST1 = TPBIS ; 1775 : TRANSI . NOMINC . LAST1 = CHBIS ; 1776 : TRANSI . 'FLUX' . LAST1 = FLUBIS ; 1777 : 'SI' SOLUB ; 1778 : TRANSI . 'PRECIPITE' . LAST1 = PRBIS ; 1779 : TRANSI . 'DISSOLUTION' . LAST1 = DISBIS ; 1780 : 'FINSI' ; 1781 : 'SI' SAUVRET ; 1782 : TRANSI . 'RETARD' . LAST1 = RETBIS ; 1783 : 'FINSI' ; 1784 : ISOR = ISOR + 1 ; 1785 : 'SINON' ; 1786 : 'QUITTER' BOU9 ; 1787 : 'FINSI' ; 1788 : 'SI' ( ISOR '>' DSOR ) ; 1789 : ISOR = DSOR ; 1790 : 'QUITTER' BOU9 ; 1791 : 'FINSI' ; 1792 : 'FIN' BOU9 ; 1793 : 'SINON' ; 1794 : * sinon, on sauvegarde tous les temps : 1795 : LAST1 = LAST1 + 1 ; 1796 : TRANSI . 'TEMPS' . LAST1 = TPS ; 1797 : TRANSI . NOMTINC . LAST1 = TP2 ; 1798 : TRANSI . NOMINC . LAST1 = CHA2 ; 1799 : TRANSI . 'FLUX' . LAST1 = FLU2 ; 1800 : 'SI' SOLUB ; 1801 : TRANSI . 'PRECIPITE' . LAST1 = PRE2 ; 1802 : TRANSI . 'DISSOLUTION' . LAST1 = DIS2 ; 1803 : 'FINSI' ; 1804 : 'SI' SAUVRET ; 1805 : TRANSI . 'RETARD' . LAST1 = RETARC ; 1806 : 'FINSI' ; 1807 : 'FINSI' ; 1808 : * Sauvegarde du dernier pas de temps si ça n'a pas déjà été fait. 1809 : * Le dernier pas de temps ne peut pas faire l'objet d'interpolation 1810 : * linéaire si le temps sauvegardé ne tombe pas pile dessus. 1811 : 'SI' (IPAS 'EGA' DCAL) ; 1812 : TPSDER = TRANSI. 'TEMPS' . LAST1 ; 1813 : 'SI' ( TPSDER 'NEG' TPSFIN EPSDT ) ; 1814 : LAST1 = LAST1 + 1 ; 1815 : TRANSI . 'TEMPS' . LAST1 = TPS ; 1816 : TRANSI . NOMTINC . LAST1 = TP2 ; 1817 : TRANSI . NOMINC . LAST1 = CHA2 ; 1818 : TRANSI . 'FLUX' . LAST1 = FLU2 ; 1819 : 'SI' SOLUB ; 1820 : TRANSI . 'PRECIPITE' . LAST1 = PRE2 ; 1821 : TRANSI . 'DISSOLUTION' . LAST1 = DIS2 ; 1822 : 'FINSI' ; 1823 : 'SI' SAUVRET ; 1824 : TRANSI . 'RETARD' . LAST1 = RETARC ; 1825 : 'FINSI' ; 1826 : 'FINSI' ; 1827 : 'FINSI' ; 1828 : * 1829 : *- Initialisations pour le pas suivant 1830 : * 1831 : PRECED = TPS ; 1832 : TRP = TP2 ; 1833 : CHRG = CHA2 ; 1834 : FLU0 = FLU2 ; 1835 : RETAR0 = RETARC ; 1836 : DELOLD = DELTAT ; 1837 : 'SI' SOLUB ; 1838 : DIS = DIS1 ; 1839 : PRCI = PRE2 ; 1840 : 'FINSI' ; 1841 : 1842 : 'MENAGE' ; 1843 : * 1844 : *============= 1845 : 'FIN' BOUTPS ; COMM 'Boucle sur le temps' ; 1846 : *============= 1847 : * 1848 : *- Indication du nombre de passages : 1849 : * 1850 : 'SI' (SOLUI 'ET' ('EGA' MET0 'PENALISATION')) ; 1851 : 'MESS' 'Nombre total de matrices calculées : ' NBIT ; 1852 : 'FINSI' ; 1853 : 'SI' (SOLUI 'ET' ('EGA' MET0 'PF DISSOLUTION')) ; 1854 : 'MESS' 'Nombre total d itérations supplémentaires de solubilité : ' 1855 : NBIT ; 1856 : 'FINSI' ; 1857 : 'SI' SOLUP ; 1858 : 'MESS' 'Nombre total d itérations supplémentaires de solubilité : ' 1859 : NBIT; 1860 : 'FINSI' ; 1861 : 'SI' RNONLIN ; 1862 : 'MESS' 'Nombre total d itérations supplémentaires de ' 1863 : 'retard non linéaire : ' NBITRNL; 1864 : 'FINSI' ; 1865 : 1866 : 'FINP' ; 1867 : 1868 : 1869 : 1870 : 1871 : 1872 :
© Cast3M 2003 - All rights reserved.
Disclaimer