$$$$ EXEC NOTICE GOUNAND 11/12/21 21:15:01 7232 DATE 11/12/21 Procedure EXEC Voir aussi : EQEX MODE -------------- DOMA KCHT KRES EXEC TAB1 ; DESCRIPTION : ____________ I/ SCALAR DIFFUSION/CONVECTION This procedure allows to compute the transport (diffusion/convection) of a scalar, in a transient or a steady state regime, with various boundary conditions (prescribed values, fluxes, exchanges, sources). The computation can be linear or not. The non linearities are resolved by a fixed point method. The description of the equation to be solved is done using EQEX operator which creates the table TAB1. The scalar can represent as well as a temperature or a concentration. The algorithm parameters are defined in TAB1 table using EQEX operator. I.1/ Explicit transient computation. The time step is limited by a stability condition. It can be imposed or automatically computed. The non linearities, in particular for the physical properties can be solved by a new computation at each time step in a 'user' procedure. The steady state can be reached as the asymptotic limit of the transient (if it exists). *.Example I.1 :......................................................... * Procedure computing a physical property depending on the temperature. 'DEBP' CALCUL ; 'ARGU' RX*TABLE ; iarg = rx . 'IARG' ; rv = rx . 'EQEX' ; * Input for the CALCUL procedure 'SI' ( 'NON' ( 'EGA' iarg 1)) ; 'MESS' 'Procedure CALCUL : nombre d arguments incorrect ' iarg ; 'QUIT' CALCUL ; 'FINSI' ; 'SI' ( 'EGA' ('TYPE' rx . 'ARG1') 'MOT ') ; TN = rv . 'INCO' . (rx . 'ARG1') ; 'SINON' ; 'MESS' 'Procedure CALCUL : type argument invalide ' ; 'QUIT' CALCUL ; 'FINSI' ; * The temperature is given in Kelvin T = TN + 273. ; *Viscosity : Sutherland law : Kg/m/s MU = 1.648*(T**1.5) * ('INVE' (T + 0.648)) *Conductivity : Sutherland law : W/m/oC LB = 1.368*(T**1.5) * ('INVE' (T + 0.368)) *J/kg/oC CP = 1015 ; * Prandtl number Pr = MU * CP * ('INVE' LB) * The Reynolds number is given by : Re = 400. ; Pe = Re * Pr ; rv . 'INCO' . 'IPE' = 'INVE' Pe ; * The last instruction creates void objects to satisfy the EXEC procedur as2 ama1 = 'KOPS' 'MATRIK' ; 'FINPROC' as2 ama1 ; * We create a RV table describing the physical problem. * We choose an explicit algorithm (OPTI 'EFM1') * The time step is automaticaly computed * (key word 'DELTAT' for DFDT) * We will perform 200 time steps. * The Peclet number is computed in the CALCUL procedur. RV = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' 200 'ZONE' $mt 'OPER' CALCUL 'TN' 'OPTI' 'EFM1' 'SUPG' 'ZONE' $mt 'OPER' 'TSCA' 'IPE' 'UN' 0. 'INCO' 'TN' 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1. 'TN' 'DELTAT' 'INCO' 'TN' 'CLIM' 'TN' 'TIMP' entree 0. 'TN' 'TIMP' paroi 1. ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' .'UN'= 'KCHT' $mt 'VECT' 'SOMMET' (1. 0.) ; rv . 'INCO' .'TN'= 'KCHT' $mt 'SCAL' 'SOMMET' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example I.1 ...................................................... I.2/ Direct computation of a steady state. - We can try to find directly the steady state (if it exists) with internal iterations to solve the non linearities. *.Example I.2 :......................................................... * We create a RV table describing the physical problem. * We chooze an implicit algorithm (OPTI 'EF' 'IMPL') * We will perform 10 iterations with a relaxation coefficient OMEGA=0.5 * The Peclet number is computed in the CALCUL procedur as in example 1. RV = 'EQEX' 'OMEGA' 0.5 'NITER' 10 'ITMA' 0 'ZONE' $mt 'OPER' CALCUL 'TN' 'OPTI' 'EF' 'SUPG' 'IMPL' 'ZONE' $mt 'OPER' 'TSCA' 'IPE' 'UN' 0. 'INCO' 'TN' 'CLIM' 'TN' 'TIMP' entree 0. 'TN' 'TIMP' paroi 1. ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' . 'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (1. 0.) ; rv . 'INCO' . 'TN' = 'KCHT' $mt 'SCAL' 'SOMMET' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example I.2 ...................................................... I.3/ Implicit transient computation - We can perform an implicit transient computation with or without internal iterations for each time step. *.Example I.3 :......................................................... * We create a RV table describing the physical problem. * We chooze an implicit algorithm (OPTI 'EF' 'IMPL') first order in time * or better a second order in time (Crank Nicolson). * ('OPTI' 'EF' 'SUPG' 'SEMI' 0.5) * We will perform 10 time steps without internal iterations. * The Peclet number is computed in the CALCUL procedur as in example 1. dt = 1. ; RV = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' 10 'ZONE' $mt 'OPER' CALCUL 'TN' 'OPTI' 'EF' 'SUPG' 'SEMI' 0.5 'ZONE' $mt 'OPER' 'TSCA' 'IPE' 'UN' 0. 'INCO' 'TN' 'CLIM' 'TN' 'TIMP' entree 0. 'TN' 'TIMP' paroi 1. 'OPTI' 'EF' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1. 'TN' dt 'INCO' 'TN' ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' . 'UN'= 'KCHT' $mt 'VECT' 'SOMMET' (1. 0.) ; rv . 'INCO' . 'TN'= 'KCHT' $mt 'SCAL' 'SOMMET' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example I.3 ...................................................... II/ INCOMPRESSIBLE NAVIER STOKES The procedure performs the resolution of the Navier_Stokes equations by a finite element method (FE) using primitive variables (pressure - velocity), for an incompressible flow, or slightly compressible, for a transient regime or a steady state. Several boundary conditions can be used (prescribed velocity, prescribed pressure, momentum source pressure drop ... etc. The system can be linear or not. The non linearities are solved by a fixed point method. The description of the equation to be solved is done by the operator EQEX which creates a table TAB1. Three algorithms are availlable to solve the velocity-pressure system. - A semi-explicit algorithm : implicit for the pressure, explicit for the velocity and eventually for all other convected scalar quantities. (Gresho) - An implicit algorithm : Direct resolution of the velocity-pressure system (Taylor-Hood). This algorithm allows a direct search of the steady state. - A projection method : This algorithm splits the resolution into two steps, first the convection/diffusion is solved, then the velocity field is projected on divergence free space (Chorin Temam). The parameters for the algorithm are defined in the table TAB1. The non linearities, in particular on the physical properties can be solved as previously mentioned, by a CALCUL procedur, called at each time step and/or at each iteration. The steady state can be obtained as the asymptotic limit of the transient. II.1/ Semi explicit computation. - One can perform a transient semi-explicit computation. The pressure is implicit and the velocity explicit. The time step is limited by a stability condition (CFL) or Fourier condition. The time step can be automatically computed (Key word 'DELTAT' as third argument of DFDT) or prescribed (instead of the keyword 'DELTAT' give the value). Generally the algorithm needs an important amount of time steps. The construction of the TAB1 table as to be done into two steps. - First we describe the explicit equations. - Then we construct a second table using again the EQEX operator which describes the operator linked to the pressure. The normal or tangential velocities conditions are included in this case because they are treated by Lagrange multipliers. In the contrary with the FE formulation (weak form) a precribed pressure is linked to the momentum equation. - At last the second table is placed at the 'POISSON' entry of the first table. The entry rv.'CALPRE' = VRAI of the first table indicates that the pressure matrix is computed at each time step. This is necessary with the A.L.E. formulation. rv.'CALPRE' = FAUX or absence of that entry leads to a computation once for all of the pressure matrix. The positionning of the variable rv.'DETMAT' to VRAI indicates that the MATRIK objects will be droped at the end of the procedure. *.Example II.1 :........................................................... * Lid driven square cavity * !!!! WARNING : The cavity is closed (V.n prescribed all around the boundary) * Then it is necessary to impose the pressure at a point. ro = 400. ; mu = 1. ; RV = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' 500 'ALFA' 0.5 'OPTI' 'EFM1' 'SUPG' 'ZONE' $mt 'OPER' 'NS' (mu/ro) 'INCO' 'UN' 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN' ; RV = 'EQEX' RV 'CLIM' 'UN' 'UIMP' CD 1. 'UN' 'VIMP' CD 0. 'UN' 'UIMP' DA 0. 'UN' 'VIMP' DA 0. 'UN' 'UIMP' AB 0. 'UN' 'VIMP' AB 0. 'UN' 'UIMP' BC 0. 'UN' 'VIMP' BC 0. ; * The choice of an iterative method (Bicg stab + preconditioning MILU0) * allow greater meshes (See INFO KRES ; ) rv. 'METHINV' . 'TYPINV' = 3 ; rv. 'METHINV' . 'IMPINV' = 0 ; rv. 'METHINV' . 'NITMAX' = 400 ; rv. 'METHINV' . 'PRECOND' = 3 ; rv. 'METHINV' . 'RESID' = 1.e-8 ; rv. 'METHINV' . 'FCPRECT' = 1 ; rv. 'METHINV' . 'FCPRECI' = 1 ; betastab=1.e2 ; RVP = 'EQEX' 'OPTI' 'EF' 'CENTRE' 'ZONE' $mt 'OPER' 'KBBT' -1. betastab 'INCO' 'UN' 'PRES' 'CLIM' 'PRES' 'TIMP' bcp 0. ; rvp . 'METHINV' . 'TYPINV' = 2 ; rvp . 'METHINV' . 'IMPINV' = 0 ; rvp . 'METHINV' . 'NITMAX' = 300 ; rvp . 'METHINV' . 'PRECOND' = 3 ; rvp . 'METHINV' . 'RESID' = 1.e-8 ; rvp . 'METHINV' . 'FCPRECT' = 100 ; rvp . 'METHINV' . 'FCPRECI' = 100 ; rv . 'POISSON' = rvp ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' .'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ; rv . 'INCO' .'PRES' = 'KCHT' $mt 'SCAL' 'CENTRE' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example II.1 ..................................................... II.2/ Implicit computation. - One can perform an implicit transient computation with or without internal iterations at each time step. *.Example II.2 :........................................................... * Lid driven square cavity * !!!! WARNING : The cavity is closed (V.n prescribed all around the boundary) * Then it is necessary to impose the pressure at a point. ro = 400. ; mu = 1. ; dt = 5. ; RV = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' 20 'OPTI' 'EF' 'IMPL' 'SUPG' 'ZONE $mt 'OPER' 'LAPN' mu 'INCO' 'UN' 'ZONE $mt 'OPER' 'KONV' ro 'UN' mu dt 'INCO' 'UN' 'OPTI' 'EF' 'CENTREE' 'ZONE $mt 'OPER' 'DFDT' ro 'UN' dt 'INCO' 'UN' 'OPTI' 'EF' 'CENTREP1' 'ZONE $mt 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PRES' ; RV = 'EQEX' RV 'CLIM' 'PRES' 'TIMP' bcp 0. 'UN' 'UIMP' CD 1. 'UN' 'VIMP' CD 0. 'UN' 'UIMP' DA 0. 'UN' 'VIMP' DA 0. 'UN' 'UIMP' AB 0. 'UN' 'VIMP' AB 0. 'UN' 'UIMP' BC 0. 'UN' 'VIMP' BC 0. ; * The choice of an iterative method (Bicg stab + preconditioning MILU0) * allow greater meshes (See INFO KRES ; ) rv . 'METHINV' . 'TYPINV' = 3 ; rv . 'METHINV' . 'IMPINV' = 0 ; rv . 'METHINV' . 'NITMAX' = 400; rv . 'METHINV' . 'PRECOND' = 3 ; rv . 'METHINV' . 'RESID' = 1.e-8 ; rv . 'METHINV' . 'FCPRECT' = 1 ; rv . 'METHINV' . 'FCPRECI' = 1 ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' . 'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ; rv . 'INCO' . 'PRES' = 'KCHT' $mt 'SCAL' 'CENTREP1' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example II.2 ..................................................... II.3/ Direct computation of a steady state - One can perform a direct computation of a steady state with internal iterations to solve the non linearities. *.Example II.3 :........................................................... * Lid driven square cavity * !!!! WARNING : The cavity is closed (V.n prescribed all around the boundary) * Then it is necessary to impose the pressure at a point. ro = 400. ; mu = 1. ; RV = 'EQEX' 'OMEGA' 0.7 'NITER' 10 'ITMA' 0 'OPTI' 'EF' 'IMPL' 'SUPG' 'ZONE' $mt 'OPER' 'LAPN' mu 'INCO' 'UN' 'ZONE' $mt 'OPER' 'KONV' ro 'UN' mu 'INCO' 'UN' 'OPTI' 'EF' 'CENTREP1' 'ZONE' $mt 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PRES' ; RV = 'EQEX' RV 'CLIM' 'PRES' 'TIMP' bcp 0. 'UN' 'UIMP' CD 1. 'UN' 'VIMP' CD 0. 'UN' 'UIMP' DA 0. 'UN' 'VIMP' DA 0. 'UN' 'UIMP' AB 0. 'UN' 'VIMP' AB 0. 'UN' 'UIMP' BC 0. 'UN' 'VIMP' BC 0. ; * The choice of an iterative method (Bicg stab + preconditioning MILU0) * allow greater meshes (See INFO KRES ; ) rv . 'METHINV' . 'TYPINV' = 3 ; rv . 'METHINV' . 'IMPINV' = 0 ; rv . 'METHINV' . 'NITMAX' = 400; rv . 'METHINV' . 'PRECOND' = 3 ; rv . 'METHINV' . 'RESID' = 1.e-8 ; rv . 'METHINV' . 'FCPRECT' = 1 ; rv . 'METHINV' . 'FCPRECI' = 1 ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' .'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ; rv . 'INCO' .'PRES' = 'KCHT' $mt 'SCAL' 'CENTREP1' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example II.3 ..................................................... II.4/ Projection Method - One can perform an implicit transient calculation where the solution of the momentum equation and of the continuity equation are split. The time step is no longer limited by a stability condition like a CFL condition. However in practice it must be choosen to obtain the desired accuracy. The algorithm needs much less time steps than the explicit one. It is the most economical for a transient computation even for a steady state. - In practice it looks like the semi explicit one except the operators can be implicit. The construction of the TAB1 table as to be done into two steps. - First we describe the explicit equations. - Then we construct a second table using again the EQEX operator which describes the operator linked to the pressure. The normal or tangential velocities conditions are included in this case because they are treated by Lagrange multipliers. In the contrary with the FE formulation (weak form) a precribed pressure is linked to the momentum equation. - At last the second table is placed at the 'PROJ' entry of the first table. The entry rv.'CALPRE' = VRAI of the first table indicates that the pressure matrix is computed at each time step. This is necessary with the A.L.E. formulation. rv.'CALPRE' = FAUX or absence of that entry leads to a computation once for all of the pressure matrix. *.Example II.4 :........................................................... * Lid driven square cavity * !!!! WARNING : The cavity is closed (V.n prescribed all around the boundary) * Then it is necessary to impose the pressure at a point. ro = 400. ; mu = 1. ; dt = 1. ; RV = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' 500 'ALFA' 0.5 'OPTI' 'EF' 'IMPL' 'SUPG' 'ZONE' $mt 'OPER' 'NS' (mu/ro) 'INCO' 'UN' 'OPTI' 'EFM1' 'CENTREE' 'ZONE' $mt 'OPER' 'DFDT' 1. 'UN' dt 'INCO' 'UN' 'CLIM' 'UN' 'UIMP' CD 1. 'UN' 'VIMP' CD 0. 'UN' 'UIMP' DA 0. 'UN' 'VIMP' DA 0. 'UN' 'UIMP' AB 0. 'UN' 'VIMP' AB 0. 'UN' 'UIMP' BC 0. 'UN' 'VIMP' BC 0. ; * The choice of an iterative method (Bicg stab + preconditioning MILU0) * allow greater meshes (See INFO KRES ; ) rv . 'METHINV' . 'TYPINV' = 3 ; rv . 'METHINV' . 'IMPINV' = 0 ; rv . 'METHINV' . 'NITMAX' = 400 ; rv . 'METHINV' . 'PRECOND' = 3 ; rv . 'METHINV' . 'RESID' = 1.e-8 ; rv . 'METHINV' . 'FCPRECT' = 1 ; rv . 'METHINV' . 'FCPRECI' = 1 ; betastab=1.e2 ; RVP = 'EQEX' 'OPTI' 'EF' 'CENTRE' 'ZONE' $mt 'OPER' 'KBBT' -1. betastab 'INCO' 'UN' 'PRES' 'CLIM' 'PRES' 'TIMP' bcp 0. ; rvp . 'METHINV' . 'TYPINV' = 2 ; rvp . 'METHINV' . 'IMPINV' = 0 ; rvp . 'METHINV' . 'NITMAX' = 300; rvp . 'METHINV' . 'PRECOND' = 3 ; rvp . 'METHINV' . 'RESID' = 1.e-8 ; rvp . 'METHINV' . 'FCPRECT' = 100 ; rvp . 'METHINV' . 'FCPRECI' = 100 ; rv . 'PROJ' = RVP ; rv . 'INCO' = 'TABLE' 'INCO' ; rv . 'INCO' . 'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ; rv . 'INCO' . 'PRES' = 'KCHT' $mt 'SCAL' 'CENTRE' 0. ; EXEC RV ; * The computed fields are stored in table rv . 'INCO' *.End example II.4 ..................................................... Remarks : _________ 1) The logic variable rv.'DETMAT' set to VRAI indicates that the MATRIK objects will be destroyed at the end of the procedure. 2) Two LOGIQUE type variables resp. rv.'STOPITER' and rv.'STOPPDT', initially set to FAUX, can be set to VRAI (by a user-defined procedure for example) in order to stop resp. the non-linearity solution loop and the time-stepping loop. rv . 'NUITER' is the value of the current non-linear iteration. 3) It is possible to use an algebraic incremental projection method to approximately solve the incompressible NS system while keeping the same syntax as in the implicit transient case. One only needs to add a table at the 'GPROJ' index in rv : rv . 'GPROJ' = 'TABLE' ; rv . 'GPROJ' . 'NOMVIT' = 'UN' ; rv . 'GPROJ' . 'NOMPRES' = 'CHAINE' 'PRES' ; where 'UN' is the speed unknown name and 'PRES' the pressure unknown name. It is also possible to give a table at the index : rv . 'GPROJ' . 'METHINV' to specify options for the pressure matrix linear solver. (cf. KRES documentation) One can also specify if a simple or double projection method is wanted : rv . 'GPROJ' . 'dblproj' = FAUX ou VRAI (default value : VRAI)
© Cast3M 2003 - All rights reserved.
Disclaimer