* fichier : cyltest.dgibi ************************************************************************ * Section : Fluides Non Stationnaire ************************************************************************ ************************************************************************ * * * Ecoulement autour d'un cylindre circulaire * * Résolution des équations de Navier Stokes laminaires instationnaires * * * ************************************************************************ *'OPTION' 'TRAC' 'PSC'; 'OPTION' 'TRAC' 'X' ; 'OPTION' isov suli; GRAPH = VRAI ; GRAPH = FAUX ; COMPLET = FAUX; *COMPLET = VRAI; 'SI' (COMPLET); * Temps final de calcul : tfin = 2.; maxit = 2000; DT0 = 0.05 ; * Définitions des nombres d'éléments: nx1 = 10; nx2 = 12; nx3 = 25; * Définitions des densites: dbeg = 0.1; dend = 0.5; 'SINON'; * Temps final de calcul : tfin = 2.; maxit = 200; DT0 = 0.5 ; * Définitions des nombres d'éléments: nx1 = 6; nx1 = 3; nx2 = 4; nx3 = 14; nx3 = 7 ; * Définitions des densites: dbeg = 0.4; dend = 0.5; dend = 1. ; 'FINSI'; *===================== * Domaine de calcul 2D: *===================== 'OPTION' 'DIME' 2 'ELEM' 'QUA8'; DISP = 'CENTRE'; ***************************************** * * * Procedure calcul * * * ***************************************** 'DEBP' CALCUL; * Calcul de la convergence. * Evolution de la vitesse en un point quelconque. 'ARGUMENT' RVX*'TABLE'; dd = RV.PASDETPS.'NUPASDT'; nn = dd/2; nf=2 ; nfg=50 ; nn = dd/nf; 'SI' ((dd-(nf*nn)) 'EGA' 0); res = 1.E-6; * Champs de vitesse aux instants n et n-1: un = RV.INCO.'UN'; unm1 = RV.INCO.'UNM1'; * Extraction des composantes de la vitesse aux pas de temps n et n-1: unx = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' un); unm1x = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UX' unm1); uny = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' un); unm1y = 'KCHT' $domtot 'SCAL' SOMMET ('EXCO' 'UY' unm1); * Sauvegarde de la vitesse en un point fixe du domaine à l'instant n: z = 'REDU' unx RV.INCO.'PoiObs'; zz = 'PROG' ('EXTR' z 'SCAL' RV.INCO.'PoiObs'); RV.INCO.'Uvit'= RV.INCO.'Uvit' 'ET' zz; * Calcul de la convergence: errx = 'KOPS' unx '-' unm1x; erry = 'KOPS' uny '-' unm1y; elix = 'MAXIMUM' errx 'ABS'; eliy = 'MAXIMUM' erry 'ABS'; elix = 'LOG'(elix + 1.0E-10)/('LOG' 10.); eliy = 'LOG'(eliy + 1.0E-10)/('LOG' 10.); erx = 'PROG' elix; ery = 'PROG' eliy; RV.INCO.'ERX' = RV.INCO.'ERX' et erx; RV.INCO.'ERY' = RV.INCO.'ERY' et ery; it = 'PROG' RV.PASDETPS.'NUPASDT'; temps = 'PROG' RV.PASDETPS.'TPS'; RV.INCO.'IT' = RV.INCO.'IT' et it; RV.INCO.'tps' = RV.INCO.'tps' et temps; RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' (RV.INCO.'UN'); PRES = RV.INCO.'PRESSION'; PRESNM1 = RV.INCO.'PRESNM1'; 'SI' ((dd-(nfg*nn)) 'EGA' 0); * Tracé du champ de pression PN = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE' ('NOMC' 'SCAL' (RV.INCO.'PRESSION'))); 'TRACER' PN domtot 'TITRE' 'Champ de pression'; * Tracé du champ de pression lambda PT = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE' ('NOMC' 'SCAL' (RV.INCO.'PTILDE'))); * 'TRACER' PT domtot 'TITRE' 'Champ de pression lambda'; * Tracé de la pression sur le cylindre presc = 'REDU' PN £cercle; precer = 'EVOL' (rouge) 'CHPO' presc cercle ; * 'DESSIN' precer 'TITRE' 'Valeurs de pression sur le cercle'; * Tracé du champ de vitesse VX = 'EXCO' 'UX' (RV.INCO.'UN'); VY = 'EXCO' 'UY' (RV.INCO.'UN'); modu = ((VX**2) '+' (VY**2))**(0.5); umax = 'MAXIMUM' modu; CHUN = 'VECTEUR' (RV.INCO.'UN') (1. '/' (1. * umax)) UX UY BLEU; 'TRACER' CHUN £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de vitesse'; * Tracé de la vitesse au point d'observation evv = 'EVOL' 'MANU' 'Temps' (RV.INCO.'tps') 'Valeur de U' (RV.INCO.'Uvit'); * 'DESSIN' evv 'TITRE' 'Evolution de la vitesse au point (0.8 0.5)'; * Tracé de la convergence evconv = 'EVOL' 'MANU' 'Nombre d itérations' (RV.INCO.'IT') 'log(err)' (RV.INCO.'ERY'); * 'DESSIN' evconv 'TITRE' 'Convergence sur la vitesse'; 'FINSI'; * Convergence en pression : norme l2 * mess 'Convergence en pression : norme l2'; 'SI' ('EGA' (RV.PASDETPS.'NUPASDT') 1); errpl2 = 1.; 'SINON'; DPDT = PRES - PRESNM1; DP2 = (nomc dpdt 'SCAL') * (nomc dpdt 'SCAL'); * DP2 = 'KOPS' DPDT 'PSCA' DPDT; DP2E = 'SOMT' DP2; errpl2 = DP2E '/' volt; errpl2 = errpl2**0.5; 'MESSAGE' 'convergence sur la pression' errpl2; 'FINSI'; errp = 'PROG' errpl2; RV.INCO.'errp'= RV.INCO.'errp' 'ET' errp; RV.INCO.'PRESNM1' = 'COPIER' RV.INCO.'PRESSION'; 'SI' (errpl2 < res); 'MESSAGE' 'La pression est résolue à 10E-6'; 'QUITTER' CALCUL; 'FINSI'; 'FINSI'; as2 ama1 = 'KOPS' 'MATRIK'; 'RESPRO' as2 ama1; 'FINP'; *============================ * Construction du maillage: *============================ * Construction des points: p0 = 0. 0.; p1 = -6. 0.; p2 = -0.5 0.; p3 = 0. 0.5; p4 = 0.5 0.; p5 = 15. 0.; p6 = 15. 6.; p7 = 3. 6.; p8 = -3. 6.; p9 = -6. 6.; * construction des points de la partie supérieure: p1p2 = p1 'DROIT' ((-1)*nx1) dini 1. dfin 0.1 p2; p2p3 = 'CERCLE' (nx2/2) p2 p0 p3; p3p4 = 'CERCLE' (nx2/2) p3 p0 p4; p4p5 = p4 'DROIT' ((-1)*nx3) dini 0.1 dfin 1.2 p5; p5p6 = p5 'DROIT' dini dbeg dfin dend p6; p6p7 = p6 'DROIT' ((-1)*nx3) dini 1.2 dfin 0.1 p7; p7p8 = p7 'DROIT' nx2 p8; p8p9 = p8 'DROIT' ((-1)*nx1) dini 0.1 dfin 1. p9; p9p1 = p9 'DROIT' dini dend dfin dbeg p1; * Création du contour: peri1 = p1p2 et p2p3 et p3p4 et p4p5; peri2 = p6p7 et p7p8 et p8p9; s1 = 'DALL' peri1 p5p6 peri2 p9p1 'PLAN'; * construction des points de la partie inférieure: p10 = -6. -6.; p11 = -3. -6.; p12 = 3. -6.; p13 = 15. -6.; p14 = 0. -0.5; p1p10 = p1 'DROIT' dini dbeg dfin dend p10; p10p11 = p10 'DROIT' ((-1)*nx1) dini 1. dfin 0.1 p11; p11p12 = p11 'DROIT' nx2 p12; p12p13 = p12 'DROIT' ((-1)*nx3) dini 0.1 dfin 1.2 p13; p13p5 = p13 'DROIT' dini dend dfin dbeg p5; p5p4 = 'INVERSE' p4p5; p4p14 = 'CERCLE' (nx2/2) p4 p0 p14; p14p2 = 'CERCLE' (nx2/2) p14 p0 p2; p2p1 = 'INVERSE' p1p2; * Création du contour: peri3 = p10p11 et p11p12 et p12p13; peri4 = p5p4 et p4p14 et p14p2 et p2p1; s2 = 'DALL' peri3 p13p5 peri4 p1p10 'PLAN'; cercle = (p2p3 et p3p4 et p4p14 et p14p2); 'ELIMINER' 1.E-3 s1 s2; * Création du domaine total: domtot = s1 et s2; DOM = 'CONTOUR' domtot; *============================================== * Changement des éléments du maillage en QUAF: *============================================== £domtot = 'CHANGER' domtot 'QUAF'; £cercle = 'CHANGER' cercle 'QUAF'; £entree = 'CHANGER' (p9p1 et p1p10) 'QUAF'; £sortie = 'CHANGER' (p5p6 et p13p5) 'QUAF'; £peri2 = 'CHANGER' peri2 'QUAF'; £peri3 = 'CHANGER' peri3 'QUAF'; £gau = 'CHANGER' (p9p1 et p1p10) 'QUAF'; £droi = 'CHANGER' (p5p6 et p13p5) 'QUAF'; £limite = 'CHANGER' (peri2 et p9p1 et p1p10 et peri3) 'QUAF'; 'ELIM' 1.E-3 (£domtot et £cercle et £entree); NN = 'NBEL' domtot; 'MESSAGE' 'Nombre d éléments du maillage' NN; *======================================= * Formulation du domaine Navier Stokes: *======================================= $domtot = 'MODE' £domtot 'NAVIER_STOKES' 'MACRO'; $entree = 'MODE' £entree 'NAVIER_STOKES' 'MACRO'; $sortie = 'MODE' £sortie 'NAVIER_STOKES' 'MACRO'; mdomtot = 'DOMA' $domtot maillage; modelt = 'MODELE' domtot 'THERMIQUE'; vol = 'DOMA' $domtot 'VOLUME'; volt = 'SOMT' vol; *==================================== * Création des tables de résolution: *==================================== rey = 300.; nu = 1./rey; arelax = 1.; RV = 'EQEX' 'OMEGA' 1. 'NITER' 1 'ITMA' maxit 'FIDT' 1 'ZONE' $domtot 'OPER' CALCUL; RV = 'EQEX' RV 'OPTI' 'EFM1' 'CENTREE' 'IMPL' * 'OPTI' 'EF' 'CENTREE' 'ZONE' $domtot 'OPER' 'DFDT' 1. 'UN' 'DT' 'INCO' 'UN'; RV = 'EQEX' RV 'OPTI' 'EF' 'CENTREE' 'IMPL' 'ZONE' $domtot 'OPER' 'NS' 1. 'UN' 'INV_RE' 'INCO' 'UN' * 'ZONE' $domtot 'OPER' 'KONV' 1. 'UN' 'INV_RE' 'INCO' 'UN' * 'ZONE' $domtot 'OPER' 'LAPN' 'INV_RE' 'INCO' 'UN' ; * Implantation des conditions aux limites: RV = 'EQEX' RV 'CLIM' 'UN' 'UIMP' £cercle 0. 'UN' 'VIMP' £cercle 0. 'UN' 'UIMP' £entree 1. 'UN' 'VIMP' £entree 0. 'UN' 'VIMP' £peri2 0. 'UN' 'VIMP' £peri3 0. ; *RV = 'EQEX' RV * 'CLIM' 'UN' 'UIMP' £cercle 0. 'UN' 'VIMP' £cercle 0. * 'UN' 'UIMP' £limite 1. 'UN' 'VIMP' £limite 0.; RV.'METHINV'.TYPINV=3 ; *RV.'METHINV'.TYPINV=1 ; RV.'METHINV'.IMPINV=0 ; RV.'METHINV'.NITMAX=300 ; RV.'METHINV'.PRECOND=3 ; RV.'METHINV'.RESID =1.e-8 ; RV. 'METHINV' . 'FCPRECT'=100 ; RV. 'METHINV' . 'FCPRECI'=100 ; RVP = 'EQEX' 'OPTI' 'EF' 'IMPL' INCOD DISP 'ZONE' $domtot 'OPER' 'KBBT' -1. 1. 'INCO' 'UN' 'PRES' ; RV.'PROJ'= RVP ; RVP.'METHINV'.TYPINV =1 ; RVP.'METHINV'.IMPINV =0 ; RVP.'METHINV'.NITMAX =150; RVP.'METHINV'.PRECOND =3 ; RVP.'METHINV'.RESID =1.e-12; RVP.'METHINV'.'FCPRECT'=100 ; RVP.'METHINV'.'FCPRECI'=100 ; * Création de la table des inconnues et initialisation: RV.INCO = 'TABLE' INCO; RV.INCO.'UN' = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3); RV.INCO.'UNM1' = 'KCHT' $domtot 'VECT' 'SOMMET' (1.E-3 1.E-3); RV.INCO.'DT' = DT0 ; RV.INCO.'INV_RE'= nu; RV.INCO.'EVOLP' = 'PROG' 0.; RV.INCO.'PRES' = 'KCHT' $domtot 'SCAL' DISP 0.; RV.INCO.'PRESNM1' = 'KCHT' $domtot 'SCAL' DISP 'COMP' 'PRES' 0.; RV.INCO.'PTILDE' = 'KCHT' $domtot 'SCAL' DISP 0.; RVP.INCO = RV.INCO; * Détermination d'un point d'observation: obs = mdomtot 'POINT' 'PROC' (0.8 0.5); RV.INCO.'PoiObs'= obs; * Initialisation des inconnues complémentaires: RV.INCO.'Uvit' = 'PROG' 0.; RV.INCO.'ERX' = 'PROG' 0.; RV.INCO.'ERY' = 'PROG' 0.; RV.INCO.'errp' = 'PROG' 0.; RV.INCO.'IT' = 'PROG' 0; RV.INCO.'tps' = 'PROG' 0.; RV.INCO.'entree'= £entree; RV.INCO.'cercle'= £cercle; RV.'TYPROJ'='VPI2'; EXEC RV; * ------------ Controle calcul ---------------------------- Si (NON COMPLET) ; PN = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE' ('NOMC' 'SCAL' (RV.INCO.'PRESSION'))); presc = 'REDU' PN £cercle; precer = 'EVOL' (rouge) 'CHPO' presc cercle ; tefp= extr precer 'ORDO' ; list tefp ; refp = PROG .38960 .15815 -.34055 -.47886 -.42818 -.20337 -6.11799E-03 -.12634 -.26610 -.47590 -.58950 -1.0184 -1.3341 -1.1573 -.73313 -4.62908E-02 .38960; errpc = abs ( tefp - refp ) maxi ; mess ' errpc=' errpc ; Si (errpc > 0.005 ); erreur 5 ; Finsi ; Finsi ; * ------------ Post traitement ---------------------------- Si GRAPH ; VX = 'EXCO' 'UX' (RV.INCO.'UN'); VY = 'EXCO' 'UY' (RV.INCO.'UN'); modu = ((VX**2) '+' (VY**2))**(0.5); umax = 'MAXIMUM' modu; CHUN = 'VECTEUR' (RV.INCO.'UN') (1. '/' (1. * umax)) UX UY BLEU; 'TRACER' CHUN £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de vitesse'; * Tracé du champ de pression 'TRACER' PN £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de pression'; * Tracé du champ de pression lambda PT = 'ELNO' $domtot ('KCHT' $domtot 'SCAL' 'CENTRE' ('NOMC' 'SCAL' (RV.INCO.'PTILDE'))); 'TRACER' PT £domtot ('CONTOUR' £domtot) 'TITRE' 'Champ de pression lambda'; * Tracé de la pression sur le cylindre 'DESSIN' precer 'TITRE' 'Valeurs de pression sur le cercle'; * Tracé de la vitesse au point d'observation evv = 'EVOL' 'MANU' 't' (RV.INCO.'tps') 'Valeur de U' (RV.INCO.'Uvit'); 'DESSIN' evv 'TITRE' 'Evolution de la vitesse au point (0.8 0.5)'; * Tracé de la convergence evconv = 'EVOL' 'MANU' 'Iteration' (RV.INCO.'IT') 'log(err)' (RV.INCO.'ERY'); 'DESSIN' evconv 'TITRE' 'Convergence Vitesse '; evconvp = 'EVOL' 'MANU' 'Iteration' (RV.INCO.'IT') 'errp' (RV.INCO.'errp'); 'DESSIN' evconvp 'TITRE' 'Convergence pression'; * Calcul des lignes de courant: y = 'COOR' 2 (RV.INCO.'entree'); psiim = 'KOPS' 1. '*' y; psii = 'NOMC' 'psi' psiim 'NATURE' 'DISCRET'; vitesse = RV.INCO.'UN'; tourb = 'KOPS' vitesse 'ROT' $domtot; RK = 'EQEX' $domtot 'OPTI' 'EF' 'SEMI' 0.5 *rk = 'EQEX' $domtot 'OPTI' 'EF' 'IMPL' 'ZONE' $domtot 'OPER' 'LAPN' 1. 'INCO' 'PSI' 'ZONE' $domtot 'OPER' 'FIMP' tourb 'INCO' 'PSI' 'CLIM' 'PSI' 'TIMP' (RV.INCO.'cercle') 0. 'CLIM' 'PSI' 'TIMP' (RV.INCO.'entree') psii; RK.'INCO' = table 'INCO'; RK.'INCO'.'PSI' = 'KCHT' $domtot 'SCAL' 'SOMMET' 0.; EXEC RK; * Traçage des lignes de courant: psi =RK.'INCO'.'PSI'; *'OPTION' trac ps; 'OPTION' isov ligne; 'LISTE' ('MAXIMUM' PSI); 'LISTE' ('MINIMUM' PSI); 'TRACER' psi £domtot ('CONTOUR' £domtot) ('TITRE' 'Lignes de courant, Re='rey 'à t='RV.PASDETPS.'TPS'); Finsi ; * 'FIN' du jeu de données. fin;