* fichier : coudep.dgibi
* Ecoulement dans un Coude
* Test 3D pression continue methode de projection
* teste en 3D VNIMP FPU NS
* Ce test pose des problemes en QUAF ??
* Il faut un minimum de mailles pour converger a fort Reynolds
*
opti dime 2 elem QUA8 ;
GRAPH=VRAI ;
GRAPH=FAUX ;
COMPLET = VRAI ;
COMPLET = FAUX ;
KPRES='MSOMMET' ;
Si COMPLET ;
R=0.5 ; nl = 4 ;n1= 5 ; lgd1=1.; lgd2=2.; n2=5 ;
NU=1.e-5 ;
NUT=5.e-3;
YP=5.e-3 ;
DT=1. ;
Sinon ;
R=0.5 ; nl = 2 ;n1= 2 ; lgd1=1.; lgd2=2.; n2=2 ;
NU=1.e-5 ;
NUT=5.e-2;
YP=5.e-2 ;
DT=1. ;
Finsi ;
P0=0 0 ; p1= R 0 ;
p2 = p1 tour 90 p0 ;
p3 = p2 tour 90 p0 ;
p4 = p3 tour 90 p0 ;
mct= p1 c p0 n1 p2 c p0 n1 p3 c p0 n1 p4 c p0 n1 p1;
ct1= p1 c p0 n1 p2;
ct2= p2 c p0 n1 p3;
ct3= p3 c p0 n1 p4;
ct4= p4 c p0 n1 p1;
*ct2= p2 d n1 p3;
*ct4= p4 d n1 p1;
*mct= p1 c p0 n1 p2 d n1 p0 d n1 p1;
*mc= surf mct ;
mc= daller ct1 ct2 ct3 ct4 ;
*trace mc ;
opti dime 3 elem seg2 ;
q1= p0 plus (0 0 lgd1);
ax1 = q1 plus (1 0 0) ;
ax2 = q1 plus (1 1 0) ;
q2= q1 tour 90. ax1 ax2 ;
q3= q2 plus (lgd2 0 0);
lg= p0 d nl q1 c ax1 n2 q2 d nl q3 ;
*lg= p0 d nl q1 ;
opti elem CU20 ;
*trace lg ;
lg1= chan lg 'POI1' ;
nbp= (nbno lg1) - 2 ;
pi1= p0 ;
si1= mc ;
pi2= poin 2 lg1 ;
vh1=pi2 moins pi1 ;
si2= si1 plus vh1 ;
Entree = si1 ;
vc=si1 volu 2 si2;
*trace cache vc ;
repeter bloc1 nbp ;
pi1=pi2;
si1=si2 ;
pi2= poin (&bloc1 + 2) lg1 ;
vh2=pi2 moins pi1 ;
vhr=vh2 ;
axe =vh1 pvec vhr ;
alfa= atg (norm(axe / (vh1 psca vhr) ));
Si ( alfa >EG 1.) ;
si2= si1 tour alfa pi1 (pi1 plus axe) ;
Sinon ;
si2= si1 ;
Finsi ;
si2=si2 plus vh2 ;
vh1=vh2 ;
vc = vc et (si1 volu 2 si2);
fin bloc1 ;
Sortie = si2 ;
oeil = 0. -1.e5 -100. ;
*trace cache vc oeil ;
Sortie= Sortie coul rouge ;
Entree= Entree coul rouge ;
Mvc = chan vc QUAF ;
Mentree = chan entree QUAF ;
Msortie = chan sortie QUAF ;
ELIM (Mvc et Mentree et Msortie) 1.e-5 ;
DEBPROC TEST DISCR*MOT ;
$vc = mode Mvc 'NAVIER_STOKES' DISCR ;
$sortie = mode Msortie 'NAVIER_STOKES' DISCR ;
sortip=doma $sortie KPRES;
$entree = mode Mentree 'NAVIER_STOKES' DISCR ;
entree = doma $entree maillage ;
Menvlp = doma $vc 'ENVELOPPE' ;
Mparoi = diff Menvlp (Mentree et Msortie) ;
doma $vc 'IMPR' ;
* On extrait le champ des normales de la surface d'entree
* en prenant garde a l'orientation des normales
nentr= doma $entree 'NORMALEV';
Mcte = doma $entree 'ENVELOPPE';
$cte= mode Mcte 'NAVIER_STOKES' DISCR ;
cte=doma $cte maillage ;
* On extrait les composantes nx ny etc
nx= exco 'UX' nentr;
ny= exco 'UY' nentr;
nz= exco 'UZ' nentr;
* on cree les componsantes cartesiennes du champ de vitesse
* qui seront imposees en condition limite type Dirichlet
U0=1. ;
usx=u0*nx ;
usy=u0*ny ;
usz=u0*nz ;
$envlp = mode Menvlp 'NAVIER_STOKES' DISCR ;
envlp=doma $envlp maillage ;
$paroi = mode Mparoi 'NAVIER_STOKES' DISCR ;
paroi=doma $paroi maillage ;
rv= eqex 'ITMA' 5 'NITER' 1 'OMEGA' 1.
OPTI EF 'SUPG' 'IMPL' KPRES
'ZONE' $paroi 'OPER' 'FPU' 1. 'UN' NU 'UET' YP 'INCO' 'UN'
'ZONE' $vc 'OPER' 'NS' 1. 'UN' NUT 'INCO' 'UN'
OPTI EF 'CENTREE'
'ZONE' $vc 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
'CLIM'
'UN' UIMP entree usx 'UN' VIMP entree usy 'UN' WIMP entree usz
* 'UN' UIMP paroi 0. 'UN' VIMP paroi 0. 'UN' WIMP paroi 0.
;
rv.'METHINV'.TYPINV=3 ;
rv.'METHINV'.IMPINV=0 ;
rv.'METHINV'.NITMAX=100;
rv.'METHINV'.PRECOND=3 ;
rv.'METHINV'.RESID =1.e-10 ;
rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $vc vect sommet (0. 0. 0.) ;
rv.'INCO'.'PRES'=kcht $vc scal KPRES 0. ;
RVP= EQEX 'OPTI' 'EF' KPRES
'ZONE' $vc 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
'ZONE' $paroi 'OPER' 'VNIMP' $vc 0. 'INCO' 'UN' 'PRES'
'CLIM' 'PRES' TIMP sortip 0.
;
rvp.'METHINV'.TYPINV=3 ;
rvp.'METHINV'.IMPINV=0 ;
rvp.'METHINV'.NITMAX=100;
rvp.'METHINV'.PRECOND=3 ;
rvp.'METHINV'.RESID =1.e-10 ;
RV.'PROJ'= RVP ;
exec rv ;
un=rv.inco.'UN' ;
un1=redu un envlp ;
ung=vect un1 0.1 ux uy uz jaune ;
Si GRAPH ;
trace ung envlp oeil ;
Finsi;
q= dbit un $envlp ;
list q ;
qe=dbit un $entree ;
qs=dbit un $sortie ;
qp=dbit un $paroi ;
mess 'qe=' qe 'qs=' qs 'qp=' qp 'qs+qp=' (qs+qp) ;
yplus=yp*(rv.inco.'UET')/nu ;
mess 'Mini yplus ' (mini yplus) ' Maxi yplus ' (maxi yplus) ;
pn=exco rv.inco.'PRESSION' 'PRES' ;
mtp= doma $vc 'MMAIL' ;
Si GRAPH ;
trace pn mtp ;
Finsi;
mess 'Mini Pres ' (mini pn) ' Maxi Pres ' (maxi pn) ;
FINPROC ;
*test 'QUAF' ;
test 'MACRO';
test 'LINE' ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer