* fichier : back_proj_3.dgibi **************** CAS TEST : back_proj_3.dgibi ******************** * * Ce test permet de vérifier le bon fonctionnement des opérateurs * utilisés pour la résolution des équations de NAVIER_STOKES en EF * par un algorithme de projection. * * Le cas étudié est celui d'un écoulement laminaire dans un canal * en présence d'une marche descendante. On teste la position du * point de réattachement par rapport au pied de la marche. * * Référence : K. Morgan, J. Périaux anf F. Thomasset, editors, * Analysis of laminar flow over a backward facing step, Vol9 of * Notes on numerical Fluid Mechanics, Vieweg, 1984. * *------------------------------------------------------------------- * Auteur : F.Dabbene (LTMF) 03/03 *------------------------------------------------------------------- * 'SAUT' 'PAGE' ; 'OPTION' 'DIME' 2 'ELEM' 'QUA8' 'ECHO' 0 ; * *- Pilotage du calcul * * DISCR : Element en vitesse (LINE/MACRO/QUAF) * KPRES : Element en pression (CENTRE/CENTREP2/MSOMMET) * KSUPG : Méthode de décentrement (CENTREE/SUPG/SUPGDC) * KMASS : Matrice masse lumpée ou non (EFM1/EF) * * DT : Valeur du pas de temps * (En cas de problème de convergence, augmenter DT et/ou * prendre KMASS en EFM1) * EPS0 : Test d'arret * RAF : Taille de maille définie par 0.1*raf * * FREQ0 : Fréquence d'évaluation du résidu * GRAPH : Booleen pour les tracés * COMPLET : Booleen volume de calcul * DISCR = 'MACRO' ; KPRES = 'CENTRE' ; KSUPG = 'SUPG' ; KMASS = 'EFM1' ; GRAPH = FAUX ; COMPLET = FAUX ; 'SI' COMPLET ; DT = 3. ; FREQ0 = 20 ; EPS0 = 1.D-5 ; RAF = 1. ; 'SINO' ; DT = 3. ; FREQ0 = 20 ; EPS0 = 1.D-2 ; RAF = 3. ; 'FINSI' ; * *--------------------------------------------------------------------------- * Recherche du point de réattachement (point où dUx/dy=0) *--------------------------------------------------------------------------- * 1/ Après avoir calculé le gradient de Ux, on ne conserve que les valeurs * sur la frontière qui nous intéresse sous la forme d'une évolution. * 2/ On borne l'évolution au voisinage du point de réattachement afin * d'avoir une variation monotone sur l'intervale de dUx/dy. * 3/ On recherche par interpolation le zero de la fonction (dUx/dy(s),s) * La valeur obtenu est l'abscisse curviligne cherchée * * Remarques : * (i) 1/ permettrait de calculer le coeff de frottement à la paroi : * il suffirait de diviser EV1 par le bon coefficient (Re/2 ici) * (ii) Il est impératif que dUx/dy soit monotone sur l'intervale 2/ afin * qu'il y ait unicité du zero (principe des valeurs intermédiaires) *--------------------------------------------------------------------------- * 'DEBPROC' attac ; * 1/ Ux = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UX' 0. RV.'INCO'.'UN' ; DUxDY = 'EXCO' 'UY' ('KOPS' Ux 'GRADS' $DOMTOT) 'SCAL' ; EV1 = 'EVOL' 'CHPO' DUxDY BOTTOM ; * 2/ EV2 = 'EXTR' EV1 'APRE' 5. ; EV3 = 'EXTR' EV2 'AVAN' 10. ; * 3/ LX1 = 'EXTR' EV3 'ABSC' ; LY1 = 'EXTR' EV3 'ORDO' ; LY1 = 'ORDO' LY1 ; Ymin = 'MINI' LY1 ; Ymax = 'MAXI' LY1 ; Delta = Ymin * Ymax ; 'SI' (Delta < 0.) ; Y0 = 0. ; 'SINON' ; 'MESS' 'Fonction 'NON' monotone' ; Y0 = Ymin ; * 'ERRE''Fonction 'NON' monotone' ; 'FINSI' ; V1 = 'IPOL' Y0 LY1 LX1 ; 'FINP' V1 ; * *================================================================== * Calcul du résidu basé sur la composante horizontale de la vitesse * et arrêt suivant un critère transmis *================================================================== * E/ : RVX : TABLE : TABLE des données créées par EQEX * ARG1 : Fréquence d'impression * ARG2 : Critère d'arrêt * ARG3 : Booleen de tracé * /S : MAT1 : MATRIK : Objet vide * /S : CHP1 : CHPO : Objet vide *------------------------------------------------------------------ * MAT1 et CHP1 permettent d'assurer la compatibilité des opérateurs * de discrétisation avec les procédures personnelles *------------------------------------------------------------------ 'DEBPROC' residu rvx*table ; RV = RVX . 'EQEX' ; FREQ = RVX . 'ARG1' ; EPS0 = RVX . 'ARG2' ; GRAPH = RVX . 'ARG3' ; NITER = RV . 'NITER' ; DD = RV . 'PASDETPS' . 'NUPASDT' ; NN = DD '/' FREQ ; L0 = 'EGA' (DD '-' (FREQ*NN)) 0 ; 'SI' L0 ; RANG0 = RV . 'PASDETPS' . 'NUPASDT' ; TIME0 = RV . 'PASDETPS' . 'TPS' ; UN0 = 'EXCO' 'UX' RV . 'INCO' . 'UN' 'SCAL' ; UNM0 = 'EXCO' 'UX' RV . 'INCO' . 'UN2' 'SCAL' ; ERR0 = ('MAXIMUM' ('ABS' (UN0 '-' UNM0))) '+' 1.D-20 ; ERR10 = ('LOG' ERR0 ) '/' ('LOG' 10.) ; 'MESSAGE' 'Résidu en vitesse suivant X au pas' RANG0 '(t=' TIME0 ') :' ERR0 ':' ERR10 ; RV . 'INCO' . 'IT' = RV . 'INCO' . 'IT' 'ET' ('PROG' RANG0) ; RV . 'INCO' . 'TI' = RV . 'INCO' . 'TI' 'ET' ('PROG' TIME0) ; RV . 'INCO' . 'ER' = RV . 'INCO' . 'ER' 'ET' ('PROG' ERR10) ; V1 = attac ; RV . 'INCO' . 'POSI' = RV . 'INCO' . 'POSI' 'ET' ('PROG' V1) ; Y1 = ('LOG' EPS0) '/' ('LOG' 10) ; Y2 = 0. ; 'SI' GRAPH ; EV1 = 'EVOL' 'MANU' (RV . 'INCO' . 'IT')(RV . 'INCO' . 'ER') ; 'DESSIN' EV1 'YBOR' Y1 Y2 'NCLK' ; 'FINSI' ; 'SI' ((ERR10 < Y1) 'ET' (DD > ('MAXI' ('LECT' 10 FREQ)))) ; RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ; 'FINSI' ; 'FINSI' ; RV . 'INCO' . 'UN2' = 'COPIER' RV . 'INCO' . 'UN' ; mat1 chp1 = 'KOPS' 'MATRIK' ; 'FINP' mat1 chp1 ; * *========== * Maillage *========== * * *----------------------------------------------------------------------- * * * < L1 X L2 - L1 > * * 7___________6__________________________________________________5 * ^ ^ * INLET H1 * v H2 OUTLET * 1___________2 * | v * 3__________________________________________________4 * BOTTOM * * *----------------------------------------------------------------------- * * * L1 : Longueur de la section d'entrée (avant la marche) * L2 : Longueur de la totalité du dispositif * H1 : Hauteur de la section d'entrée * H2 : Hauteur de la section de sortie * d1 : Dimension caractéristique d'une maille L1 = 3.0 ; L2 = 22.0 ; H1 = 1.0 ; H2 = 1.5 ; d1 = 0.1 * raf ; * * H2-H1 : Hauteur de la marche servant à l'adimensionnalisation HDIM = H2 - H1 ; L1 = L1 / HDIM ; L2 = L2 / HDIM ; H1 = H1 / HDIM ; H2 = H2 / HDIM ; d1 = d1 / HDIM ; * * Points du maillage p1 = 0. (H2-H1) ; p2 = L1 (H2-H1) ; p3 = L1 0. ; p4 = L2 0. ; p5 = L2 H2 ; p6 = L1 H2 ; p7 = 0. H2 ; p8 = L2 (H2-H1) ; * * Section d'éntrée p1p2 = p1 'DROIT' 'DINI' d1 'DFIN' d1 p2 ; p2p6 = p2 'DROIT' 'DINI' d1 'DFIN' d1 p6 ; p6p7 = p6 'DROIT' 'DINI' d1 'DFIN' d1 p7 ; p7p1 = p7 'DROIT' 'DINI' d1 'DFIN' d1 p1 ; mesh1 = 'DALL' p1p2 p2p6 p6p7 p7p1 'PLAN' ; * * Section de sortie p6p2 = 'INVE' p2p6 ; p2p3 = p2 'DROIT' 'DINI' d1 'DFIN' d1 p3 ; p3p4 = p3 'DROIT' 'DINI' d1 'DFIN' d1 p4 ; p4p8 = p4 'DROIT' 'DINI' d1 'DFIN' d1 p8 ; p8p5 = p8 'DROIT' 'DINI' d1 'DFIN' d1 p5 ; p5p6 = p5 'DROIT' 'DINI' d1 'DFIN' d1 p6 ; mesh2 = 'DALL' (p6p2 'ET' p2p3) p3p4 (p4p8 'ET' p8p5) p5p6 'PLAN' ; * * *======================= * Modèles NAVIER_STOKES *======================= * * * Définition des équations * $DOMTOT : Modèle volumique défini sur le maillage complet * * Conditions aux limites en vitesse : * $INLET : Modèle surfacique défini à l'entrée fluide (Poiseuille) * $WALL : Modèle surfacique défini sur les murs (adhérence en paroi) * * Conditions aux limites en pression : * $OUTLET : Modèle surfacique défini à la sortie fluide (sert à imposer * la pression de sortie en pression continue MSOMMET) * * Post-traitement * $BOTTOM : Modèle surfacique défini sur le plancher après la marche * (sert à évaluer la position du point de réattachement) * DOMTOT = 'CHAN' 'QUAF' (mesh1 ET mesh2) ; $DOMTOT = 'MODE' DOMTOT 'NAVIER_STOKES' DISCR ; $INLET = 'MODE' p7p1 'NAVIER_STOKES' DISCR ; $OUTLET = 'MODE' (p4p8 'ET' p8p5) 'NAVIER_STOKES' DISCR ; $BOTTOM = 'MODE' p3p4 'NAVIER_STOKES' DISCR ; $WALL = 'MODE' (p1p2 ET p2p3 ET p3p4 ET p5p6 ET p6p7) 'NAVIER_STOKES' DISCR ; * * Elimination ad hoc * (En 2D, il faut éliminer les points centres des modèles surfaciques * avec les points faces des modèles volumiques à cause des MACROs) FDOMTOT = 'DOMA' $DOMTOT 'FACE' ; CINLET = 'DOMA' $INLET 'CENTRE' ; COUTLET = 'DOMA' $OUTLET 'CENTRE' ; CBOTTOM = 'DOMA' $BOTTOM 'CENTRE' ; CWALL = 'DOMA' $WALL 'CENTRE' ; 'ELIM' (FDOMTOT 'ET' CINLET 'ET' COUTLET 'ET' CWALL 'ET' CBOTTOM) EPS0 ; * * On écrase les anciens maillages afin d'éviter toute ambiguitée DOMTOT = 'DOMA' $DOMTOT 'MAILLAGE' ; INLET = 'DOMA' $INLET 'MAILLAGE' ; OUTLET = 'DOMA' $OUTLET 'MAILLAGE' ; BOTTOM = 'DOMA' $BOTTOM 'MAILLAGE' ; WALL = 'DOMA' $WALL 'MAILLAGE' ; * * Maillage pour d'éventuelles conditions aux limites en pression 'SI' ('EGA' kpres 'MSOMMET') ; OUTLETP = 'DOMA' $outlet KPRES ; 'SINON' ; OUTLETP = 'DOMA' $domtot KPRES ; 'FINSI' ; * *=========================== * Description des équations *=========================== * * Grandeurs adimensionnées Umax = 1.0 ; Re = 150. ; * * Profil de vitesse parabolique à l'entrée YINLET = 'COOR' 2 INLET ; YMAX = 'MAXI' YINLET ; YMIN = 'MINI' YINLET ; UIN = (YINLET '-' YMAX) '*' (YINLET '-' YMIN) ; UIN = UIN '*' (-4.0*Umax/((YMAX-YMIN)*(YMAX-YMIN))) ; UIN = 'NOMC' 'UX' UIN 'NATU' 'DISCRET' ; VIN = 'KCHT' $INLET 'SCAL' 'SOMMET' 'COMP' 'UY' 0. ; * * Description du système en vitesse-pression RV = 'EQEX' $DOMTOT 'ITMA' 5000 'ALFA' 1. 'FIDT' 20000 'OPTI' 'EF' 'IMPL' KSUPG KPRES 'ZONE' $DOMTOT 'OPER' residu FREQ0 EPS0 GRAPH 'ZONE' $DOMTOT 'OPER' 'NS' 1. 'UN' (1./Re) 'INCO' 'UN' * Conditions aux limites par défaut * 'ZONE' $OUTLET 'OPER' 'TOIM' (0. 0.) 'INCO' 'UN' 'OPTI' KMASS KSUPG 'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'UN' DT 'UN' (1./Re) 'INCO' 'UN' ; 'SI' ('EGA' kpres 'MSOMMET') ; RVP = 'EQEX' 'OPTI' 'EF' KPRES 'IMPL' 'ZONE' $DOMTOT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES' 'CLIM' 'PRES' 'TIMP' OUTLETP 0. ; 'SINON' ; RVP = 'EQEX' 'OPTI' 'EF' KPRES 'IMPL' 'ZONE' $DOMTOT 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES' * En milieu ouvert, pression discontinue il est inutile * d'imposer la pression en un point * 'CLIM' 'PRES' 'TIMP' (OUTLETP 'ELEM' 1) 0. ; 'FINSI' ; * * Description des conditions aux limites RV = 'EQEX' RV 'CLIM' 'UN' 'UIMP' WALL 0. 'UN' 'VIMP' WALL 0. 'UN' 'UIMP' INLET UIN 'UN' 'VIMP' INLET VIN ; * * Déclaration des inconnues et initialisations (table INCO) RV . 'INCO' = 'TABLE' 'INCO' ; RV . 'INCO' . 'UN' = 'KCHT' $DOMTOT 'VECT' 'SOMMET' (0. 0.) ; RV . 'INCO' . 'PRES' = 'KCHT' $DOMTOT 'SCAL' KPRES 0. ; * * Champs supplémentaires pour la procédure residu RV . 'INCO' . 'IT' = 'PROG' ; RV . 'INCO' . 'TI' = 'PROG' ; RV . 'INCO' . 'ER' = 'PROG' ; RV . 'INCO' . 'POSI' = 'PROG' ; RV . 'INCO' . 'UN2' = 'KCHT' $DOMTOT 'VECT' 'SOMMET' (1.D-3 1.D-3) ; * * Algorithme de résolution : méthode de projection RV . 'PROJ' = RVP ; * EXEC RV ; * *================= * Post-traitement *================= * * Localisation du point de réattachement * Tracés de la vitesse et de la pression * CNT1 = WALL ; NLI0 = 14 ; * * Point de réattachement (point où dUx/dy=0) V1 = attac ; 'SAUT' 10 'LIGNE' ; 'MESS' 'ABSCISSE DU POINT DE REATTACHEMENT :' V1 ; 'SI' GRAPH ; EV1 = 'EVOL' 'MANU' RV . 'INCO' . 'TI' RV . 'INCO' . 'POSI' ; 'DESS' EV1 'MIMA' 'GRIL' 'TITR' 'Localisation reattachement = f(t)' ; * * Vitesse un = RV . 'INCO' . 'UN' ; vun = 'VECT' UN 0.1 'UX' 'UY' 'JAUNE' ; trace vun DOMTOT CNT1 'TITR' 'Vitesse' ; * * Pression pe = 'EXCO' (RV . 'INCO' . 'PRESSION') 'PRES' 'SCAL' ; 'SI' ('EGA' kpres 'MSOMMET') ; mp = 'DOMA' $DOMTOT 'MMAIL' ; trace pe mp ('CONTOUR' mp) ; 'SINO' ; pn = 'ELNO' $DOMTOT ('KCHT' $DOMTOT 'SCAL' kpres 0. pe) kpres; trace pn domtot CNT1 NLI0 'TITR' 'Pression' ; 'FINSI' ; 'FINSI' ; * *======================== * Test de non régression *======================== * 'SI' COMPLET ; TEST = ('ABS' (V1 '-' 6.1639)) < 0.0005 ; 'SINON' ; TEST = ('ABS' (V1 '-' 5.5586)) < 0.0005 ; 'FINSI' ; 'MESSAGE' 'Element Vitesse : ' DISCR ; 'MESSAGE' 'Element Pression : ' KPRES ; 'MESSAGE' 'Décentrement : ' KSUPG ; 'MESSAGE' 'Matrice masse : ' KMASS ; 'SI' TEST ; 'ERREUR' 0 ; 'SINON' ; 'ERREUR' 5 ; 'FINSI' ; * 'FIN' ;