* fichier : recirc.dgibi
*
* Exemple test : recirculation dans une cavite semi circulaire
* Eclt laminaire incompressible
* On impose une vitesse tg sur la surface libre
* V normale nulle sur le fond
*
* Post : debit sur le contour , fonction de courant , pression
* Teste vitesse normale imposee en 2D plan algorithme Taylor Hood
*
* Voir aussi RECIRCP.DGIBI pour algorithme de projection
opti dime 2 elem QUA8 ;
DISCR = 'QUAF';
KPRES= 'CENTREP1' ;
GRAPH=FAUX ;
*GRAPH=VRAI ;
nbit = 5 ;
p0=0. 0. ;
p1= -1. 0.;
p2= 1. 0. ;
p03= 0. -1. ;
p3= 1. -1. ;
p4= -1. -1. ;
n1= 5 ; n2= 5 ;
ab= p1 d n1 p2 ;
slib= (chan ab POI1) ;
nab= nbel slib ;
slib= elem slib (lect 2 pas 1 (nab - 1));
bc= p2 c p0 n2 p03 c p0 n2 p1;
*bc= p2 d n2 p3 d n2 p4 d n2 p1 ;
cnt=ab et bc ;
mt= cnt surf ;
Mmt=chan mt quaf ;
Mbc=chan bc quaf ;
elim (Mmt et Mbc) 1.e-5 ;
$mt = mode Mmt 'NAVIER_STOKES' DISCR;
$bc = mode Mbc 'NAVIER_STOKES' DISCR;
* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt KPRES ;
bcp=elem prep1 POI1 (lect 1) ;
*
*NU= 1. / 61. ;
NU= 1.e-2 ;
U0=1. ;
rv = eqex 'ITMA' 0 'NITER' nbit 'OMEGA' 0.8
OPTI EF 'CENTREE' 'IMPL' KPRES
'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'INCO' 'UN'
'ZONE' $MT 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PRES'
'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
'CLIM'
'UN' UIMP slib 1. 'UN' VIMP ab 0.
'PRES' TIMP bcp 0. ;
rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
rv.'INCO'.'VN'=kcht $bc scal sommet 0. ;
exec rv ;
un = rv.inco.'UN' ;
ung= vect un 0.11 ux uy jaune ;
Si GRAPH ;trace ung Mmt ; Finsi ;
cnt1=chan cnt seg2 ;
Mcnt1= chan cnt1 quaf;
$cnt=model Mcnt1 'NAVIER_STOKES' line ;
q=dbit (rv.inco.'UN') $cnt ;
Mess ' Debit sur le contour ' q ;
si ( q > 1.e-15 ) ; erreur 5 ; finsi ;
*
* Calcul de la fonction de courant
*
sw = 'KOPS' un 'ROT' $mt ;
rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
'ZONE' $mt 'OPER' 'LAPN' 1.0D0 'INCO' 'PSI'
'ZONE' $mt 'OPER' 'FIMP' sw 'INCO' 'PSI'
'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;
opti isov SULI ;
exec rk ;
psi=rk.inco.'PSI' ;
Si GRAPH ;trace psi Mmt cnt ; Finsi ;
mpsi=(mini psi);
dpsi= abs (mpsi + .26820) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
si ( dpsi > 1.e-5 ) ; erreur 5 ; finsi ;
*
* Calcul de la pression
*
pn= elno (rv.inco.'PRES') $mt KPRES;
Si GRAPH ; trace pn Mmt ; Finsi ;
*
* Deuxieme cas : Meme probleme geometrie rectangulaire
* MACRO CENTRE
*
*
*
*
*
*
DISCR = 'MACRO';
KPRES= 'CENTRE' ;
p0=0. 0. ;
p1= -1. 0.;
p2= 1. 0. ;
p03= 0. -1. ;
p3= 1. -1. ;
p4= -1. -1. ;
n1= 5 ; n2= 5 ;
ab= p1 d n1 p2 ;
slib= (chan ab POI1) ;
nab= nbel slib ;
slib= elem slib (lect 2 pas 1 (nab - 1));
bc= p2 d n2 p3 d n2 p4 d n2 p1 ;
cnt=ab et bc ;
mt= cnt surf ;
Mmt=chan mt quaf ;
Mbc=chan bc quaf ;
elim (Mmt et Mbc) 1.e-5 ;
$mt = mode Mmt 'NAVIER_STOKES' DISCR;
$bc = mode Mbc 'NAVIER_STOKES' DISCR;
* La cavité est fermée il faut imposer la pression en 1 point !
prep1=doma $mt KPRES ;
bcp=elem prep1 POI1 (lect 1) ;
*
NU= 1.e-2 ;
U0=1. ;
rv = eqex 'ITMA' 0 'NITER' nbit 'OMEGA' 0.8
OPTI EF 'CENTREE' 'IMPL' KPRES
'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'INCO' 'UN'
'ZONE' $MT 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PRES'
'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
'CLIM'
'UN' UIMP slib 1. 'UN' VIMP ab 0.
'PRES' TIMP bcp 0. ;
rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
rv.'INCO'.'VN'=0. ;
exec rv ;
un = rv.inco.'UN' ;
ung= vect un 0.11 ux uy jaune ;
Si GRAPH ;trace ung Mmt ; Finsi ;
cnt1=chan cnt seg2 ;
Mcnt1= chan cnt1 quaf;
$cnt=model Mcnt1 'NAVIER_STOKES' line ;
q=dbit (rv.inco.'UN') $cnt ;
Mess ' Debit sur le contour ' q ;
si ( q > 1.e-15 ) ; erreur 5 ; finsi ;
sw = 'KOPS' un 'ROT' $mt ;
rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
'ZONE' $mt 'OPER' 'LAPN' 1.0D0 'INCO' 'PSI'
'ZONE' $mt 'OPER' 'FIMP' sw 'INCO' 'PSI'
'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
rk . 'INCO' = 'TABLE' 'INCO' ;
rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;
opti isov SULI ;
exec rk ;
psi=rk.inco.'PSI' ;
mess ' MINI PSI = ' (mini psi) ;
Si GRAPH ;trace psi Mmt cnt ; Finsi ;
mpsi=(mini psi);
dpsi= abs (mpsi + .22) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
si ( dpsi > 1.e-3 ) ; erreur 5 ; finsi ;
*
* Calcul de la pression
*
pn= elno (rv.inco.'PRES') $mt KPRES;
Si GRAPH ; trace pn Mmt ; Finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer