Accueil | Cast3M

* 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 ;