$$$$ EXEC NOTICE GOUNAND 22/11/22 21:15:01 11504 DATE 22/11/22 Procedure EXEC Voir aussi : EQEX MODE -------------- DOMA KCHT KRES EXEC TAB1 ; Section : Fluides Resolution FRAN=================================================================== Objet : _______ I/ TRANSPORT/DIFFUSION D'UN SCALAIRE Cette procedure permet d'effectuer un calcul du transport (convection/diffusion), d'un scalaire passif, en transitoire ou en regime permanent, en presence de conditions aux limites variees, valeur imposee, flux, echange avec le milieu exterieur, source etc. Le calcul peut etre lineaire ou non. Les non linearites sont resolues par une methode de point fixe. La description de l'equation a resoudre se fait a l'aide de l'operateur EQEX qui cree la table TAB1. Le scalaire peut aussi bien etre la temperature qu'une concentration ou toute grandeur physique intensive. Les parametres de l'algorithme sont definis dans la table TAB1 a l'aide de EQEX. I.1/ Calcul transitoire explicite. Le pas de temps est soumis a une contrainte de stabilite. Il peut etre impose ou calcule automatiquement. Les non linearites, notamment sur les proprietes physiques,peuvent etre resolues en les actualisant dans une procedure de calcul appelee a chaque pas de temps. Le regime permanent peut etre obtenu comme limite asymptotique du transitoire. *.Exemple I.1 :......................................................... * Procedure calculant une propriete physique dependant de la temperature 'DEBP' CALCUL ; 'ARGU' RX*TABLE ; iarg = rx . 'IARG' ; rv = rx . 'EQEX' ; *Lecture des arguments de l'operateur calcul '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' ; * La temperature est exprimee en Kelvin T = TN + 273. ; *Viscosite dynamique : loi de Sutherland : Kg/m/s MU = 1.648*(T**1.5) * ('INVE' (T + 0.648)) *conductivite : loi de Sutherland : W/m/oC LB = 1.368*(T**1.5) * ('INVE' (T + 0.368)) *J/kg/oC CP = 1015 ; * Nombre de Prandtl Pr = MU * CP * ('INVE' LB) * le Reynolds est donne Re = 400. ; Pe = Re * Pr ; rv . 'INCO' . 'IPE' = 'INVE' Pe ; * La derniere instruction cree des objet vides * pour satisfaire la procedure EXEC as2 ama1 = 'KOPS' 'MATRIK' ; 'FINPROC' as2 ama1 ; * On cree la table RV decrivant le probleme physique * On choisit un algorithme explicite (OPTI 'EFM1') * Le pas de temps est calcule automatiquement * (mot cle 'DELTAT' sur DFDT) * On fera 200 pas de temps * Le Peclet est calcule dans la procedure CALCUL 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple I.1 ...................................................... I.2/ Calcul direct d'un regime permanent. - On peut faire une recherche directe d'un regime permanent avec des iterations internes pour resoudre les non-linearites. *.Exemple I.2 :......................................................... * On cree la table RV decrivant le probleme physique * On choisit un algorithme IMPLICITE (OPTI 'EF' 'IMPL') * On fera 10 iterations avec un facteur de relaxation de OMEGA=0.5 * Le Peclet est calcule dans la procedure CALCUL decrite dans * l'exemple 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple I.2 ...................................................... I.3/ Calcul transitoire implicite. - On peut faire un calcul transitoire implicite avec ou sans iterations internes a chaque pas de temps. *.Exemple I.3 :......................................................... * On cree la table RV decrivant le probleme physique * On choisit un algorithme IMPLICITE (OPTI 'EF' 'IMPL') ordre 1 en * temps ou mieux un schema en temps d'ordre 2 (Crank Nicolson) * ('OPTI' 'EF' 'SUPG' 'SEMI' 0.5) * On fera 10 pas de temps sans iteration interne * Le Peclet est calcule dans la procedure CALCUL decrite dans * l'exemple 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple I.3 ...................................................... II/ NAVIER STOKES INCOMPRESSIBLE La procedure permet de resoudre les equations de Navier_Stokes par une methode d'elements finis (EF) en variables primitives (vitesse et pression), pour un ecoulement incompressible, ou faiblement compre- ssible, en regime permanent ou en transitoire. Les conditions limites peuvent etre variees, vitesse imposee, pression imposee, source de quantite de mouvement, perte de charge ... etc. Le calcul peut etre lineaire ou non. Les non linearites sont resolues par une methode de point fixe. La description de l'equation a resoudre se fait a l'aide de l'operateur EQEX qui cree la table TAB1. Trois algorithmes sont proposes pour resoudre le systeme vitesse- pression. - Un algorithme semi-explicite : implicite pour la pression, explicite pour la vitesse et eventuellement toutes les autres quantites scalaires transportees (a la Gresho). - Un algorithme implicite : resolution directe du systeme en vitesse- pression (a la Taylor-Hood). Cet algorithme permet aussi de faire une recherche directe d'un regime permanent. - Une methode de projection : cet algorithme consiste en deux etapes : une premiere etape de transport diffusion et une seconde de projection sur un espace a divergence nulle (a la Chorin Temam) . Les parametres de l'algorithme sont definis dans la table TAB1 Les non linearites, notamment sur les proprietes physiques,peuvent etre resolues, comme precedemment, par l'intermediaire d'une methode de point fixe, en actualisant dans une procedure les non linearites. de calcul appelee a chaque pas de temps ou chaque iteration. Le regime permanent peut etre obtenu comme limite asymptotique du transitoire II.1/ Calcul semi-explicite - On peut faire un calcul transitoire semi explicite. La pression est implicite, la vitesse explicite. De ce fait le pas de temps est soumis a une contrainte de stabilite de type CFL liee a la convection ou a la diffusion. Le pas de temps peut etre impose ou calcule automatiquement (mot cle 'DELTAT' en 3eme argument de DFDT). L'algorithme necessite en general un grand nombre de pas de temps. La construction de la table se fait en deux etapes. - Une premiere etape consiste a creer une table decrivant la partie explicite. - Dans la deuxieme etape on construit la table decrivant les operateurs lies a la pression. Les conditions de vitesse normale ou vitesse tangente en font partie. Elles sont traitees a l'aide de multiplicateurs de Lagrange. Par contre dans nos formulations elements finis (formulation faible) la condition de pression imposee n'en fait pas partie. La pression est imposee comme les contraintes visqueuses sur l'equation de quantite de mouvement. - Enfin on place la deuxieme table a l'indice 'POISSON' de la premiere. L'indice rv.'CALPRE' = VRAI de la premiere table indique que l'on recalcule a chaque pas de temps la matrice de pression. Ceci est necessaire en formulation A.L.E. rv.'CALPRE' = FAUX ou l'absence de cet indice fait que la matrice de pression sera calculee une fois pour toute. Le positionnement de la variable rv.'DETMAT' a VRAI indique que les objets MATRIK presents dans la table rv seront supprimes. *.Exemple II.1 :........................................................... * Cas de la cavite carree a paroi defilante * !!!! ATTENTION : La cavite etant fermee (V.n impose sur toute la frontiere) * il est necessaire d'imposer la pression en un 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. ; * Le choix d'une methode iterative (Bicg stab + preconditionement MILU0) * permet de passer des cas un peu plus gros. (Faire 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple II.1 ..................................................... II.2/ Calcul implicite. - On peut faire un calcul transitoire implicite avec ou sans iterations internes a chaque pas de temps. *.Exemple II.2 :........................................................... * Cas de la cavite carree a paroi defilante * !!!! ATTENTION : La cavite etant fermee (V.n impose sur toute la frontiere) * il est necessaire d'imposer la pression en un 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. ; * Le choix d'une methode iterative (Bicg stab + preconditionement MILU0) * permet de passer des cas un peu plus gros. (Faire 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple II.2 ..................................................... II.3/ Calcul direct d'un regime permanent. - On peut faire une recherche directe d'un regime permanent avec des iterations internes pour resoudre les non-linearites. *.Exemple II.3 :........................................................... * Cas de la cavite carree a paroi defilante * !!!! ATTENTION : La cavite etant fermee (V.n impose sur toute la frontiere) * il est necessaire d'imposer la pression en un 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. ; * Le choix d'une methode iterative (Bicg stab + preconditionement MILU0) * permet de passer des cas un peu plus gros. (Faire 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple II.3 ..................................................... II.4/ Methode de projection - On peut faire un calcul transitoire implicite en decouplant la reso- lution de l'equation de quantite de mouvement de la resolution de l'eq- uation de continuite. Le pas de temps n'est plus soumis a une contrainte de stabilite de type CFL. Cependant en pratique on ne peut pas prendre un pas de temps arbitrairement grand. L'algorithme necessite beaucoup moins de pas de temps que l'algorithme semi explicite. A l'usage il s'avere etre le plus economique pour un calcul transitoire voire meme un calcul permanent. - La mise en oeuvre est assez similaire a la difference que les opera- teurs peuvent etre implicites. La construction de la table se fait en deux etapes. - Une premiere etape consiste a creer une table decrivant l'equation de quantite de mouvement. - Dans la deuxieme etape on construit la table decrivant les operateurs lies a la projection (div U = 0) Les conditions de vitesse normale ou vitesse tangente en font partie, elles sont traitees a l'aide de multiplicateurs de Lagrange. Comme dans l'algorithme semi explicite la condition de pression imposee n'en fait pas partie. La pression est imposee comme les contraintes visqueuses sur l'equation de quantite de mouvement. - Enfin on place la deuxieme table a l'indice 'PROJ' de la premiere. L'indice rv.'CALPRE' = VRAI de la premiere table indique que l'on recalcule a chaque pas de temps la matrice de pression. Ceci est necessaire en formulation A.L.E. rv.'CALPRE' = FAUX ou l'absence de cet indice fait que la matrice de pression sera calculee une fois pour toute. *.Exemple II.4 :........................................................... * Cas de la cavite carree a paroi defilante * !!!! ATTENTION : La cavite etant fermee (V.n impose sur toute la frontiere) * il est necessaire d'imposer la pression en un 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. ; * Le choix d'une methode iterative (Bicg stab + preconditionement MILU0) * permet de passer des cas un peu plus gros. (Faire 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 ; * Les resultats se trouvent dans la table rv . 'INCO' *.Fin exemple II.4 ..................................................... Remarques : ___________ 1) Le positionnement de la variable rv.'DETMAT' a VRAI indique que les objets MATRIK presents dans la table rv seront supprimes. 2) Deux variables de type LOGIQUE resp. rv.'STOPITER' et rv.'STOPPDT' initialisees a FAUX peuvent etre mises a VRAI (par exemple a l'aide d'une procedure utilisateur) afin de stopper resp. la boucle de resolution de la non-linearite et la boucle sur les pas de temps. L'indice rv . 'NUITER' contient le numero de l'iteration non lineaire courante. 3) Il est possible d'utiliser une methode de projection algebrique incrementale pour resoudre le systeme NS incompressible de maniere approchee en gardant la meme syntaxe qu'en calcul transitoire implicite et en rajoutant une table a l'indice 'GPROJ' : rv . 'GPROJ' = 'TABLE' ; rv . 'GPROJ' . 'NOMVIT' = 'UN' ; rv . 'GPROJ' . 'NOMPRES' = 'CHAINE' 'PRES' ; ou 'UN' est le nom de l'inconnue vitesse et 'PRES' le nom de l'inconnue pression. Il est egalement possible de donner une table a l'indice : rv . 'GPROJ' . 'METHINV' precisant les options pour l'inversion de la matrice de pression. (cf. notice KRES) On peut également préciser si on veut une méthode de simple ou double projection : rv . 'GPROJ' . 'dblproj' = FAUX ou VRAI (VRAI par défault) ANGL=================================================================== 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. Setting the variable rv.'DETMAT' to VRAI indicates that the MATRIK objects in the rv TABLE will be deleted. *.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) Setting the logical variable rv.'DETMAT' to VRAI indicates that the MATRIK objects in the rv TABLE will be deleted. 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)