$$$$ 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' Description : _____________ This procedure calculates the transient flow throught a saturated / non saturated porous media. Darcy's equation apply in the saturated zones, and Richards' in the unsaturated zone (writen here with head in meters). saturated zone (Pw - Pg > 0) Se dh/dt = -div(U) ; U = -Ks grad(h) h = P/rho*g + z or h = (Pw-Pg) + rho*g*z (Pa) where, Se, storativity coefficient (1/m), Ks, saturation permeability (m/s), Pw, water pressure (Pa), Pg, gaz pressure (Pa), h, hydraulic head (m), U, Darcy's velocity (m/s). unsaturated zone (Pw - Pg < 0) C(P) dh/dt = -div(U) ; U = -K(P) (grad h) ou, C(P) capacite capillaire (m-1) et K(P), permeabilite (m/s) The resolutions uses mixed finite element discretisation. Water content is deduced from pression through procedure HT_PRO. Permeability is deduced from saturation by procedure KR_PRO. Capillary capacity is deduced from Water content and pressure. pression. Automatic mecanisms are available to increase the convergence of this non-linear problem : capacity and permeability relaxation, element averaging of physical properties, automatique time-step refinement, penalisation. Contents : __________ Modeling options and parameters must be enterred in TAB1 with the following indexes. ___________________________________________________________________ | | | Indexes definition | | | ------------------------------------------------------------------- | | |---------------------------------------- | | Physical, geometrical et material data: | |---------------------------------------- | | | |'SOUSTYPE' 'DARCY_TRANSATUR' (type MOT) | | | |'MODELE' model object (type MMODEL, see MODE) | | | |'LOI_PERMEABILITE' table structure (type TABLE) containing | | permeability parameters used by KR_PRO procedure. | | 3 subtype indexes can be used : | | | | - subtype : | | 'PUISSANCE', 'MUALEM', 'BURDINE', 'MUALEM_BURDINE' | | 'BROOKS_COREY', 'EXPONENTIELLE', 'LOGARITHMIQUE' | | : indexes are those used in | | KR_PRO standard procedure | | | | - subtype 'PERSONNELLE' : indexes are those used in | | KR_PRO user procedure, or in another whose name is | | specified at the 'NOM_PROCEDURE' index | | | | - subtype 'MULTIZONE' : each index receives a table | | structure defining permeability parameters in parts | | of the modeled field (subtypes | | 'PUISSANCE', 'MUALEM', 'BURDINE', 'MUALEM_BURDINE' | | 'BROOKS_COREY', 'EXPONENTIELLE', 'LOGARITHMIQUE' | | or | | 'PERSONNELLE'). Indexes are free. These tables must | | contain the model structure at the index 'MODELE', | | these sub-models building a partition of the overall| | domain. | | | |'LOI_SATURATION' table structure (type TABLE) containing | | saturation parameters used by HT_PRO procedure. | | 3 subtype indexes can be used : | | | | - subtype 'VAN_GENUCHTEN', | | 'EXPONENTIELLE', 'LOGARITHMIQUE' | | : indexes are those used in | | HT_PRO standard procedure | | | | - subtype 'PERSONNELLE' : indexes are those used in | | HT_PRO user procedure, or in another whose name is | | specified at the 'NOM_PROCEDURE' index | | | | - subtype 'MULTIZONE' : each index receives a table | | structure defining saturation parameters in parts | | of the modeled field (subtypes 'VAN_GENUCHTEN', | | 'EXPONENTIELLE', 'LOGARITHMIQUE' | | or | | 'PERSONNELLE'). Indexes are free. These tables must | | contain the model structure at the index 'MODELE', | | these sub-models building a partition of the overall| | domain. | | | | sub-indexes contained for each sub-law : | | | | 'PREC' value of the minimal water pressure | | separation for capillary pressure calculation | | (optional, float, same unit as head) | | | | for the other (defaut) sub-indexes, see HT_PRO | | | |'COEFEMMA' Specific storativity when saturated (m-1) | | (FLOTTANT or 'CHPOINT', comp. 'SCAL', nature DISCRET| | default = 0) | | | |'FORCE_GRAVITE' (type 'FLOTTANT', default=1000*9.81) gravity | | value (N/m3), equal to water specific mass | | multiplied by gravity acceleration. Without this | | index, gravity forces are ignored. | | | |'CONVERSION_CHARGE' (type 'FLOTTANT', default = 1.) conversion | | factor for head to Pascals | | Pressure head in Pascals -> 1. | | Pressure head in meters -> rho_w * g | | | |'AXE_G' VECTEUR giving the direction of gravity forces | | (optional, defaut = 0. 1. in 2D, 0. 0. 1. in 3D) | | Useful only with 'FORCE_GRAVITE' index. | | | |--------------------- | | Initial conditions : | |--------------------- | | | |'TEMPS' TABLE containing the initial time at the 0 index | | (s) ('FLOTTANT', default = 0.) | | | |'TRACE_CHARGE' TABLE containing initial head's trace at the | | 0 index ('CHPOINT', Face, Comp 'TH') | | | |'CHARGE' TABLE containing initial head at the 0 index | | ('CHPOINT', Center, Comp 'H') | | | |'FLUX' TABLE containing initial integrated diffusive flux | | at the 0 index (Type 'CHPO' Face, Comp 'FLUX') | | without this index, the diffusive flux is not saved| | | |------------------------------------ | | Boundary Conditions / Chargements : | |------------------------------------ | | | |'TRACE_IMPOSE' Values of imposed heads on boundaries | | | |'FLUX_IMPOSE' Values of imposed fluxes integrated on each face | | (Type CHARGEMENT de CHPO Face) - | | | |'FLUXTOT_IMP' Values of the imposed total flux integrated on face| | (Type CHARGEMENT de CHPO Face, comp. 'H') | | | |'MIXTES' Table : - index C contains imposed flux values | | integrated on each face | | (Type CHARGEMENT of CHPO Face, | | comp. 'H' default 0.) | | - indices A and B are float | | | | TH mixed conditions writes : | | C = A * flux + B * Head | | | |'SOURCE' imposed water source integrated on each cell for a | | time unit | | (option, default=0) (type 'CHARGEMENT' 'SOUR') | | | |------------------------------------------------------------------| | Capilarity capacity and permeability homogeneisation method | |------------------------------------------------------------------| | | |'HOMOGENEISATION' two possible values | | (type 'MOT', optional, defaut 'CENTRE') : | | | | 'TRACE' : C (capil. capacity) and K (permeability) | | (only for MHFE) deduced from traces given by MHFE | | | | 'DECENTRE': C (capil. capacity) and K (permeability)| | deduced from pressure head trace. Stable, | | average precision | | | | 'CENTRE' : C et K deduced from mean pressure head in| | cells. may be unstable, godd precision | | | |-------------------------------- | | Resolution algorithm (PICARD) : | |-------------------------------- | | | |'ITMAX' maximum number of iteration (default=40) | | | |'NITER' objective number of iterations before reaching | | convergence. A slow (resp. fast) convergence rate | | induces a decrease (resp. increase) in 'CFL' factor | | (FLOTTANT, default = 0.5) | | | |'RESIDU_MAX' convergence criterion (default = 1.D-4) | | = maximal acceptable relative volumic error on mass | | conservation | | | |'SOUS_RELAXATION' sub-relaxation coefficient for capacity and | | permeability, between 0. and 1. (default = 1.) | | sub-relaxation coefficient value lowers to 0.5 after| | ITMAX/2 iterations. | | After convergence, sub-relaxation coefficient is | | refreshed to its initial value | | | |'XI' artificial storativity used to improve convergence | | The new equation is | | (Xi + C) dh/dt = -div(U) + Xi dh/dt | | First member Xi term is implicite and second member | | terms is explicite. | | If this index is missing, the code will | | modify the time step if convergence is not reached. | | ('FLOTTANT', default = 0.) | | | |'DIVISION_DT' Logicat indicating time-step division when non | | convergency, and use of the following two indexes | | | |'COFDIV' time step is divided| | by this amount (type 'FLOTTANT', default = 0.8) and | | the calculation is started over | | | |'NMAXDT' maximal number of this time step division | | (type 'ENTIER', default = 6) | | | |-------------------- | |numerical parmeters : | |-------------------- | | | | 'LUMP' FAUX if no mass lumping, VRAI otherwise. | | VRAI on rectangular mesh ONLY!! | | And permeability tensor must be orthotrope | | enables monotonicity of the scheme | | | | 'TYPDISCRETISATION' 'VF' if VF and 'EFMH' if EFMH | | | | 'CONSERVATIF' conservative scheme | | | | 'METHINV' TABLE OF PARAMETERS for solver KRES, cf KRES | | might be partially filled | | two important indices : | | 'TYPINV' COMPULSARY, 1 for direct 3 for BiCGSTAB | | 'PRECOND' COMPULSARY, 3 for ILU0, 5 for ILUT | | | |------------------------ | | Time discretisation : | |------------------------ | | | |'THETA' Relaxation coefficient for the diffusive fluxes. | | Value between 0. and 1. (theta-method) | | (FLOTTANT - default 1.) | | Unstability can appear if theta<1/2. Here is the | | list of current theta's scheme : | | Euler's explicit scheme : 0. | | Crank-Nicholson's scheme : 1/2 | | Galerkin's scheme : 2/3 | | Euler's implicit scheme : 1. | | | | | |'TEMPS_CALCULES' values of the calculated times (optional, type | | LISTREEL). | | | | !!! If this index is missing, the code will adopt an !!! | | !!! adaptative time step. Then user must gives !!! | | !!! the following data : !!! | | | |'NPAS' number of time steps, ignored if 'TEMPS_CALCULES' | | (ENTIER, defaut = until TEMPS_FINAL) | | | |'TEMPS_FINAL' final time ('FLOTTANT', necessary if not 'NPAS') | | | |'DT_INITIAL' initial time step ('FLOTTANT', optional) | | | |'CFL' initial factor on the adaptative time step. | | It is then modified during the resolution in order | | to respect the objective number of iterations NITER | | (FLOTTANT, defaut = 0.5) | | | |------------------ | | Desired results : | |------------------ | | | | 'TEMPS_SAUVES' list of the indexes to be saved. | | ('LISTREEL', default = tous) | | | | 'MESSAGES' (type ENTIER, 0 to 2) level of detail in displayed | | messages (0=few, 1=normal, 2=all) | | | ------------------------------------------------------------------- | Output : | |--------- | | | |'TEMPS' TABLE containing les temps sauvegardes ('FLOTTANT') | | | |'TRACE_CHARGE' TABLE containing hydraulic head traces (m or Pa) | | (CHPOINT face, comp. 'TH', nature diffus) | | | |'CHARGE' TABLE containing hydraulic head (m or Pa) | | (CHPOINT centre, comp. 'H', nature diffus) | | | |'FLUX' TABLE containing hydraulic flux (m3/s) | | (CHPOINT face, comp. 'FLUX', nature discret) | | (optional, only if its index 0 is defined | | | |'SUCCION' TABLE containing succion Pw - Pg (m or Pascal) | | (CHPOINT centre, comp. 'SCAL') | | | |'PRESSION' TABLE containing hydraulic pressure (m or Pascal) | | (CHPOINT centre, comp. 'SCAL') | | | |'SATURATION' TABLE containing saturation (w.d. between 0 and 1) | | (CHPOINT centre, comp. 'SCAL') | | | |'TH2O' TABLE containing water content | | (s.d. between 0 and porosity) | | (CHPOINT centre, comp. 'SCAL') | | | |'TEMPS_EFFECTUES' LISTREEL containing the effective calculated | | times. Usefull with an automatique time step and | | specified TEMPS_SAUVES | | | |---------------------------------------------- | | for PASAPAS if hydro-mechanical calculation : | |---------------------------------------------- | | (if 'HYDRO_MECA' index present and VRAI) | | | |'SATU' saturation (centre ou face) | |'PRES' hydraulic pressure (centre, Pa) | |'CHRG' hydraulic head (centre, meters or Pa) | |'TRH' hydraulic head trace (face, meters ou Pa) | | | |_________________________________________________________________| ___________________________________________________________________ | | | Results tables are indexed by integers varying from 0 to N | | At index 0 are stored the initial values, and following index | | yield the corresponding values at stored times TEMPS.I . | | The fields used to resume calculation are those stored at last | | index. | |_________________________________________________________________| Example : in order to list the saturation CHPOINT calculated for the value of the evolution parameter 2.5, the following will have to be coded : LIST ( PECHE TAB1 'SATURATION' 2.5 ) ; or if we know that the index i of the saved time table TEMPS contains 2.5, we can code : LIST ( TAB1 . 'CHARGE' . i ) ; Remarks : _________ 1 - Since the results are located in TAB1, no named objects are created by this procedure. 2 - When time step is automatic, the CFL coefficient is decreased (resp. increased) according to the excess (resp. deficit) of the number of iterations required to reach convergence in comparison with the objective NITER. 3 - When the domain is fully saturated, the time step becomes very high. The solved equation becoming linear, the convergency of the iterative algorithm is immediate. 4 - The value of certain parameters are displayed during calculation (if 'MESSAGE' > 0) : - residual error on mass conservation - hydraulic head trace extrema - permeability minimum - capillary capacity minimum - number of desaturated cells 5 - Once you have left DARCYSAT you may re-enter it by defining new values for TEMPS_CALCULES or TEMPS_FINAL and by calling once again DARCYSAT with the same operands as for the first call. 6 - Computations can be carried out with other unit systems (particularly, using something else than Pascal as pressure unit). The user must then insure himself the consistency in the parameter units and fill in the 'CONVERSION_CHARGE' index. 7 - The saturation law, as a function of capillary pressure is defined by the HT_PRO procedur, the one relating water permeability to saturation is defined by the KR_PRO procedur. The user can either modify them, or build another one (other ones), provided their name(s) is stored in the sub-tables 'LOI_SATURATION' and 'LOI_PERMEABILITE'. $$$$
© Cast3M 2003 - All rights reserved.
Disclaimer