* fichier : hy1.dgibi
*$$$ HY1
** Exemple HY1
** ___________
**
** --- 2 OCTOBRE 1990 ---
**
** CANAL LONGUEUR 10. LARGEUR 1.
** test cas isotherme NS et TOIMP
**
** On considère l'ecoulement de poiseuille dans un canal plan vertical
** compris entre les plans x=0 et x=1. la hauteur est de 10.
** On maintient entre l'entrée et la sortie une différence de pression.
** En regime établi le profil de vitesse ne depend que de x et la
** pression que de y : V(x)=4Vm X(1-X) Vm vitesse max (x=1/2) et
** Vmoy=2/3Vm Vm=-1/8/Mu DP/DY Mu viscosité dynamique
** T0=Mu DV/DX =4Mu Vm et le Reynolds Re=1/12/Mu/Mu DP/DY
** (RO=1 Diam=1)
** On a choisi Mu=5.e-2 et DP/DY=0.6 => Re=20 Vm=1.5 Vmoy=1
** Avec ces valeurs le temps d'établissement du régime est dominé par
** la diffusion. Il faut 12s pour arriver au régime établi (Vmoy=0.99)
** soit approximativement 800 pas de temps avec des QUA8 et le double avec
** des TRI6.
** Si on réduit la viscosité (d'un facteur 10 par exemple) la vitesse
** augmente pour une même perte de charge et le temps d'établissement
** du régime est dominé‚par la convection.
** Dans le cas ou le canal est coudé (IKAS=1) on peut noter une
** différence sur les isobares lorsque la viscosité est faible
** (Mu=5.e-3). Ceci est du à l'effet d'inertie dans le coude.
**
** On peut tester aussi le cas ou on impose la tension en paroi.
**
GRAPH = 'N' ;
COMPLET = FAUX ;
SI ( COMPLET ) ;
nbit3=2150 ;
err1 =1.e-5 ;
qe13= 0.9922268 ;
qe14= 0.9922544 ;
dq1 =1.e-14 ;
dpes1 = 5.6997 ;
erpe1 = 1.e-3 ;
nbit4=1000 ;
tfin=12. ;
SINON ;
nbit3=180 ;
err1 =1.e-5 ;
err1 =1.e-5 ;
qe14 =0.3989106 ;
qe13 =0.3956187 ;
dq1 =1.e-14 ;
dpes1 = 5.6997 ;
erpe1 = 1.e-3 ;
nbit4=80 ;
tfin=1. ;
FINSI ;
typelt=mot qua8 ;
nbe=7 ; nbv=10 ;
angle=0. ;
nitma=nbit4 ;
nitma=100 ;
tfinal=tfin ;
DEBPROC HY1 ;
ARGU TYPELT*MOT NBE*ENTIER NBV*ENTIER ANGLE*FLOTTANT
NITMA*ENTIER TFINAL*FLOTTANT GRAPH*MOT ;
OPTION DIME 2 ELEM TYPELT ;
p1=0 0.;
p2=1. 0.;
VH=0 10;
entree= p1 d nbe p2 ;
c1= 5. 0. ;
ikas=0 ;
*mess 'ikas= 0 --> DROIT (option par defaut) ikas=1 --> COURBE ? ';
*obtenir ikas*entier ;
si (EGA ikas 1) ;
q1=p1 tour c1 -90. ;
q2=p2 tour c1 -90. ;
pp1=p1 c c1 q1 nbv;
pp2=p2 c c1 q2 nbv;
sortie=entree tour c1 -90 ;
sinon ;
q1=p1 plus vh ;
q2=p2 plus vh ;
pp1=p1 d q1 nbv;
pp2=p2 d q2 nbv;
sortie=entree plus vh ;
finsi ;
pp1=inve pp1 ;
sortie=inve sortie;
elim 0.0001 (sortie et pp1 et pp2 et entree );
cnt=entree et pp2 et sortie et pp1 ;
*mt=surf cnt ;
mt=daller entree pp2 sortie pp1 ;
*pp1=inve pp1 ; entree=inve entree ; sortie=inve sortie ;
*pp2=inve pp2 ;
*mt=daller entree pp1 sortie pp2 ;
*entree=inve entree ;
mt=mt tour c1 angle ;
entree=entree tour c1 angle ;
sortie=sortie tour c1 angle ;
pp1=pp1 tour c1 angle ;
pp2=pp2 tour c1 angle ;
elim (pp1 et pp2 et entree et sortie et mt) 0.0001 ;
bell = chan mt quaf ;
entree = chan entree quaf ;
sortie = chan sortie quaf ;
pp1 = chan pp1 quaf ;
pp2 = chan pp2 quaf ;
elim (mt et entree et sortie et pp1 et pp2) 1.e-4 ;
$bell = mode bell 'NAVIER_STOKES' MACRO ; doma $bell 'IMPR' ;
$entree = mode entree 'NAVIER_STOKES' MACRO ; doma $entree 'IMPR' ;
$sortie = mode sortie 'NAVIER_STOKES' MACRO ; doma $sortie 'IMPR' ;
$pp1 = mode pp1 'NAVIER_STOKES' MACRO ; doma $pp1 'IMPR' ;
$pp2 = mode pp2 'NAVIER_STOKES' MACRO ; doma $pp2 'IMPR' ;
nble=(nbel (doma $entree 'MAILLAGE') ) - 1 ;
entref=elem (doma $entree 'MAILLAGE') (lect 2 pas 1 nble) ;
*?$entref=mode entref 'NAVIER_STOKES' MACRO ;doma $entref 'IMPR' ;
mu=5.E-2 ;
ro=1 ;
nu=mu/ro ;
to = kcht $entree vect centre ( 0. 6.) ;
tos = kcht $sortie vect centre ( 0. 0.) ;
tt1 = kcht $pp1 vect centre (-0.03 0.) ;
tt2 = kcht $pp2 vect centre (-0.03 0.) ;
dt=1.e-2 ;
rv=eqex $bell 'DUMP' ALFA 0.3 ITMA nitma 'TFINAL' tfinal
ZONE $BELL OPER NS NU INCO 'UN'
OPTI 'CENTREE'
ZONE $BELL OPER DFDT 1. 'UN' 'DELTAT ' INCO 'UN'
ZONE $ENTREE OPER TOIMP TO INCO 'UN'
ZONE $SORTIE OPER TOIMP TOS INCO 'UN'
**ZONE $ENTREE OPER TOIMP TO 'UN' nu INCO 'UN'
**ZONE $SORTIE OPER TOIMP TOS 'UN' nu INCO 'UN'
* ZONE $PP1 OPER TOIMP TT1 INCO 'UN'
* ZONE $PP2 OPER TOIMP TT2 INCO 'UN' ;
;
rvp= eqpr $bell
zone $PP1 oper VNIMP 0.
zone $bell oper PRESSION 0.
zone $PP1 oper VTIMP 0.
zone $PP2 oper VNIMP 0.
zone $PP2 oper VTIMP 0.
* zone $ENTREF oper VNIMP 1.
* zone $ENTREF oper VTIMP 0.
;
rv.'INCO'=table 'INCO' ;
rv.'INCO'.'UN' = kcht $bell VECT SOMMET (0 0) ;
rv.pression=rvp ;
lh= (noeu 10) et (noeu 20) et (noeu 30) et (noeu 40) et
(noeu 50) et (noeu 60) ;
lj= (manu poi1 ((doma $bell 'MAILLAGE') poin proc( 0.5 0.5) ) ) ;
his = khis 'UN' 1 lh 'UN' 2 (lh et lj) ;
rv.'HIST'=his ;
exec rv ;
pn=elno (kcht $bell scal centre (rvp.pression)) $bell ;
ung1 = vect 0.5 (rv.'INCO'.'UN') ux uy jaune ;
qe=dbit (rv.'INCO'.'UN') $entree ;
qs=dbit (rv.'INCO'.'UN') $sortie ;
dq=(abs qs )-(abs qe) ;
mess (' Bilan : dq=') dq ;
pet=somt (redu pn entree) / (nbno entree) ;
pst=somt (redu pn sortie) / (nbno sortie) ;
dpes=pet - pst ;
mess ' pet pst dpes = ' pet pst dpes ;
si ('EGA' graph 'O' );
dessin (his.'TABD') (his.'1UN') ;
dessin (his.'TABD') (his.'2UN') ;
trace ung1 bell ;
trace pn bell (cont bell) ;
finsi ;
FINPROC RV dq qe dpes ;
type=mot tri6 ;
* mess 'Type d élément TRI6 ou QUA8 ? ' ;
* obtenir TYPE*mot ;
angle=0. ;
* obtenir angle*flottant ;
nbe=7 ; nbv=10 ;
type=mot qua8 ;
rv dq qe dpes = hy1 TYPE nbe nbv angle nbit4 tfin graph ;
erq = abs (qe - qe14) ;
mess ' Erreur sur le debit ' erq ;
si ( erq > err1 ) ; erreur 5 ; finsi ;
dq=abs dq ;
si ( dq > dq1 ) ; erreur 5 ; finsi ;
* On teste quand meme le delta P qui viendrait d'une erreur de discretisation
ddp=abs ( dpes - dpes1) ;
mess ' ddp erpe1 ' ddp erpe1 ;
si ( ddp > erpe1 ) ; erreur 5 ; finsi ;
*opti donn 5 ;
type=mot tri6 ;
rv dq qe dpes = hy1 TYPE nbe nbv angle nbit3 tfin graph ;
erq = abs (qe - qe13) ;
mess ' Erreur sur le debit ' erq ;
si ( erq > err1 ) ; erreur 5 ; finsi ;
dq=abs dq ;
si ( dq > dq1 ) ; erreur 5 ; finsi ;
FIN ;
© Cast3M 2003 - All rights reserved.
Disclaimer