Afficher cette notice en
1 : $$$$ DARCYSAT NOTICE CHAT 11/09/12 21:15:42 7124 2 : DATE 11/09/12 3 : Procedure DARCYSAT 4 : ------------------ Voir aussi : HT_PRO 5 : KR_PRO 6 : PECHE 7 : TRACHIS 8 : TRACHIT 9 : 10 : TAB1.'SOUSTYPE'.'MODELE' 11 : .'LOI_PERMEABILITE'.'LOI_SATURATION'.'HOMOGENEISATION' 12 : .'COEFEMMA'.'BLOCAGES_DARCY'.'CHARGEMENT'.'FORCE_GRAVITE' 13 : .'CONVERSION_CHARGE'.'TEMPS'.'TRACE_CHARGE'.'CHARGE'.('FLUX') 14 : .'TEMPS_CALCULES'.'TEMPS_SAUVES'.'TEMPS_FINAL' 15 : .'DT_INITIAL'.'NPAS'.'CFL' 16 : .'NITER'.'ITMAX'.'RESIDU_MAX'.'SOUS_RELAXATION'.'XI' 17 : .'DIVISION_DT'.'COFDIV'.'NMAXDT'.'MESSAGE' 18 : 19 : 20 : Objet : 21 : _______ 22 : 23 : Cette procedure permet de simuler un transitoire d'ecoulement 24 : en zones saturee et non saturee d'un milieu poreux. L'ecoulement 25 : est decrit en zone saturee par l'equation de Darcy et en zone non 26 : saturee par l'equation de Richard's (ecrite ici avec h en metres) 27 : 28 : zone saturee (Pw - Pg > 0) 29 : 30 : Se dh/dt = -div(U) ; U = -Ks grad(h) ; 31 : h = (Pw-Pg)/rho*g + z (m) ou h = (Pw-Pg) + rho*g*z (Pa) 32 : avec, Se, coefficient d'emmagasinement (1/m), 33 : Ks, permeabilite a saturation (m/s), 34 : Pw, pression d'eau (Pa), 35 : Pg, pression de gaz (Pa), 36 : h, charge (m) et 37 : U, vitesse de Darcy (m/s). 38 : 39 : zone non saturee (Pw - Pg < 0) 40 : 41 : C(P) dh/dt = -div(U) ; U = -K(P) (grad h) 42 : 43 : avec, C(P) capacite capillaire (1/m) 44 : K(P), permeabilite (m/s) 45 : 46 : La resolution est effectuee en trace de charge par la methode EFMH. 47 : La teneur en eau se deduit de la pression par la procedure HT_PRO. 48 : La permeabilite se deduit de la saturation par la procedure KR_PRO. 49 : La capacite capillaire est deduite de la teneur en eau et de la 50 : pression. 51 : 52 : Le probleme etant non lineaire, des mecanismes automatiques ou non 53 : ont ete mis en place pour fiabiliser la convergence du systeme : 54 : relaxation de la capacite et de la conductivite hydraulique et 55 : homogeneisation des caracteristiques decentree sur les elements, 56 : changement automatique de pas de temps ou penalisation sur le 57 : terme transitoire. 58 : 59 : Commentaire : 60 : _____________ 61 : 62 : En entree, TAB1 sert a definir les options et les parametres du 63 : calcul. Les indices de la table TAB1 sont des mots (a coder tels 64 : quels) dont voici la description : 65 : 66 : ___________________________________________________________________ 67 : | | 68 : | Indice Definition | 69 : | | 70 : ------------------------------------------------------------------- 71 : | | 72 : |------------------------------------------------- | 73 : | Donnees physiques, geometriques et materielles : | 74 : |------------------------------------------------- | 75 : | | 76 : |'SOUSTYPE' 'DARCY_TRANSATUR' (type MOT) | 77 : | | 78 : |'MODELE' Objet modele (MMODEL cree par MODE) | 79 : | | 80 : |'LOI_PERMEABILITE' Objet table contenant les valeurs des | 81 : | parametres de la loi de permeabilite calculee par la| 82 : | procedure KR_PRO. | 83 : | 3 indices 'SOUSTYPE' possibles : | 84 : | | 85 : | - soustypes au choix | 86 : | 'PUISSANCE', 'MUALEM', 'BURDINE', 'MUALEM_BURDINE' | 87 : | 'BROOKS_COREY', 'EXPONENTIELLE', 'LOGARITHMIQUE' | 88 : | : la table doit contenir les | 89 : | indices designant les parametres utilises dans la | 90 : | procedure standard KR_PRO. | 91 : | | 92 : | - soustype 'PERSONNELLE' : la table doit contenir | 93 : | les indices de la loi construite par l'utilisateur | 94 : | dans sa propre procedure KR_PRO, ou dans une autre | 95 : | dont le nom est stipule au sous-indice | 96 : | 'NOM_PROCEDURE'. | 97 : | | 98 : | - soustype 'MULTIZONE' : les indices doivent | 99 : | designer les differentes zones, les valeurs de ces | 100 : | indices etant alors des tables de soustype | 101 : | 'PUISSANCE', 'MUALEM', 'BURDINE', 'MUALEM_BURDINE' | 102 : | 'BROOKS_COREY', 'EXPONENTIELLE', 'LOGARITHMIQUE' | 103 : | documentees dans la notice de KR_PRO | 104 : | ou'PERSONNELLE' dont les indices | 105 : | designent les parametres des lois affectees a chaque| 106 : | zone. Ces sous-tables doivent aussi contenir le | 107 : | modele de la zone a l'indice 'MODELE', ces sous- | 108 : | modeles constituant une partition du domaine general| 109 : | | 110 : | pour les autres sous-indices (defaut), voir KR_PRO | 111 : | | 112 : |'LOI_SATURATION' Objet table contenant les valeurs des | 113 : | parametres de la loi analytique de teneur en eau | 114 : | calculee par la procedure HT_PRO. | 115 : | 3 indices 'SOUSTYPE' possibles : | 116 : | | 117 : | - soustypes aux choix 'VAN_GENUCHTEN', | 118 : | 'EXPONENTIELLE', 'LOGARITHMIQUE' | 119 : | : la table doit contenir les | 120 : | indices designant les parametres utilises dans la | 121 : | procedure standard HT_PRO. | 122 : | | 123 : | - soustype 'PERSONNELLE' : la table doit contenir | 124 : | les indices de la loi construite par l'utilisateur | 125 : | dans sa propre procedure HT_PRO, ou dans une autre | 126 : | dont le nom est stipule au sous-indice | 127 : | 'NOM_PROCEDURE'. | 128 : | | 129 : | - soustype 'MULTIZONE' : les indices doivent | 130 : | designer les differentes zones, les valeurs de ces | 131 : | indices etant alors des tables de soustype | 132 : | 'VAN_GENUCHTEN', 'EXPONENTIELLE', 'LOGARITHMIQUE' | 133 : | ou 'PERSONNELLE' dont les indices designent les | 134 : | parametres des lois affectees a chaque zone. | 135 : | Ces sous-tables doivent aussi contenir le modele de | 136 : | la zone a l'indice 'MODELE', ces sous-modeles | 137 : | constituant une partition du domaine general | 138 : | | 139 : | pour les autres sous-indices (defaut), voir HT_PRO | 140 : | | 141 : |'COEFEMMA' coefficient d'emmagasinement en zone saturee (m-1) | 142 : | (FLOTTANT ou CHPOINT, comp. 'SCAL', nature | 143 : | discrete, defaut = 0) | 144 : | | 145 : |'FORCE_GRAVITE' (type FLOTTANT, par ex. = 1000*9.81) valeur de la| 146 : | force de gravite (N/m3) egale a la masse volumique | 147 : | de l'eau multipliee par l'acceleration de la | 148 : | pesanteur. En absence de ce mot cle, la gravite | 149 : | n'est pas prise en compte. | 150 : | | 151 : |'CONVERSION_CHARGE' (type 'FLOTTANT', defaut = 1.) facteur de | 152 : | convertion de la charge en Pascals. | 153 : | Pour une pression qui s'exprime en Pascals, | 154 : | cette valeur est sans dimension et vaut 1. | 155 : | Pour une pression qui s'exprimerait en metres, | 156 : | cette valeur serait en Pa/m, et vaudrait rho_w g. | 157 : | | 158 : |'AXE_G' VECTEUR donnant la direction des forces de gravite | 159 : | (facultatif, defaut = 0. 1. en 2D, 0. 0. 1. en 3D) | 160 : | Utile seulement avec l'indice 'FORCE_GRAVITE'. | 161 : | | 162 : |----------------------- | 163 : | Conditions initiales : | 164 : |----------------------- | 165 : | | 166 : |'TEMPS' TABLE contenant a l'indice 0 la valeur du temps | 167 : | initial (s) ('FLOTTANT', defaut = 0.) | 168 : | | 169 : |'CHARGE' TABLE contenant a l'indice 0 la charge initiale (m) | 170 : | (CHPOINT centre, composante 'H') | 171 : | | 172 : |'FLUX' TABLE contenant a l'indice 0 le flux diffusif | 173 : | initial (m3/s) (CHPOINT face, composante 'FLUX') | 174 : | En absence de cet indice le flux n'est pas | 175 : | sauvegarde | 176 : | | 177 : | | 178 : |--------------------------------------- | 179 : | Conditions aux limites / chargements : | 180 : |--------------------------------------- | 181 : | | 182 : |'TRACE_IMPOSE' Valeurs des traces imposees (charge) | 183 : | | 184 : |'FLUX_IMPOSE' Valeurs des flux imposes integres par face | 185 : | (Type CHARGEMENT de CHPO Face) - | 186 : | | 187 : | | 188 : |'FLUXTOT_IMP' Valeurs des flux totaux imposes integres par face | 189 : | (Type CHARGEMENT de CHPO Face, comp. 'H') | 190 : | | 191 : |'MIXTES' Table : - indice C contient les valeurs des flux | 192 : | mixtes imposes integres par face | 193 : | (Type CHARGEMENT de CHPO Face, | 194 : | comp. 'H' defaut 0.) | 195 : | - indices A et B sont des reels | 196 : | | 197 : | la condition mixte s'ecrit | 198 : | C = A * flux diffusif + B * Concentration | 199 : | | 200 : |'SOURCE' source d'eau imposee par maille et unite de temps | 201 : | (option, defaut=0) (type 'CHARGEMENT' 'SOUR') | 202 : | | 203 : |------------------------------------------- | 204 : | Homogeneisation des capacite capillaire et | 205 : | permeabilite relative sur un element : | 206 : |------------------------------------------- | 207 : | | 208 : |'HOMOGENEISATION' mot 'CENTRE' (defaut) ou 'DECENTRE' | 209 : | | 210 : | 'TRACE' : C (capacite capil.) et K (permeabilite) | 211 : | (uniquement avec les EFMH) calculees a partir des | 212 : | traces donnees par les EFMH | 213 : | | 214 : | 'DECENTRE' : C (capacite capil.) et K (permeabilite)| 215 : | calculees par moyenne arithmetique des valeurs | 216 : | aux faces evaluees par interpolation. | 217 : | Permet une meilleure precision lors | 218 : | de la precense de chocs | 219 : | | 220 : | 'CENTRE' : C et K calcules avec la pression moyenne | 221 : | elementaire. moins precis en cas de choc mais | 222 : | le point fixe converge souvent plus vite | 223 : | | 224 : |------------------------------------ | 225 : | Algorithme de resolution (PICARD) : | 226 : |------------------------------------ | 227 : | | 228 : |'ITMAX' nombre maximum d'iterations (defaut = 40) | 229 : | | 230 : |'NITER' Nombre d'iterations recherche pour atteindre la | 231 : | convergence d'un temps calcule : une convergence | 232 : | trop lente induit une diminution automatique du | 233 : | facteur sur le pas de temps (indice 'CFL') et | 234 : | inversement (ENTIER, defaut = 10) | 235 : | | 236 : |'RESIDU_MAX' critere adimensionnel de convergence | 237 : | = quantite volumique relative maximale acceptable | 238 : | de fluide ajoute au systeme pendant un pas de temps.| 239 : | Le residu est construit sur le systeme non relaxe | 240 : | ('FLOTTANT', defaut = 1.D-4) | 241 : | | 242 : |'SOUS_RELAXATION' sous-relaxation sur la capacite et la | 243 : | permeabilite. La valeur de ce parametre doit etre | 244 : | comprise entre 0 et 1. Valeur inferieure a 1 | 245 : | conseillee en presence d'oscillations. | 246 : | La sous-relaxation passe automatiquement a 0,5 au | 247 : | dela de 'ITMAX'/2 iterations (difficulte de | 248 : | convergence averee). | 249 : | Une fois la convergence atteinte, la sous-relaxation| 250 : | retrouve sa valeur initiale pour le temps calcule | 251 : | suivant. | 252 : | (FLOTTANT, defaut = 1 : pas de relaxation) | 253 : | (rappel : la sous-relaxation effectue une moyenne | 254 : | des valeurs d'une variable obtenues sur deux | 255 : | iterations successives) | 256 : | | 257 : |'COFDIV' en l'absence de convergence le pas de temps est | 258 : | diminue de ce facteur ('FLOTTANT', defaut=0.8), et | 259 : | le calcul est relance | 260 : | | 261 : |'NMAXDT' nombre maximum de fois que ce processus de division | 262 : | de pas de temps est renouvelle (ENTIER, defaut = 6) | 263 : | | 264 : |-------------------------- | 265 : | Discretisation en temps : | 266 : |-------------------------- | 267 : | | 268 : |'THETA' Coefficient d'implicitation compris entre 0. et 1. | 269 : | (theta-methode diffusion) (FLOTTANT - defaut 1.) | 270 : | Possibilite de non-convergence lorsque theta<1/2 | 271 : | Valeurs de theta generalement utilisees : | 272 : | Schema de Euler explicite : 0. | 273 : | Schema de Crank-Nicholson : 1/2 | 274 : | Schema de Galerkin : 2/3 | 275 : | Schema de Euler implicite : 1. | 276 : | | 277 : |'TEMPS_CALCULES' LISTREEL des temps a calculer (optionnel) | 278 : | | 279 : | !!! En absence de l'indice 'TEMPS_CALCULES', le pas de !!! | 280 : | !!! temps est calcule automatiquement. Dans ce cas les !!! | 281 : | !!! indices suivants doivent etre renseignes : !!! | 282 : | | 283 : |'NPAS' nombre de pas de temps | 284 : | (ENTIER, defaut = jusqu'a TEMPS_FINAL) | 285 : | | 286 : |'TEMPS_FINAL' temps final (FLOTTANT, obligatoire si NPAS absent) | 287 : | | 288 : |'DT_INITIAL' pas de temps initial ('FLOTTANT', optionel) | 289 : | | 290 : |'CFL' Facteur initial sur le pas de temps automatique. | 291 : | Ce facteur est modifie ensuite au cours de la | 292 : | resolution pour respecter le critere demande sur la | 293 : | rapidite de convergence NITER | 294 : | (FLOTTANT, defaut = 0.5) | 295 : | | 296 : |-------------------- | 297 : |Donnees numeriques : | 298 : |-------------------- | 299 : | | 300 : | 'LUMP' FAUX SI pas de mass lumping, VRAI sinon. | 301 : | VRAI seulement sur des maillages de rectangles et | 302 : | parallelepipedes rectangles et tenseur de dissusion| 303 : | orthotrope. Permet de rendre les schemas monotone | 304 : | pour la diffusion-instationnaire | 305 : | | 306 : | 'DECENTR' VRAI si diffusion numerique pour Peclet = 2, permet| 307 : | de stabiliser (en explicite) voire rendre monotone | 308 : | le schema de convection. | 309 : | FAUX si schema sans convection, ou en implicite et | 310 : | absence d'oscillations - plus precis | 311 : | | 312 : | 'TYPDISCRETISATION' 'VF' si VF et 'EFMH' si EFMH | 313 : | | 314 : | 'CONSERVATIF' Schema conservatif (l'option par defaut ne | 315 : | l'est pas). Cette option peut donner des resultats plus precis | 316 : | dans certains cas. | 317 : | | 318 : | 'METHINV' TABLE DE PARAMETRE du solveur KRES, cf KRES | 319 : | peut etre remplie partiellement | 320 : | deux indices importants : | 321 : | 'TYPINV' OBLIGATOIRE 1 pour direct 3 pour BiCGSTAB | 322 : | 'PRECOND' obligatoire 3 pour ILU0, 5 pour ILUT | 323 : | conseil TYPINV = 1 en 2D ou petits calculs, 3 sinon| 324 : | conseil PRECOND = 3 sauf si problemes mettre 5 | 325 : | DERNIER CONSEIL : si message du type convergence | 326 : | breakdown, Pivot nul ... mettre la tolerence | 327 : | BCGSTOL a la precision machine 1.D-300, cf notice | 328 : | de KRES | 329 : | | 330 : |--------------------------- | 331 : | Sauvegarde des resultats : | 332 : |--------------------------- | 333 : | | 334 : |'TEMPS_SAUVES' liste des indices des temps a sauver. | 335 : | ('LISTREEL', defaut = tous) | 336 : | | 337 : |'MESSAGES' (type ENTIER, 0 a 2) niveau de detail des messages | 338 : | imprimes a l'execution (0=peu, 1=normal, 2=tous) | 339 : | | 340 : ------------------------------------------------------------------- 341 : | Sortie : | 342 : ---------- | 343 : | | 344 : |'TEMPS' TABLE contenant les temps sauvegardes ('FLOTTANT') | 345 : | | 346 : |'CHARGE' TABLE contenant a la charge hydrauliquea (m ou Pa) | 347 : | (CHPOINT centre, comp. 'H', nature diffus) | 348 : | | 349 : |'FLUX' TABLE contenant les flux hydrauliques (m3/s) | 350 : | (CHPOINT face, comp. 'FLUX', nature discret) | 351 : | (optionnel, seulement si l'indice 0 a ete defini | 352 : | | 353 : |'SUCCION' TABLE contenant la succion Pw - Pg (m ou Pascal) | 354 : | (CHPOINT centre, comp. 'SCAL') | 355 : | | 356 : |'PRESSION' TABLE contenant la pression (m ou Pascal) | 357 : | (CHPOINT centre, comp. 'SCAL') | 358 : | | 359 : |'SATURATION' TABLE contenant la saturation (s.d. entre 0 et 1) | 360 : | (CHPOINT centre, comp. 'SCAL') | 361 : | | 362 : |'TH2O' TABLE contenant la teneur en eau | 363 : | (s.d. entre 0 et la porosite) | 364 : | (CHPOINT centre, comp. 'SCAL') | 365 : | | 366 : |'TEMPS_EFFECTUES' LISTREEL contenant les temps effectivement | 367 : | calcules. Utile avec un pas de temps automatique et | 368 : | les TEMPS_SAUVES specifies | 369 : | | 370 : |_________________________________________________________________| 371 : 372 : 373 : En sortie, TAB1 permet de retrouver les resultats. Ceux-ci sont mis 374 : dans des tables dont les indices sont des entiers ( 0 1 2 ... N ) 375 : correspondants aux numeros de sauvegarde des resultats, l'indice 0 376 : designant le temps initial. 377 : 378 : Exemple : pour lister le CHPOINT de charge calcule pour la 379 : valeur du parametre d'evolution 2.5, il faudra coder : 380 : LIST ( PECHE TAB1 CHARGE 2.5 ) ; 381 : ou 382 : 383 : si on sait que l'indice i de la table TEMPS contient 384 : la valeur du parametre d'evolution 2.5, on peut coder 385 : LIST ( TAB1 . CHARGE . i ) ; 386 : 387 : des profils ou des evolutions peuvent etre obtenus 388 : avec les procedures TRACHIS et TRACHIT 389 : 390 : Remarques : 391 : ___________ 392 : 393 : 1 - Les resultats etant reperes dans TAB1. Il n'y a pas d'objets 394 : nommes crees par cette procedure. 395 : 396 : 2 - En presence d'un pas de temps automatique, un nombre d'iterations 397 : demande (indice 'NITER') trop important induit un CFL grand 398 : (> 0.5). Le CFL est fourni en commentaire a chaque pas de temps 399 : sii le niveau de message est superieur ou egal a 1 400 : 401 : 3 - Lorsque le milieu devient sature dans sa totalite, le pas de 402 : temps automatique devient tres grand. L'equation resolue 403 : redevenant lineaire, la convergence du processus iteratif est 404 : immediate. 405 : 406 : 4 - Des messages permettent de suivre l'evolution d'un certain 407 : nombres de parametre au cours du processus iteratif 408 : (si 'MESSAGE' > 0) : 409 : - valeur du critere de convergence 410 : - maximum et minimum de la trace de charge 411 : - minimum de la permeabilite 412 : - maximum de la capacite capillaire 413 : - nombre de mailles non saturees 414 : 415 : 5 - Une fois sorti de DARCYSAT on peut y re-entrer en definissant 416 : soit le temps final soit le nombre de pas de temps et en invoquant 417 : de nouveau DARCYSAT avec les memes operandes que lors du premier 418 : appel. 419 : 420 : 6 - Les calculs peuvent etre effectues avec un autre systeme d'unite 421 : que celui propose dans la notice (en particulier les charges en 422 : Pascal). L'utilisateur doit alors veiller lui-meme a la coherence 423 : des unites des differents parametres qu'il utilise et renseigner 424 : l'indice 'CONVERSION_CHARGE'. 425 : 426 : 7 - La loi de saturation en fonction de la pression capillaire se 427 : trouve dans la procedure HT_PRO, celle de permeabilite dans la 428 : procedure KR_PRO. L'utilisateur peut les modifier ou construire 429 : ses propres procedures sur le meme moule, a condition de preciser 430 : leur(s) nom(s) dans les tables stoquees aux indices 431 : 'LOI_SATURATION' et 'LOI_PERMEABILITE' respectivement. 432 :
© Cast3M 2003 - Tous droits réservés.
Mentions légales