* fichier : recircp.dgibi
opti dime 2 elem QUA8 ;
DISCR = 'QUAF';
KPRES= 'CENTREP1' ;
GRAPH=FAUX ;
*GRAPH=VRAI ;
nbit = 5 ;
*
* 1/ Exemple test : canal courbe
* Eclt laminaire incompressible
* On impose une vitesse a un bout
* V normale nulle sur les parois
* On indique comment imposer une Vitesse en entree d'un domaine
* (Vitesse normale non nulle et Vitesse tangentielle nulle !!)
* et comment imposer un glissement sur une paroi (Vitesse normale nulle)
* Post : debit sur le contour , Vitesses ,pression
* Teste : vitesse normale imposee en 2D plan
* algorithme de projection transitoire
*
DISCR='QUAF' ;
KPRES='CENTREP1' ;
R=4. ; n1=5 ; n2=20 ;
p1=0 0 ; p2= 1 0 ; pc= R 0 ;
entree = p1 d n1 p2 ;
sortie = entree tour pc (-60.) ;
sortie = inve sortie;
q2 = poin 1 sortie ;
q1 = sortie poin final ;
pp1=p1 c n2 pc q1 ;
pp2=p2 c n2 pc q2 ;
paroi= pp1 et pp2 ;
cnt= entree et pp2 et sortie et (inve pp1) ;
mt= surf cnt;
Mmt= chan quaf mt ;
Mcnt= chan quaf cnt ;
Mentree= chan quaf entree;
Msortie= chan quaf sortie;
Mparoi = chan quaf paroi ;
elim (Mmt et Mcnt et Mentree et Msortie et Mparoi) 1.e-5 ;
$mt = mode Mmt 'NAVIER_STOKES' DISCR ;
mt= doma $mt maillage ;
$entree= mode Mentree 'NAVIER_STOKES' DISCR ;
$sortie= mode Msortie 'NAVIER_STOKES' DISCR ;
$paroi = mode Mparoi 'NAVIER_STOKES' DISCR ;
NU= 1.e-2 ;
U0=1. ;
DT=2. ;
* Comment imposer un debit d'entree
* En general il faut imposer la composante normale a
* la vitesse debitante et les composantes tangentielles a ZERO !
* On n'utilise donc pas VNIMP :
* On extrait le champ des normales de la surface d'entree
* en prenant garde a l'orientation des normales
nsor=doma $sortie 'NORMALEV';
* On extrait les composantes nx ny etc
nx= exco 'UX' nsor ;
ny= exco 'UY' nsor ;
* on cree les componsantes cartesiennes du champ de vitesse
* qui seront imposees en condition limite type Dirichlet
usx=u0*nx ;
usy=u0*ny ;
* Sur les parois la condition de glissement s'obtient
* effectivement en imposant une vitesse normale nulle a l'aide de VNIMP
rv = eqex 'ITMA' 1 'NITER' 1 'OMEGA' 1.
OPTI EF 'CENTREE' 'IMPL' KPRES
'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'INCO' 'UN'
'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
'CLIM'
'UN' UIMP sortie usx 'UN' VIMP sortie usy
;
rv.'INCO'= TABLE 'INCO' ;
rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
RVP= EQEX 'OPTI' 'EF' KPRES
'ZONE' $mt 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
'ZONE' $paroi 'OPER' 'VNIMP' $mt 0. 'INCO' 'UN' 'PRES'
;
* rvp. 'METHINV' . 'FCPRECT'=100 ;
RV.'PROJ'= RVP ;
exec rv ;
* affichage vitesse
un=rv.inco.'UN' ;
ung=vect un 0.1 ux uy jaune ;
Si GRAPH ;trace ung mt;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 ;
dq= (abs q) - 1.99528E-04 ;
si (dq > 1.e-5 ) ; erreur 5 ; finsi ;
* affichage pression
* On extrait la bonne composante de la pression
* il y a aussi les multiplicateurs de Lagrange vitesse normale
pe=exco rv.inco.'PRESSION' 'PRES' ;
* On reaffecte les bon SPG au CHPOINT
pe= kcht $mt scal KPRES pe ;
* On transforme le CHPOINT en CHPOINT aux sommets pour tracer les iso
pn= elno $mt pe KPRES ;
Si GRAPH ;trace pn mt ;Finsi ;
mpn=(maxi pn); mess (mini pn) ;
dpn= abs (mpn - 1.7757);
mess ' MAXI PN = ' mpn 'dpn=' dpn ;
si ( dpn > 1.e-4 ) ; erreur 5 ; finsi ;
*
* 2/ 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
* Voir RECIRC.DGIBI pour l'algorithme Taylor Hood
* Post : debit sur le contour , fonction de courant , pression
* Teste : vitesse normale imposee en 2D plan
* algorithme de projection transitoire
*
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. ;
DT=1. ;
rv = eqex 'ITMA' 10 'NITER' 1 'OMEGA' 1.
OPTI EF 'CENTREE' 'IMPL' KPRES
'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'INCO' 'UN'
'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
'CLIM'
'UN' UIMP slib 1. 'UN' VIMP ab 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. ;
RVP= EQEX 'OPTI' 'EF' KPRES
'ZONE' $mt 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
CLIM PRES TIMP bcp 0. ;
* rvp. 'METHINV' . 'FCPRECT'=100 ;
RV.'PROJ'= RVP ;
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 + .20777) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
si ( dpsi > 1.e-3 ) ; erreur 5 ; finsi ;
*
* Calcul de la pression
*
pnc=kcht $mt scal KPRES (exco rv.inco.'PRESSION' 'PRES') ;
pn= (elno pnc $mt KPRES) ;
Si GRAPH ; trace pn Mmt ; Finsi ;
mpn=(maxi pn);
dpn= abs (mpn - .32650);
mess ' MAXI PN = ' mpn 'dpn=' dpn ;
si ( dpn > 5.e-4 ) ; erreur 5 ; finsi ;
*
* Troisieme 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. ;
DT=1. ;
rv = eqex 'ITMA' 10 'NITER' 1 'OMEGA' 1.
OPTI EF 'CENTREE' 'IMPL' KPRES
'ZONE' $MT 'OPER' 'NS' 1. 'UN' NU 'INCO' 'UN'
'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
'CLIM'
'UN' UIMP slib 1. 'UN' VIMP ab 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. ;
RVP= EQEX 'OPTI' 'EF' KPRES
'ZONE' $mt 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
CLIM PRES TIMP bcp 0. ;
RV.'PROJ'= RVP ;
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 + .18) ;
mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
si ( dpsi > 1.e-3 ) ; erreur 5 ; finsi ;
*
* Calcul de la pression
*
pnc=kcht $mt scal KPRES (exco rv.inco.'PRESSION' 'PRES') ;
pn= (elno pnc $mt KPRES) ;
Si GRAPH ; trace pn Mmt ; Finsi ;
mpn=(maxi pn);
dpn= abs (mpn - .27);
mess ' MAXI PN = ' mpn 'dpn=' dpn ;
si ( dpn > 5.e-2 ) ; erreur 5 ; finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer