* fichier : canalKL.dgibi
* Ce cas teste le modèle K-L bas Reynolds (Wolfshtein Yap)
* 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
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*--- Relecture d'un profil de vitesse asymptotique de Référence --------
Si (NON COMPLET);
zsr=prog 0.
2.22209E-03 9.51930E-03 2.15866E-02 3.81009E-02 5.87392E-02
8.31785E-02 0.11110 0.14217 0.17607 0.21248
0.25108 0.29155 0.33355 0.37676 0.42088
0.46556 0.51049 0.55534 0.59980 0.64353
0.68623 0.72755 0.76718 0.80480 0.84008
0.87270 0.90234 0.92868 0.95138 0.97014
0.98462 0.99450 0.99946 1. ;
usr=prog 0.
0.25980 0.60654 0.74057 0.81674 0.87078
0.91336 0.94860 0.97849 1.0042 1.0263
1.0454 1.0617 1.0752 1.0858 1.0935
1.0979 1.0988 1.0962 1.0903 1.0812
1.0692 1.0544 1.0369 1.0164 0.99266
0.96512 0.93291 0.89458 0.84739 0.78511
0.68969 0.48982 6.49844E-02 0. ;
Sinon;
option acquerir 'ucanal_data.dgibi';
acque aa*mot;
acque aa*mot;
acque aa*mot;
acque aa*mot;
acque aa*mot;
acque aa*mot;
aa=1.;
acque z*flottant u*flottant v*flottant;
usr=prog u; zsr = prog z;
REPETER BLOC1 501;
acque z*flottant u*flottant v*flottant;
usr=usr et (prog u);
zsr=zsr et (prog z);
FIN BLOC1 ;
Finsi;
evusr = evol manu zsr usr;
su=extr 1 (somme evusr);
ep = maxi zsr;
um=su/ep;
mess ' Vitesse moyenne ' um;
usr = usr * (1./um);
zsr = zsr * (1./ep);
evusr = evol manu zsr usr ;
Si GRAPH ; dess evusr TITR ' Profil de vitesse de référence'; Finsi;
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*--- Profil de Reichardt -----------------------------------------------
Ka = 0.41;
yplusr= prog 1. pas 0.5 1500.;
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;
Si GRAPH ;dess evr1 'LOGX' TITR 'Profil de Reichardt'; Finsi;
*--- Profil de Reichardt -----------------------------------------------
*-----------------------------------------------------------------------
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 =400 ;
DT = 0.01;
H=1. ;
L=40.*H;
dcl=H/120;
dm=H/20;
dcx=H/15.;
dmx=H/10.;
Sinon;
NBIT =7;
DT = 0.1;
H=1. ;
L=10.*H;
dcl=H/20;
dm=H/5;
dcx=H/10.;
dmx=H/5.;
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 ));
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);
lm = (lm1*ilm1)+(lm2*ilm2) + lma;
dparoi=(coor 2 mt);
dparoi=dparoi*(H - dparoi);
dparoi=dparoi+1.e-10;
Si GRAPH ;
dess (evol chpo dparoi sortie) TITR ' Distance à la paroi';
trace mt TITR ' Maillage ';
Finsi ;
ye= coor 2 entre1;
Rey = 10000.;
U0 = 1. ;
Ufr = U0 ;
*Ufr = ye*(H - ye)*U0*1.5*4. ;
uux=(cos ang) * U0;
uuy=(sin ang) * U0;
NU = 1./Rey;
NU = 3.14E-5;
YP = 1.e-2 ;
fnu=1.e10;
fnu=Nu ;
RV= EQEX 'OMEGA' 1. 'NITER' 1 ITMA NBIT 'FIDT' 1
'OPTI' 'EF' 'IMPL' KSUPG KPRES
ZONE $mt 'OPER' 'KEPSILON' 1. 'UN' 'NU' DT 'INCO' 'KN' 'EN'
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.
KN TIMP entre1 ((0.02*Ufr)**2.)
KN TIMP (par1 et par2) 0. EN TIMP (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'.'KN' = kcht $mt scal sommet 1.e-5 ;
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'.'NU'=NU;
rv.'LTPS' = PROG ;
rv.'Tps' = 1. ;
rv.'DT' = DT;
rv.'NUPAT' =1;
mess ' EXECUTION '; exec rv ;
un = rv.inco.'UN' ;
Si GRAPH;
unx =exco un 'UX';
dudy=exco (kops 'GRADS' unx $mt) 'UY';
uet=(nu*(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';
mess ' Nb de pts ' (nbno entre1);
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*(1./ rv.inco.'NU');
yplus2= yp2*uet2*(1./ rv.inco.'NU');
evyplus1=evol chpo yplus1 par1;
evyplus2=evol chpo yplus2 (inve par2);
dess (evyplus1 et evyplus2) TITR ' Y+ le long des deux parois';
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 1 (extr evuet2 'ORDO');
mess ' uet1 et uet2=' uet1 uet2;
ypls1= evol chpo ((1./nu) * uet1 * (coor 2 sort1)) (inve sort1);
ypls2= evol chpo ((1./nu) * 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');
evupl1=evol manu xypls1 yypls1;
evupl2=evol manu (extr ypls2 'ORDO')(extr upl2 'ORDO');
evupl1=evupl1 coul rouge;
evupl2=evupl2 coul verte;
TAB1 = TABLE ;
tab1.1 = 'MOT' 'TIRC';
tab1.2 = 'MOT' 'TIRR';
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 (evupl2 et evupl1 et evr1) 'LOGX' 'XBOR' 1. 500.
TITR ' U+ en sortie' 'GRIL' ;
nutsnu=(rv.inco.'MUF')*(1./NU);
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';
Finsi ;
qs=dbit un $sortie IMPR;
um=(1./H)*qs;
unn=(1./um)*un;
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 ;
dess (evus et evusr) TITR 'U en sortie';
dess evdu TITR '|U - uref|**2';
Finsi;
du = extr 1 (somme evdu);
errn = (abs (du - 5.89880E-03));
mess ' delta u = ' du ' errn = ' errn;
Si (NON COMPLET);
Si(errn > 0.0011); erreur 5 ; Finsi ;
Finsi;
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer