* fichier : canalBu.dgibi
* Ce cas teste le modèle longueur de melange de Buleev
* sur l'écoulement turbulent dans un canal plan.
opti trace 'PSC';
mess ' Configuration Courte ' ; COMPLET = FAUX ;
mess ' Configuration SANS TRACE' ; GRAPH = FAUX ;
*mess ' Configuration Complete ' ; COMPLET = VRAI ;
*mess ' Configuration POST ' ; GRAPH = VRAI ;
*-----------------------------------------------------------------------
opti dime 2 elem TRI6;
opti dime 2 elem QUA8;
*-----------------------------------------------------------------------
* |Triangles | Quadrangles
* ---------------------------
* LINE | | OK
* LINB | |
* MACRO| ? | OK
* QUAF | | OK
*-----------------------------------------------------------------------
DISCR=QUAF;
DISCR=MACRO;
*DISCR=LINE;
KPRES=CENTREP1;
*KPRES=MSOMMET ;
*KPRES=CENTRE ;
KSUPG=SUPG;
KSUPG=CENTREE;
KSUPG=SUPG ;
epsi = 1.e-5;
******=================================================**********
******=================================================**********
******=================================================**********
******=================================================**********
******=================================================**********
Si COMPLET;
NBIT =250;
DT = 0.05;
H=1. ;
L=30.*H;
dcl=H/120;
dm=H/20;
dcl=H/160;
dm=H/25;
dcl=H/1200;
dm=H/40;
dcx=H/15.;
dmx=H/10.;
Sinon;
NBIT =7;
*NBIT =100;
DT = 0.1;
H=1. ;
L=10.*H;
dcl=H/20;
dm=H/5;
dcx=H/10.;
dmx=H/5.;
Finsi;
*................................................................
DEBPROC CALCUL;
ARGU RVX*'TABLE';
RV = RVX.'EQEX' ;
un = rv.'INCO'.'UN' ;
* Prt = PRODT UN $mt (GB T) ;
* Prt = PRODT UN $mt ;
Prt = (exco (kops 'GRADS' (exco un 'UX') $mt) 'UY') abs;
rv.'INCO'.'PRODT'=Prt;
Si ( EGA LMEL 1) ;
mut = Ro * lm1 * lm1 * (Prt**0.5);
Finsi;
Si ( EGA LMEL 2 ) ;
mut = Ro * lm2 * lm2 * Prt ;
Finsi;
Si ( EGA LMEL 3 ) ;
cm2 = 7.1*7.1;
uaxe= exco un 'UX';
rett= ro*uaxe*lm3*(inve rv.'INCO'.'MU');
rett=rett - cm2;
irett= masq rett 'SUPERIEUR' 0.;
mut = rett * irett * (rv.inco.'MU' * (1./cm2));
Finsi;
rv.'INCO'.'MUF'= mut + rv.'INCO'.'MU';
as2 ama1 = 'KOPS' 'MATRIK' ;
* mess ' FIN CALCUL' ;
'FINPROC' as2 ama1 ;
******=================================================**********
******=================================================**********
******=================================================**********
******=================================================**********
******=================================================**********
P0 = 0. 0. ;
P1 = L 0. ;
P2 = L H ;
P3 = 0. H ;
PS3 = (P0 PLUS P3)*0.5;
entree= P0 d dini dcl dfin dm ps3 d dini dm dfin dcl P3;
*entree= P0 d 12 P3;
sortie=(entree plus (L 0.)) inve;
par1 = p0 d dini dcx dfin dmx p1 ;
*par1 = p0 d dini dcl dfin dcl p1 ;
par2 = (par1 plus (0. H)) inve;
pc= 1. 1.;
ang= 43.;
ang= 60.;
ang= 0. ;
entree=entree TOUR ang pC;
sortie=sortie TOUR ang pC;
par1 =par1 TOUR ang pC;
par2 =par2 TOUR ang pC;
elim (entree et par1 et par2 et sortie) epsi;
mt = daller entree par2 sortie par1;
v1=1. 0. ;
v2=0. -1.;
v3=-1. 0.;
v4=0. 1. ;
vt1= tour v1 ang p0;
vt2= tour v2 ang p0;
vt3= tour v3 ang p0;
vt4= tour v4 ang p0;
list vt1;
list vt2;
list vt3;
list vt4;
Mmt = chan mt QUAF ;
Mentree= chan entree QUAF ;
Msortie= chan sortie QUAF ;
Mpar1 = chan par1 QUAF ;
Mpar2 = chan par2 QUAF ;
elim (Mmt et Mentree et Msortie et Mpar1 et Mpar2) epsi ;
$mt = model Mmt 'NAVIER_STOKES' DISCR;
$entree = model Mentree 'NAVIER_STOKES' DISCR;
$sortie = model Msortie 'NAVIER_STOKES' DISCR;
$par1 = model Mpar1 'NAVIER_STOKES' DISCR;
$par2 = model Mpar2 'NAVIER_STOKES' DISCR;
doma $mt 'IMPR';
doma $entree 'IMPR';
mt = doma $mt maillage;
entree = doma $entree maillage;
sortie = doma $sortie maillage;
par1 = doma $par1 maillage;
par2 = doma $par2 maillage;
axe = par1 plus (0. (H/2.));
elim (axe et mt) epsi;
entre1= chan entree poi1;
entre1= elem entre1 (lect 2 pas 1 ((nbel entre1) - 1 ));
*****************************************************************
*** Création des longueurs de mélange
* LMEL=1
ck=0.41;
lm1=0.41 *(coor 2 mt);
ilm1= lm1 masq 'INFERIEUR' (0.4 *H/2.) ;
lm2=0.41*(H - (coor 2 mt));
ilm2= lm2 masq 'INFERIEUR' (0.4 *H/2.) ;
lma=0.41 *(coor 2 axe);
lm1= (lm1*ilm1)+(lm2*ilm2) + lma;
Ef=H/2.;
* LMEL=2
yr = (coor 2 mt)*(1./Ef/2.);
lm2 = 2.*Ef*(0.14 - (0.08*(1. - yr)**2.) - (0.06*(1. - yr)**4.));
* LMEL=3
dparoi=(coor 2 mt);
dparoi=dparoi*(H - dparoi);
dparoi=dparoi+1.e-10;
zd = dparoi ;
*lm3 = H *( (Zd / H) * (1. - (Zd/H)) ) ;
lm3 = Zd ;
evlm1 = evol chpo lm1 sortie;
evlm1 = evlm1 coul rouge;
evlm2 = evol chpo lm2 sortie;
evlm2 = evlm2 coul jaune;
evlm3 = evol chpo lm3 sortie;
evlm3 = evlm3 coul verte;
evdpar=evol chpo dparoi sortie;
evdpar= evdpar coul bleus;
TAB1 = TABLE ;
tab1.1 = 'MOT' 'TIRR';
*tab1.2 = 'MOT' 'TIRC';
tab1.2 = 'MOT' 'TIRM';
tab1.3 = 'MOT' '' ;
tab1.'TITRE' = TABLE ;
tab1.'TITRE'. 1 = 'MOT' 'Karman' ;
*tab1.'TITRE'. 2 = 'MOT' 'lm2' ;
tab1.'TITRE'. 2 = 'MOT' 'Buleev' ;
tab1.'TITRE'. 3 = 'MOT' 'd paroi' ;
Si GRAPH ;
trace mt TITR ' Maillage ';
*dess (evlm1 et evlm2 et evlm3 et evdpar)
dess (evlm1 et evlm3 et evdpar)
TITR ' Longueurs de mélange' 'GRIL' lege TAB1;
Finsi ;
*****************************************************************
ye= coor 2 entre1;
U0 = 1. ;
*U0 = 5. ;
Ufr = U0 ;
*Ufr = ye*(H - ye)*U0*1.5*4. ;
uux=(cos ang) * U0;
uuy=(sin ang) * U0;
MU = 3.14E-5;
Ro = 1. ;
YP = 1.e-2 ;
*Modèle de turbulence : Buleev
LMEL=3;
Rey= Ro * Ufr *H * 2. / MU;
mess ' Reynolds = ' Rey;
RV= EQEX 'OMEGA' 1. 'NITER' 1 ITMA NBIT 'FIDT' 1
'OPTI' 'EF' 'IMPL' KSUPG KPRES 'MUVARI'
ZONE $mt 'OPER' 'CALCUL'
ZONE $mt 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN'
'OPTI' 'EFM1' 'IMPL' 'CENTREE'
ZONE $mt 'OPER' DFDT 1. 'UNM' DT INCO 'UN';
;
rv.'ALGO_KEPSILON'=mots 'KLbr';
RV= EQEX RV CLIM
* UN UIMP entre1 Ufr UN VIMP entree 1.e-4
UN UIMP entre1 uux UN VIMP entree (uuy + 1.e-4)
UN UIMP (par1 et par2) 0. UN VIMP (par1 et par2) 0.
;
rv.'METHINV'.TYPINV=3 ;
rv.'METHINV'.IMPINV=0 ;
rv.'METHINV'.NITMAX=5000;
rv.'METHINV'.PRECOND=3 ;
*rv.'METHINV'.PRECOND=2 ;
rv.'METHINV'.ILUTLFIL=4 ;
rv.'METHINV'.ILUTDTOL=0.;
rv.'METHINV'.RESID =1.e-8;
rv.'METHINV' . 'FCPRECT'=1 ;
rv.'METHINV' . 'FCPRECI'=1 ;
RVP= EQEX
'OPTI' 'EF' KPRES
ZONE $mt 'OPER' KBBT (-1.) 'INCO' 'UN' 'PRES'
;
* rvp.'METHINV'.TYPINV=TYPINVPR ;
rvp.'METHINV'.TYPINV=3 ;
rvp.'METHINV'.IMPINV=0 ;
rvp.'METHINV'.NITMAX=2500;
rvp.'METHINV'.PRECOND=3 ;
rvp.'METHINV'.ILUTLFIL=10;
rvp.'METHINV'.ILUTDTOL=0.;
rvp.'METHINV'.RESID =1.e-10 ;
rvp.'METHINV' . 'FCPRECT'=50000 ;
rvp.'METHINV' . 'FCPRECI'=50000 ;
RV.'PROJ'= RVP ;
rv.'INCO'= table inco ;
rv.'INCO'.'UN' = kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.'INCO'.'MUF' = kcht $mt scal sommet MU ;
rv.'INCO'.'Echl'= kcht $mt scal sommet dparoi ;
rv.'INCO'.'dparoi'= kcht $mt scal sommet dparoi ;
rv.'INCO'.'EN' = kcht $mt scal sommet 1.e-5 ;
rv.'INCO'.'UNM' = kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.'INCO'.'UNMM'= kcht $mt vect sommet (1.e-5 1.e-5) ;
rv.'INCO'.'PRES'= kcht $mt scal kpres 0.;
rv.'INCO'.'MU' = kcht $mt scal sommet MU;
rv.'LTPS' = PROG ;
rv.'Tps' = 1. ;
rv.'DT' = DT;
rv.'NUPAT' =1;
mess ' EXECUTION '; exec rv ;
*Si VRAI ; opti rest 'KBU'; rest ; Finsi;
un = rv.inco.'UN' ;
qs=dbit un $sortie IMPR;
um=(1./H)*qs;
unn=(1./um)*un;
*-----------------------------------------------------------------------
*--- Edition d'un profil de vitesse asymptotique de Référence ----------
*--- Modèle de Buleev
Debproc ANABU H*'FLOTTANT' Rey*'FLOTTANT' UM*'FLOTTANT' ;
H2=H*H;
p0=0. 0. ; p1=0. H. ; Pm=0. (H/2.); dpar=H/Rey*0.05; dctr=H*0.05;
de= p0 d dini dpar dfin dctr pm d dini dctr dfin dpar p1 ;
z= coor 2 de;
Dh=2.*H;
Rho=1.;
mu=um/Rey*Rho*Dh;
m=7.1;
nu=mu/rho;
uet=um;
REPETER CCP 15 ;
fp=(uet/Um)**2.*8. ;
F=fp/2/Dh*rho*(um*um);
xi = z*(1./H);
xiet = (1. - ((1. - (nu*m/h/Uet*4))**0.5))/2.;
u1=H2*F/2./mu*xi*(1. - xi);
ax=xi*(1. - xi)/xiet/(1. - xiet) + 1.e-10;
ax=m*m*F*H/Rho* (log ax);
ax =(abs ax) + 1.e-10;
ax0=(h*h*F/mu/2*xiet*(1. - xiet))**2.;
u2=(ax0+ax)**0.5;
iua=masq xi 'INFERIEUR' xiet;
iub=masq xi 'SUPERIEUR' (1. - xiet);
iu1=iua+iub;
uref = (iu1 * u1) + ((1. - iu1)*u2) ;
evuref = evol chpo uref de ;
umn=somme evuref ;
umn=(extr 1 umn)/H; Reyn=rho*umn*Dh/mu;
w=0.9;
uet=(w*(uet*um/umn)) + ((1. - w)*uet);
dumn=abs ((umn - um)/um);
mess 'Iteration ' &ccp ' Umn=' umn ' Reyn=' Reyn ' Uet=' Uet
' dumn=' dumn;
*dess evuref ;
Si (dumn < 1.e-4) ; quitter CCP; Finsi ;
Fin CCP;
* Calcul des u+g et u+d
yplus= (extr evuref 'ABSC')*uet/nu;
evuplus=evuref*(1./uet);
yy=extr evuplus ORDO ;
lyz=lect (dime yy) pas -1 1;
yyz=extr yy lyz;
evuplud = (evol manu yplus yyz) coul verte;
evuplug = (evol manu yplus yy) coul rouge ;
Resproc evuplug evuplud evuref uet;
Finproc ;
evupd1 evupg1 evuref1 uet1 = anabu 1. 2000. 1. ;
evupd2 evupg2 evuref2 uet2 = anabu 1. 1.e4 1. ;
evupd3 evupg3 evuref3 uet3 = anabu 1. 1.e5 1. ;
evupd4 evupg4 evuref4 uet4 = anabu 1. 1.e6 1. ;
evupd5 evupg5 evuref5 uet5 = anabu 1. 1.e7 1. ;
*--- Profil de Reichardt -----------------------------------------------
Ka = 0.41;
yplusr= prog 0.1 pas 0.1 30 pas 1. 100. pas 50. 1.e3
pas 500. 1.e4 pas 1.e4 2.e5 ;
unp = prog (dime yplusr)*1.;
E=9.5;
uplusr=(1./Ka)*(Log (unp + (Ka*yplusr)));
uplusr=uplusr + (7.8* (unp - (exp (yplusr*(-1./11))))
- ( (yplusr*(-1./11)) * (exp (yplusr*(-1./3.))) ) );
evr1=evol manu yplusr uplusr;
evr1=evr1 coul verte;
*--- Profil de Reichardt -----------------------------------------------
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
Si GRAPH;
TAB1 = TABLE ;
tab1.1 = 'MOT' ' ';
tab1.2 = 'MOT' 'MARQ CROI REGU';
tab1.3 = 'MOT' 'MARQ CARR REGU';
tab1.4 = 'MOT' 'MARQ LOSA REGU';
tab1.5 = 'MOT' 'MARQ PLUS REGU';
tab1.6 = 'MOT' 'MARQ TRIA REGU';
tab1.'TITRE' = TABLE ;
tab1.'TITRE'. 1 = 'MOT' 'Re 2.e3' ;
tab1.'TITRE'. 2 = 'MOT' 'Re 1.e4' ;
tab1.'TITRE'. 3 = 'MOT' 'Re 1.e5' ;
tab1.'TITRE'. 4 = 'MOT' 'Re 1.e6' ;
tab1.'TITRE'. 5 = 'MOT' 'Re 1.e7' ;
tab1.'TITRE'. 6 = 'MOT' 'Reichardt' ;
tit1=chai 'Comparaison Profils de vitesse Buleev à différents Reynds';
tit2=chai tit1 ' et Reichardt';
dess (evupd1 et evupd2 et evupd3 et evupd4 et evupd5 et evr1)
LOGX XBOR 0.1 5.e5 'GRIL' 'CARR' lege TAB1 titr tit2;
dess (evuref1 et evuref2 et evuref3 et evuref4 et evuref5)
'GRIL' 'CARR' lege TAB1 titr tit1;
*-----------------------------------------------------------------------
unx =exco un 'UX';
dudy=exco (kops 'GRADS' unx $mt) 'UY';
uet=(mu*(abs (dudy + 1.e-5)))**0.5;
uet1=redu uet par1;
uet2=redu uet par2;
evuet1= evol chpo uet1 par1;
evuet2= evol chpo uet2 (inve par2);
dess (evuet1 et evuet2) TITR ' U* le long des deux parois'
tity 'm/s' titx 'm';
mess ' Liste de U*';
yevuet1=extr evuet1 'ORDO';
yevuet2=extr evuet2 'ORDO';
dl=(dime yevuet1) ;
ll=lect (dl - 10) pas 1 dl;
list (extr yevuet1 ll) ;
list (extr yevuet2 ll) ;
f1=8.*(((1./um)*uet1)**2.);
f2=8.*(((1./um)*uet2)**2.);
fiso = 0.202*(Rey**(-0.196));
evf1=evol chpo f1 par1;
evf1=evf1 coul rouge;
evf2=evol chpo f2 (inve par2);
evf2=evf2 coul verte;
xx=extr evf1 'ABSC';
yfiso= prog (dime xx)*fiso;
evfiso=evol manu xx yfiso ;
evfiso=evfiso coul bleue ;
TAB1 = TABLE ;
tab1.1 = 'MOT' 'TIRR';
tab1.2 = 'MOT' 'TIRC';
tab1.3 = 'MOT' 'TIRM';
tab1.'TITRE' = TABLE ;
tab1.'TITRE'. 1 = 'MOT' 'Corr' ;
tab1.'TITRE'. 2 = 'MOT' 'f paroi1' ;
tab1.'TITRE'. 3 = 'MOT' 'f paroi2' ;
dess (evfiso et evf2 et evf1)
TITR ' Coefficient de frottement' 'GRIL' lege TAB1;
yp1=norm ((poin 2 entre1) moins (poin 1 entre1));
yp2=norm ((poin entre1 final) moins (poin ((nbno entre1) - 1) entre1));
mess ' YP1=' YP1 ' YP2=' YP2;
yplus1= yp1*uet1*(rv.inco.'MU'**(-1.));
yplus2= yp2*uet2*(rv.inco.'MU'**(-1.));
evyplus1=evol chpo yplus1 par1;
evyplus2=evol chpo yplus2 (inve par2);
dess (evyplus1 et evyplus2) TITR ' Y+ le long des deux parois'
tity ' ' titx 'm';
sort2=elem sortie (lect 1 pas 1 ((nbel sortie)/2));
sort1= syme sort2 'DROIT' (point 1 axe) (axe poin final);
elim (sort1 et sort2 et sortie) epsi;
uet1=extr (nbno par1) (extr evuet1 'ORDO');
uet2=extr (nbno par1) (extr evuet2 'ORDO');
mess ' uet1 et uet2=' uet1 uet2;
ypls1= evol chpo ((1./mu) * uet1 * (coor 2 sort1)) (inve sort1);
ypls2= evol chpo ((1./mu) * uet2 * (H - (coor 2 sort2))) ( sort2);
upl1=(evol chpo unx (inve sort1))*(1./uet1);
upl2=(evol chpo unx sort2)*(1./uet2);
xypls1= ORDO (extr ypls1 'ORDO');
yypls1= ORDO (extr upl1 'ORDO');
xypls1=xypls1 enlev 1;
yypls1=yypls1 enlev 1;
evupl1=evol manu xypls1 yypls1;
xypls2=extr ypls2 'ORDO';
yypls2=extr upl2 'ORDO';
xypls2=xypls2 enlev 1;
yypls2=yypls2 enlev 1;
evupl2=evol manu xypls2 yypls2 ;
evupl1=evupl1 coul rouge;
evupl2=evupl2 coul verte;
TAB1 = TABLE ;
tab1.1 = 'MOT' ' ';
tab1.2 = 'MOT' 'TIRC';
tab1.3 = 'MOT' 'TIRM';
tab1.'TITRE' = TABLE ;
tab1.'TITRE'. 1 = 'MOT' 'Reichardt' ;
tab1.'TITRE'. 2 = 'MOT' 'U+ paroi1' ;
tab1.'TITRE'. 3 = 'MOT' 'U+ paroi2' ;
dess (evr1 et evupl2 et evupl1) 'LOGX' 'XBOR' 1. 500.
TITR ' U+ en sortie' 'GRIL' lege TAB1
tity ' ' titx 'm';
nutsnu=(rv.inco.'MUF')*(1./MU);
trace nutsnu mt (cont mt) 'TITR' ' Nu effectif / nu ';
ung = vect un 0.5 'UX' 'UY' jaune;
trace ung mt 'TITR' ' Vitesses ';
pn = elno $mt rv.inco.PRESSION kpres;
trace pn mt (cont mt) TITR 'Pression';
evpaxe= evol chpo pn axe;
dess evpaxe TITR ' Pression sur l axe';
evax= evol chpo (exco un 'UX') axe ;
dess evax TITR 'U sur l axe'
tity 'm/s' titx 'm';
Finsi ;
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
*profils de référence
evupdr evupgr evusr uetr = anabu 1. Rey 1. ;
mess '***************************************************************';
mess '***************************************************************';
mess ' UET1 = ' UET1 ' UET2 = ' UET2 ' UETR = ' UETR ;
mess '***************************************************************';
mess '***************************************************************';
evus= evol chpo (exco unn 'UX') sortie;
usix = extr evusr 'ABSC';
usi = ipol usix (extr evus 'ABSC') (extr evus 'ORDO');
evusi = evol manu usix usi ;
evdu = (evusi - evusr)*(evusi - evusr);
Si GRAPH ;
tab1 = table;
tab1.1 = 'MOT' 'TIRC';
tab1.2 = 'MOT' ' ';
tab1.'TITRE' = TABLE ;
tab1.'TITRE'. 1 = 'MOT' 'Buleev' ;
tab1.'TITRE'. 2 = 'MOT' 'Buleev analytique' ;
evusr=evusr coul rouge;
dess (evusr et evus) TITR 'U en sortie' 'GRIL' LEGE TAB1
tity 'm/s' titx 'm';
dess evdu TITR '|U - uref|**2';
Finsi;
du = extr 1 (somme evdu);
errn = (abs (du - 4.06931E-03));
mess ' delta u = ' du ' errn = ' errn;
Si (NON COMPLET);
Si(errn > 0.0001); erreur 5 ; Finsi ;
Finsi;
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer