Télécharger la plaquette Cast3M


fr en

Télécharger hy1.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  hy1.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *$$$ HY1
  5.  
  6.  
  7. **  Exemple HY1
  8. **  ___________
  9. **
  10. **  --- 2 OCTOBRE 1990 ---
  11. **
  12. **   CANAL LONGUEUR 10. LARGEUR 1.
  13. **   test cas isotherme NS et TOIMP
  14. **
  15. **  On considère l'ecoulement de poiseuille dans un canal plan vertical
  16. **  compris entre les plans x=0 et x=1. la hauteur est de 10.
  17. **  On maintient entre l'entrée et la sortie une différence de pression.
  18. **  En regime établi le profil de vitesse ne depend que de x et la
  19. **  pression que de y : V(x)=4Vm X(1-X)  Vm vitesse max (x=1/2) et
  20. **  Vmoy=2/3Vm  Vm=-1/8/Mu DP/DY   Mu viscosité dynamique
  21. **  T0=Mu DV/DX =4Mu Vm  et le Reynolds Re=1/12/Mu/Mu DP/DY
  22. **  (RO=1 Diam=1)
  23. **  On a choisi Mu=5.e-2 et DP/DY=0.6 => Re=20 Vm=1.5 Vmoy=1
  24. **  Avec ces valeurs le temps d'établissement du régime est dominé par
  25. **  la diffusion. Il faut 12s pour arriver au régime établi (Vmoy=0.99)
  26. **  soit approximativement 800 pas de temps avec des QUA8 et le double avec
  27. **  des TRI6.
  28. **  Si on réduit la viscosité (d'un facteur 10 par exemple) la vitesse
  29. **  augmente pour une même perte de charge et le temps d'établissement
  30. **  du régime est dominé‚par la convection.
  31. **  Dans le cas ou le canal est coudé (IKAS=1) on peut noter une
  32. **  différence sur les isobares lorsque la viscosité est faible
  33. **  (Mu=5.e-3). Ceci est du à l'effet d'inertie dans le coude.
  34. **
  35. **  On peut tester aussi le cas ou on impose la tension en paroi.
  36. **
  37.  
  38. GRAPH = 'N' ;
  39.  
  40. COMPLET = FAUX ;
  41. SI ( COMPLET ) ;
  42.  nbit3=2150  ;
  43.  err1 =1.e-5 ;
  44.  qe13= 0.9922268 ;
  45.  qe14= 0.9922544 ;
  46.  dq1  =1.e-14 ;
  47.  dpes1 = 5.6997   ;
  48.  erpe1 = 1.e-3    ;
  49.  nbit4=1000  ;
  50.  tfin=12. ;
  51. SINON        ;
  52.  nbit3=180   ;
  53.  err1 =1.e-5 ;
  54.  err1 =1.e-5 ;
  55.  qe14 =0.3989106 ;
  56.  qe13 =0.3956187 ;
  57.  dq1  =1.e-14 ;
  58.  dpes1 = 5.6997  ;
  59.  erpe1 = 1.e-3    ;
  60.  nbit4=80    ;
  61.  tfin=1. ;
  62. FINSI ;
  63.  
  64. typelt=mot qua8 ;
  65. nbe=7 ; nbv=10 ;
  66. angle=0. ;
  67. nitma=nbit4 ;
  68. nitma=100   ;
  69. tfinal=tfin ;
  70.  
  71.  
  72.  DEBPROC HY1 ;
  73.  ARGU TYPELT*MOT NBE*ENTIER NBV*ENTIER ANGLE*FLOTTANT
  74.      NITMA*ENTIER TFINAL*FLOTTANT GRAPH*MOT ;
  75.  
  76. OPTION DIME 2 ELEM TYPELT ;
  77. p1=0  0.;
  78. p2=1. 0.;
  79. VH=0 10;
  80. entree= p1 d nbe p2 ;
  81. c1= 5. 0. ;
  82.  
  83. ikas=0 ;
  84. *mess 'ikas= 0 --> DROIT (option par defaut)  ikas=1 --> COURBE  ? ';
  85. *obtenir ikas*entier ;
  86.  
  87. si (EGA ikas 1) ;
  88.  
  89. q1=p1 tour c1 -90. ;
  90. q2=p2 tour c1 -90. ;
  91. pp1=p1 c c1 q1 nbv;
  92. pp2=p2 c c1 q2 nbv;
  93. sortie=entree tour c1 -90 ;
  94.  
  95. sinon ;
  96.  
  97. q1=p1 plus vh ;
  98. q2=p2 plus vh  ;
  99. pp1=p1 d q1 nbv;
  100. pp2=p2 d q2 nbv;
  101. sortie=entree plus vh ;
  102. finsi ;
  103.  
  104. pp1=inve pp1 ;
  105. sortie=inve sortie;
  106. elim 0.0001 (sortie et pp1 et pp2 et entree );
  107. cnt=entree et pp2 et sortie et pp1 ;
  108. *mt=surf cnt ;
  109.  mt=daller entree pp2 sortie pp1 ;
  110.  
  111. *pp1=inve pp1 ; entree=inve entree ; sortie=inve sortie ;
  112. *pp2=inve pp2 ;
  113. *mt=daller entree pp1 sortie pp2 ;
  114. *entree=inve entree ;
  115.  
  116. mt=mt tour c1 angle ;
  117. entree=entree tour c1 angle ;
  118. sortie=sortie tour c1 angle ;
  119. pp1=pp1 tour c1 angle ;
  120. pp2=pp2 tour c1 angle ;
  121.  
  122. elim (pp1 et pp2 et entree et sortie et mt) 0.0001 ;
  123.  
  124.  bell   = chan mt     quaf ;
  125.  entree = chan entree quaf ;
  126.  sortie = chan sortie quaf ;
  127.  pp1    = chan pp1    quaf ;
  128.  pp2    = chan pp2    quaf ;
  129.  elim (mt et entree et sortie et pp1 et pp2) 1.e-4 ;
  130.  
  131.  $bell   = mode bell   'NAVIER_STOKES' MACRO ; doma $bell 'IMPR' ;
  132.  
  133.  $entree = mode entree 'NAVIER_STOKES' MACRO ; doma $entree 'IMPR' ;
  134.  $sortie = mode sortie 'NAVIER_STOKES' MACRO ; doma $sortie 'IMPR' ;
  135.  $pp1    = mode pp1    'NAVIER_STOKES' MACRO ; doma $pp1    'IMPR' ;
  136.  $pp2    = mode pp2    'NAVIER_STOKES' MACRO ; doma $pp2    'IMPR' ;
  137.  
  138. nble=(nbel (doma $entree 'MAILLAGE') ) - 1 ;
  139. entref=elem (doma $entree 'MAILLAGE') (lect 2 pas 1 nble) ;
  140.  
  141. *?$entref=mode entref 'NAVIER_STOKES' MACRO ;doma $entref 'IMPR' ;
  142.  
  143.   mu=5.E-2 ;
  144.   ro=1 ;
  145.   nu=mu/ro ;
  146.  
  147.   to  = kcht $entree vect centre ( 0.  6.) ;
  148.   tos = kcht $sortie vect centre ( 0.  0.) ;
  149.  
  150.   tt1 = kcht $pp1    vect centre (-0.03 0.) ;
  151.   tt2 = kcht $pp2    vect centre (-0.03 0.) ;
  152.  
  153.   dt=1.e-2 ;
  154.   rv=eqex $bell 'DUMP' ALFA 0.3  ITMA nitma 'TFINAL' tfinal
  155.   ZONE  $BELL    OPER    NS    NU          INCO 'UN'
  156.   OPTI  'CENTREE'
  157.   ZONE  $BELL    OPER    DFDT  1. 'UN' 'DELTAT  ' INCO 'UN'
  158.   ZONE  $ENTREE  OPER    TOIMP TO          INCO 'UN'
  159.   ZONE  $SORTIE  OPER    TOIMP TOS         INCO 'UN'
  160. **ZONE  $ENTREE  OPER    TOIMP TO  'UN' nu INCO 'UN'
  161. **ZONE  $SORTIE  OPER    TOIMP TOS 'UN' nu INCO 'UN'
  162. * ZONE  $PP1     OPER    TOIMP TT1         INCO 'UN'
  163. * ZONE  $PP2     OPER    TOIMP TT2         INCO 'UN'  ;
  164.   ;
  165.  
  166.   rvp= eqpr $bell
  167.   zone $PP1    oper VNIMP 0.
  168.   zone $bell oper PRESSION  0.
  169.   zone $PP1    oper VTIMP 0.
  170.   zone $PP2    oper VNIMP 0.
  171.   zone $PP2    oper VTIMP 0.
  172. * zone $ENTREF oper VNIMP 1.
  173. * zone $ENTREF oper VTIMP 0.
  174.   ;
  175.  
  176.   rv.'INCO'=table 'INCO' ;
  177.   rv.'INCO'.'UN' = kcht $bell VECT SOMMET (0 0) ;
  178.   rv.pression=rvp ;
  179.  
  180.  lh= (noeu 10) et (noeu 20) et (noeu 30) et (noeu 40)  et
  181.  (noeu 50) et (noeu 60) ;
  182.  lj= (manu poi1 ((doma $bell 'MAILLAGE') poin proc( 0.5 0.5) ) ) ;
  183.  his = khis 'UN' 1 lh 'UN' 2 (lh et lj) ;
  184.    rv.'HIST'=his ;
  185.  
  186.     exec rv ;
  187.  
  188.  pn=elno (kcht $bell scal centre (rvp.pression)) $bell ;
  189.  ung1 = vect 0.5 (rv.'INCO'.'UN') ux uy jaune ;
  190.  
  191.  qe=dbit (rv.'INCO'.'UN') $entree  ;
  192.  qs=dbit (rv.'INCO'.'UN') $sortie  ;
  193.  dq=(abs qs )-(abs qe) ;
  194.  mess (' Bilan : dq=') dq ;
  195.  pet=somt (redu pn entree) / (nbno entree) ;
  196.  pst=somt (redu pn sortie) / (nbno sortie) ;
  197.  dpes=pet - pst ;
  198.  mess ' pet pst dpes = ' pet pst dpes ;
  199.  
  200.  si ('EGA' graph 'O' );
  201.  dessin (his.'TABD') (his.'1UN') ;
  202.  dessin (his.'TABD') (his.'2UN') ;
  203.  
  204.  trace ung1 bell ;
  205.  trace pn   bell (cont bell) ;
  206.  finsi ;
  207.  
  208. FINPROC RV dq qe dpes ;
  209.  
  210. type=mot tri6 ;
  211. * mess 'Type d élément TRI6 ou QUA8 ? ' ;
  212. * obtenir TYPE*mot  ;
  213.  
  214. angle=0. ;
  215. * obtenir angle*flottant ;
  216.  
  217. nbe=7 ; nbv=10 ;
  218.  
  219.  
  220. type=mot qua8 ;
  221. rv dq qe dpes = hy1 TYPE nbe nbv angle nbit4 tfin graph ;
  222.  erq =  abs (qe - qe14) ;
  223.  mess ' Erreur sur le debit ' erq ;
  224.  si ( erq > err1 ) ; erreur 5 ; finsi ;
  225.  dq=abs dq ;
  226.  si ( dq  > dq1  ) ; erreur 5 ; finsi ;
  227. * On teste quand meme le delta P qui viendrait d'une erreur de discretisation
  228.  ddp=abs ( dpes - dpes1) ;
  229.  mess ' ddp erpe1 '  ddp erpe1 ;
  230.  si ( ddp  > erpe1  ) ; erreur 5 ; finsi ;
  231.  
  232.  
  233. *opti donn 5 ;
  234.  
  235. type=mot tri6 ;
  236. rv dq qe dpes = hy1 TYPE nbe nbv angle nbit3 tfin graph ;
  237.  erq =  abs (qe - qe13) ;
  238.  mess ' Erreur sur le debit ' erq ;
  239.  si ( erq > err1 ) ; erreur 5 ; finsi ;
  240.  dq=abs dq ;
  241.  si ( dq  > dq1  ) ; erreur 5 ; finsi ;
  242.  
  243. FIN ;
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales