* fichier : test1fpu.dgibi
* Ce cas vise à tester la fonction de paroi intégrée dans le maillage
* sur l'écoulement turbulent dans un canal plan. On utilise le modèle
* de Buleev.
*- 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 | |
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*--- 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 reference'; Finsi;
*-----------------------------------------------------------------------
*-----------------------------------------------------------------------
*--- Profil de Reichardt -----------------------------------------------
'SI' FAUX;
Ka = 0.41;
yplus = prog 1. pas 0.5 1500.;
unp = prog (dime yplus)*1.;
E=9.5;
*uplus=(1./Ka)*(Log ((E*yplus)));
uplus=(1./Ka)*(Log (unp + (Ka*yplus)));
uplus=uplus + (7.8* (unp - (exp (yplus*(-1./11))))
- ( (yplus*(-1./11)) * (exp (yplus*(-1./3.))) ) );
b=4.;
a=1./b;
uplus2=(1./Ka)*(Log (unp + (Ka*a*yplus)));
uplus2=uplus2 + (7.8* (unp - (exp (a*yplus*(-1./11))))
- ( (a*yplus*(-1./11)) * (exp (a*yplus*(-1./3.))) ) );
uplus2=b*uplus2;
b=2.;
a=1./b;
uplus3=(1./Ka)*(Log (unp + (Ka*a*yplus)));
uplus3=uplus3 + (7.8* (unp - (exp (a*yplus*(-1./11))))
- ( (a*yplus*(-1./11)) * (exp (a*yplus*(-1./3.))) ) );
uplus3=b*uplus3;
evr1=evol manu yplus uplus;
evr2=evol manu yplus uplus2;
evr3=evol manu yplus uplus3;
dess (evr1 et evr2 et evr3) 'LOGX' ;
*opti donn 5 ;
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;
DEBPROC CALCUL;
ARGU RVX*'TABLE' ;
RV = RVX.'EQEX' ;
* Tps = rv.'Tps' ;
* DT = rv.'DT' ;
* NUPAT = rv.'NUPAT' + 1;
* Tps = Tps + DT ;
* rv.'Tps' = Tps ;
* rv.'NUPAT'=NUPAT ;
un = rv.'INCO'.'UN' ;
uux= exco un 'UX';
dudy = exco (kops 'GRADS' uux $mt) 'UY' ;
anut=(lm*lm * (abs dudy));
rv.inco.'NUT'= kcht $mt scal sommet anut;
as2 ama1 = 'KOPS' 'MATRIK' ;
finproc as2 ama1 ;
******=================================================**********
******=================================================**********
******=================================================**********
******=================================================**********
******=================================================**********
Si COMPLET;
NBIT =200 ;
H=1. ;
L=40.*H;
dcl=H/50.;
dm=H/10.;
Sinon;
NBIT =20 ;
H=1. ;
L=10.*H;
dcl=H/4.;
dm =H/5.;
Finsi;
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 dcl dfin dm 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';
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;
*lm = lm * (1. - (exp ((-1.)*(coor 1 mt))));
Si GRAPH ; trace mt; 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;
DT = 0.5;
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' 'CALCUL'
ZONE $par1 'OPER' 'FPU' 1. 'UN' 'NU' 'UET1' $mt 'NUT' 'INCO' 'UN'
ZONE $par2 'OPER' 'FPU' 1. 'UN' 'NU' 'UET2' $mt 'NUT' 'INCO' 'UN'
ZONE $mt 'OPER' 'NS' 1. 'UN' 'NUT' 'INCO' 'UN'
'OPTI' 'EFM1' 'IMPL' 'CENTREE'
ZONE $mt OPER DFDT 1. 'UNM' DT INCO 'UN';
;
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=450;
rv.'METHINV'.PRECOND=3 ;
*rv.'METHINV'.PRECOND=2 ;
rv.'METHINV'.ILUTLFIL=4 ;
rv.'METHINV'.ILUTDTOL=0.;
rv.'METHINV'.RESID =1.e-6;
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=500;
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'.'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.'INCO'.'NUT'= kcht $mt scal sommet (Nu*500.) ;
rv.'LTPS' = PROG ;
rv.'Tps' = 1. ;
rv.'DT' = DT;
rv.'NUPAT' =1;
*- trace mt;
*- trace mt;
mess ' EXECUTION '; exec rv ;
un = rv.inco.'UN' ;
Si GRAPH;
uet1 = rv.INCO.'UET1';
uet2 = rv.INCO.'UET2';
evuet1= evol chpo uet1 par1;
evuet2= evol chpo uet2 (inve par2);
dess (evuet1 et evuet2) TITR ' U* ';
yp1= rv.2FPU.'MYP' ;
yp2= rv.3FPU.'MYP' ;
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+ ';
trace rv.inco.'NUT' mt (cont mt);
ung = vect un 0.5 'UX' 'UY' jaune ;
trace ung mt ;
pn = elno $mt rv.inco.PRESSION kpres;
trace pn mt;
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 - 4.12346E-02));
mess ' delta u = ' du ' errn = ' errn;
Si(errn > 0.002); erreur 5 ; Finsi ;
FIN;
© Cast3M 2003 - All rights reserved.
Disclaimer