Afficher cette notice en
$$$$ DARCYSAT NOTICE CHAT 11/09/12 21:15:42 7124 DATE 11/09/12 Procedure DARCYSAT ------------------ Voir aussi : HT_PRO KR_PRO PECHE TRACHIS TRACHIT TAB1.'SOUSTYPE'.'MODELE' .'LOI_PERMEABILITE'.'LOI_SATURATION'.'HOMOGENEISATION' .'COEFEMMA'.'BLOCAGES_DARCY'.'CHARGEMENT'.'FORCE_GRAVITE' .'CONVERSION_CHARGE'.'TEMPS'.'TRACE_CHARGE'.'CHARGE'.('FLUX') .'TEMPS_CALCULES'.'TEMPS_SAUVES'.'TEMPS_FINAL' .'DT_INITIAL'.'NPAS'.'CFL' .'NITER'.'ITMAX'.'RESIDU_MAX'.'SOUS_RELAXATION'.'XI' .'DIVISION_DT'.'COFDIV'.'NMAXDT'.'MESSAGE' Objet : _______ Cette procedure permet de simuler un transitoire d'ecoulement en zones saturee et non saturee d'un milieu poreux. L'ecoulement est decrit en zone saturee par l'equation de Darcy et en zone non saturee par l'equation de Richard's (ecrite ici avec h en metres) zone saturee (Pw - Pg > 0) Se dh/dt = -div(U) ; U = -Ks grad(h) ; h = (Pw-Pg)/rho*g + z (m) ou h = (Pw-Pg) + rho*g*z (Pa) avec, Se, coefficient d'emmagasinement (1/m), Ks, permeabilite a saturation (m/s), Pw, pression d'eau (Pa), Pg, pression de gaz (Pa), h, charge (m) et U, vitesse de Darcy (m/s). zone non saturee (Pw - Pg < 0) C(P) dh/dt = -div(U) ; U = -K(P) (grad h) avec, C(P) capacite capillaire (1/m) K(P), permeabilite (m/s) La resolution est effectuee en trace de charge par la methode EFMH. La teneur en eau se deduit de la pression par la procedure HT_PRO. La permeabilite se deduit de la saturation par la procedure KR_PRO. La capacite capillaire est deduite de la teneur en eau et de la pression. Le probleme etant non lineaire, des mecanismes automatiques ou non ont ete mis en place pour fiabiliser la convergence du systeme : relaxation de la capacite et de la conductivite hydraulique et homogeneisation des caracteristiques decentree sur les elements, changement automatique de pas de temps ou penalisation sur le terme transitoire. Commentaire : _____________ En entree, TAB1 sert a definir les options et les parametres du calcul. Les indices de la table TAB1 sont des mots (a coder tels quels) dont voici la description : ___________________________________________________________________ | | | Indice Definition | | | ------------------------------------------------------------------- | | |------------------------------------------------- | | Donnees physiques, geometriques et materielles : | |------------------------------------------------- | | | |'SOUSTYPE' 'DARCY_TRANSATUR' (type MOT) | | | |'MODELE' Objet modele (MMODEL cree par MODE) | | | |'LOI_PERMEABILITE' Objet table contenant les valeurs des | | parametres de la loi de permeabilite calculee par la| | procedure KR_PRO. | | 3 indices 'SOUSTYPE' possibles : | | | | - soustypes au choix | | 'PUISSANCE', 'MUALEM', 'BURDINE', 'MUALEM_BURDINE' | | 'BROOKS_COREY', 'EXPONENTIELLE', 'LOGARITHMIQUE' | | : la table doit contenir les | | indices designant les parametres utilises dans la | | procedure standard KR_PRO. | | | | - soustype 'PERSONNELLE' : la table doit contenir | | les indices de la loi construite par l'utilisateur | | dans sa propre procedure KR_PRO, ou dans une autre | | dont le nom est stipule au sous-indice | | 'NOM_PROCEDURE'. | | | | - soustype 'MULTIZONE' : les indices doivent | | designer les differentes zones, les valeurs de ces | | indices etant alors des tables de soustype | | 'PUISSANCE', 'MUALEM', 'BURDINE', 'MUALEM_BURDINE' | | 'BROOKS_COREY', 'EXPONENTIELLE', 'LOGARITHMIQUE' | | documentees dans la notice de KR_PRO | | ou'PERSONNELLE' dont les indices | | designent les parametres des lois affectees a chaque| | zone. Ces sous-tables doivent aussi contenir le | | modele de la zone a l'indice 'MODELE', ces sous- | | modeles constituant une partition du domaine general| | | | pour les autres sous-indices (defaut), voir KR_PRO | | | |'LOI_SATURATION' Objet table contenant les valeurs des | | parametres de la loi analytique de teneur en eau | | calculee par la procedure HT_PRO. | | 3 indices 'SOUSTYPE' possibles : | | | | - soustypes aux choix 'VAN_GENUCHTEN', | | 'EXPONENTIELLE', 'LOGARITHMIQUE' | | : la table doit contenir les | | indices designant les parametres utilises dans la | | procedure standard HT_PRO. | | | | - soustype 'PERSONNELLE' : la table doit contenir | | les indices de la loi construite par l'utilisateur | | dans sa propre procedure HT_PRO, ou dans une autre | | dont le nom est stipule au sous-indice | | 'NOM_PROCEDURE'. | | | | - soustype 'MULTIZONE' : les indices doivent | | designer les differentes zones, les valeurs de ces | | indices etant alors des tables de soustype | | 'VAN_GENUCHTEN', 'EXPONENTIELLE', 'LOGARITHMIQUE' | | ou 'PERSONNELLE' dont les indices designent les | | parametres des lois affectees a chaque zone. | | Ces sous-tables doivent aussi contenir le modele de | | la zone a l'indice 'MODELE', ces sous-modeles | | constituant une partition du domaine general | | | | pour les autres sous-indices (defaut), voir HT_PRO | | | |'COEFEMMA' coefficient d'emmagasinement en zone saturee (m-1) | | (FLOTTANT ou CHPOINT, comp. 'SCAL', nature | | discrete, defaut = 0) | | | |'FORCE_GRAVITE' (type FLOTTANT, par ex. = 1000*9.81) valeur de la| | force de gravite (N/m3) egale a la masse volumique | | de l'eau multipliee par l'acceleration de la | | pesanteur. En absence de ce mot cle, la gravite | | n'est pas prise en compte. | | | |'CONVERSION_CHARGE' (type 'FLOTTANT', defaut = 1.) facteur de | | convertion de la charge en Pascals. | | Pour une pression qui s'exprime en Pascals, | | cette valeur est sans dimension et vaut 1. | | Pour une pression qui s'exprimerait en metres, | | cette valeur serait en Pa/m, et vaudrait rho_w g. | | | |'AXE_G' VECTEUR donnant la direction des forces de gravite | | (facultatif, defaut = 0. 1. en 2D, 0. 0. 1. en 3D) | | Utile seulement avec l'indice 'FORCE_GRAVITE'. | | | |----------------------- | | Conditions initiales : | |----------------------- | | | |'TEMPS' TABLE contenant a l'indice 0 la valeur du temps | | initial (s) ('FLOTTANT', defaut = 0.) | | | |'CHARGE' TABLE contenant a l'indice 0 la charge initiale (m) | | (CHPOINT centre, composante 'H') | | | |'FLUX' TABLE contenant a l'indice 0 le flux diffusif | | initial (m3/s) (CHPOINT face, composante 'FLUX') | | En absence de cet indice le flux n'est pas | | sauvegarde | | | | | |--------------------------------------- | | Conditions aux limites / chargements : | |--------------------------------------- | | | |'TRACE_IMPOSE' Valeurs des traces imposees (charge) | | | |'FLUX_IMPOSE' Valeurs des flux imposes integres par face | | (Type CHARGEMENT de CHPO Face) - | | | | | |'FLUXTOT_IMP' Valeurs des flux totaux imposes integres par face | | (Type CHARGEMENT de CHPO Face, comp. 'H') | | | |'MIXTES' Table : - indice C contient les valeurs des flux | | mixtes imposes integres par face | | (Type CHARGEMENT de CHPO Face, | | comp. 'H' defaut 0.) | | - indices A et B sont des reels | | | | la condition mixte s'ecrit | | C = A * flux diffusif + B * Concentration | | | |'SOURCE' source d'eau imposee par maille et unite de temps | | (option, defaut=0) (type 'CHARGEMENT' 'SOUR') | | | |------------------------------------------- | | Homogeneisation des capacite capillaire et | | permeabilite relative sur un element : | |------------------------------------------- | | | |'HOMOGENEISATION' mot 'CENTRE' (defaut) ou 'DECENTRE' | | | | 'TRACE' : C (capacite capil.) et K (permeabilite) | | (uniquement avec les EFMH) calculees a partir des | | traces donnees par les EFMH | | | | 'DECENTRE' : C (capacite capil.) et K (permeabilite)| | calculees par moyenne arithmetique des valeurs | | aux faces evaluees par interpolation. | | Permet une meilleure precision lors | | de la precense de chocs | | | | 'CENTRE' : C et K calcules avec la pression moyenne | | elementaire. moins precis en cas de choc mais | | le point fixe converge souvent plus vite | | | |------------------------------------ | | Algorithme de resolution (PICARD) : | |------------------------------------ | | | |'ITMAX' nombre maximum d'iterations (defaut = 40) | | | |'NITER' Nombre d'iterations recherche pour atteindre la | | convergence d'un temps calcule : une convergence | | trop lente induit une diminution automatique du | | facteur sur le pas de temps (indice 'CFL') et | | inversement (ENTIER, defaut = 10) | | | |'RESIDU_MAX' critere adimensionnel de convergence | | = quantite volumique relative maximale acceptable | | de fluide ajoute au systeme pendant un pas de temps.| | Le residu est construit sur le systeme non relaxe | | ('FLOTTANT', defaut = 1.D-4) | | | |'SOUS_RELAXATION' sous-relaxation sur la capacite et la | | permeabilite. La valeur de ce parametre doit etre | | comprise entre 0 et 1. Valeur inferieure a 1 | | conseillee en presence d'oscillations. | | La sous-relaxation passe automatiquement a 0,5 au | | dela de 'ITMAX'/2 iterations (difficulte de | | convergence averee). | | Une fois la convergence atteinte, la sous-relaxation| | retrouve sa valeur initiale pour le temps calcule | | suivant. | | (FLOTTANT, defaut = 1 : pas de relaxation) | | (rappel : la sous-relaxation effectue une moyenne | | des valeurs d'une variable obtenues sur deux | | iterations successives) | | | |'COFDIV' en l'absence de convergence le pas de temps est | | diminue de ce facteur ('FLOTTANT', defaut=0.8), et | | le calcul est relance | | | |'NMAXDT' nombre maximum de fois que ce processus de division | | de pas de temps est renouvelle (ENTIER, defaut = 6) | | | |-------------------------- | | Discretisation en temps : | |-------------------------- | | | |'THETA' Coefficient d'implicitation compris entre 0. et 1. | | (theta-methode diffusion) (FLOTTANT - defaut 1.) | | Possibilite de non-convergence lorsque theta<1/2 | | Valeurs de theta generalement utilisees : | | Schema de Euler explicite : 0. | | Schema de Crank-Nicholson : 1/2 | | Schema de Galerkin : 2/3 | | Schema de Euler implicite : 1. | | | |'TEMPS_CALCULES' LISTREEL des temps a calculer (optionnel) | | | | !!! En absence de l'indice 'TEMPS_CALCULES', le pas de !!! | | !!! temps est calcule automatiquement. Dans ce cas les !!! | | !!! indices suivants doivent etre renseignes : !!! | | | |'NPAS' nombre de pas de temps | | (ENTIER, defaut = jusqu'a TEMPS_FINAL) | | | |'TEMPS_FINAL' temps final (FLOTTANT, obligatoire si NPAS absent) | | | |'DT_INITIAL' pas de temps initial ('FLOTTANT', optionel) | | | |'CFL' Facteur initial sur le pas de temps automatique. | | Ce facteur est modifie ensuite au cours de la | | resolution pour respecter le critere demande sur la | | rapidite de convergence NITER | | (FLOTTANT, defaut = 0.5) | | | |-------------------- | |Donnees numeriques : | |-------------------- | | | | 'LUMP' FAUX SI pas de mass lumping, VRAI sinon. | | VRAI seulement sur des maillages de rectangles et | | parallelepipedes rectangles et tenseur de dissusion| | orthotrope. Permet de rendre les schemas monotone | | pour la diffusion-instationnaire | | | | 'DECENTR' VRAI si diffusion numerique pour Peclet = 2, permet| | de stabiliser (en explicite) voire rendre monotone | | le schema de convection. | | FAUX si schema sans convection, ou en implicite et | | absence d'oscillations - plus precis | | | | 'TYPDISCRETISATION' 'VF' si VF et 'EFMH' si EFMH | | | | 'CONSERVATIF' Schema conservatif (l'option par defaut ne | | l'est pas). Cette option peut donner des resultats plus precis | | dans certains cas. | | | | 'METHINV' TABLE DE PARAMETRE du solveur KRES, cf KRES | | peut etre remplie partiellement | | deux indices importants : | | 'TYPINV' OBLIGATOIRE 1 pour direct 3 pour BiCGSTAB | | 'PRECOND' obligatoire 3 pour ILU0, 5 pour ILUT | | conseil TYPINV = 1 en 2D ou petits calculs, 3 sinon| | conseil PRECOND = 3 sauf si problemes mettre 5 | | DERNIER CONSEIL : si message du type convergence | | breakdown, Pivot nul ... mettre la tolerence | | BCGSTOL a la precision machine 1.D-300, cf notice | | de KRES | | | |--------------------------- | | Sauvegarde des resultats : | |--------------------------- | | | |'TEMPS_SAUVES' liste des indices des temps a sauver. | | ('LISTREEL', defaut = tous) | | | |'MESSAGES' (type ENTIER, 0 a 2) niveau de detail des messages | | imprimes a l'execution (0=peu, 1=normal, 2=tous) | | | ------------------------------------------------------------------- | Sortie : | ---------- | | | |'TEMPS' TABLE contenant les temps sauvegardes ('FLOTTANT') | | | |'CHARGE' TABLE contenant a la charge hydrauliquea (m ou Pa) | | (CHPOINT centre, comp. 'H', nature diffus) | | | |'FLUX' TABLE contenant les flux hydrauliques (m3/s) | | (CHPOINT face, comp. 'FLUX', nature discret) | | (optionnel, seulement si l'indice 0 a ete defini | | | |'SUCCION' TABLE contenant la succion Pw - Pg (m ou Pascal) | | (CHPOINT centre, comp. 'SCAL') | | | |'PRESSION' TABLE contenant la pression (m ou Pascal) | | (CHPOINT centre, comp. 'SCAL') | | | |'SATURATION' TABLE contenant la saturation (s.d. entre 0 et 1) | | (CHPOINT centre, comp. 'SCAL') | | | |'TH2O' TABLE contenant la teneur en eau | | (s.d. entre 0 et la porosite) | | (CHPOINT centre, comp. 'SCAL') | | | |'TEMPS_EFFECTUES' LISTREEL contenant les temps effectivement | | calcules. Utile avec un pas de temps automatique et | | les TEMPS_SAUVES specifies | | | |_________________________________________________________________| En sortie, TAB1 permet de retrouver les resultats. Ceux-ci sont mis dans des tables dont les indices sont des entiers ( 0 1 2 ... N ) correspondants aux numeros de sauvegarde des resultats, l'indice 0 designant le temps initial. Exemple : pour lister le CHPOINT de charge calcule pour la valeur du parametre d'evolution 2.5, il faudra coder : LIST ( PECHE TAB1 CHARGE 2.5 ) ; ou si on sait que l'indice i de la table TEMPS contient la valeur du parametre d'evolution 2.5, on peut coder LIST ( TAB1 . CHARGE . i ) ; des profils ou des evolutions peuvent etre obtenus avec les procedures TRACHIS et TRACHIT Remarques : ___________ 1 - Les resultats etant reperes dans TAB1. Il n'y a pas d'objets nommes crees par cette procedure. 2 - En presence d'un pas de temps automatique, un nombre d'iterations demande (indice 'NITER') trop important induit un CFL grand (> 0.5). Le CFL est fourni en commentaire a chaque pas de temps sii le niveau de message est superieur ou egal a 1 3 - Lorsque le milieu devient sature dans sa totalite, le pas de temps automatique devient tres grand. L'equation resolue redevenant lineaire, la convergence du processus iteratif est immediate. 4 - Des messages permettent de suivre l'evolution d'un certain nombres de parametre au cours du processus iteratif (si 'MESSAGE' > 0) : - valeur du critere de convergence - maximum et minimum de la trace de charge - minimum de la permeabilite - maximum de la capacite capillaire - nombre de mailles non saturees 5 - Une fois sorti de DARCYSAT on peut y re-entrer en definissant soit le temps final soit le nombre de pas de temps et en invoquant de nouveau DARCYSAT avec les memes operandes que lors du premier appel. 6 - Les calculs peuvent etre effectues avec un autre systeme d'unite que celui propose dans la notice (en particulier les charges en Pascal). L'utilisateur doit alors veiller lui-meme a la coherence des unites des differents parametres qu'il utilise et renseigner l'indice 'CONVERSION_CHARGE'. 7 - La loi de saturation en fonction de la pression capillaire se trouve dans la procedure HT_PRO, celle de permeabilite dans la procedure KR_PRO. L'utilisateur peut les modifier ou construire ses propres procedures sur le meme moule, a condition de preciser leur(s) nom(s) dans les tables stoquees aux indices 'LOI_SATURATION' et 'LOI_PERMEABILITE' respectivement.
© Cast3M 2003 - Tous droits réservés.
Mentions légales