Afficher cette notice en
$$$$ 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 ; 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)
© Cast3M 2003 - Tous droits réservés.
Mentions légales