Accueil | Cast3M

* fichier :  hy4.dgibi
************************************************************************
************************************************************************
*$$$ HY4
 
 
**  Exemple HY4
**  ___________
**
**  --- 2 FEVRIER 2002 ---
**
**   CANAL LONGUEUR  3 X 10. LARGEUR 1.
**   test cas isotherme NS et FROT (Faisceau de tube)
**   Sur le tronçon du milieu on impose une forte perte de charge
**   qui a pour effet de redresser le champ de vitesse
**   Les conditions limites (adhérence en paroi) conviennent pour
**   la zone Poiseuille mais physiquement ne sont pas adaptée pour
**   le tronçon faisceau de tube. Il faudrait imposer un glissement.
**   Mais ça marche !
 
GRAPH=VRAI ;
GRAPH=FAUX ;
DISCR=QUAF ;
DISCR=LINE ;
KPRES=MSOMMET ;
 
type=mot tri6 ;
type=mot QUA8 ;
*mess 'Type d élément TRI6 ou QUA8 ? ' ;
*obtenir TYPE*mot  ;
 
OPTION DIME 2 ELEM TYPE ;
nbe=7 ; nbv=10 ;
 
p1=0  0.;
p2=1. 0.;
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 (0 10)  ;
q2=p2 plus (0 10)  ;
pp1=p1 d q1 nbv;
pp2=p2 d q2 nbv;
sortie=entree plus (0 10) ;
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 ;
*bell=surf cnt ;
bell=daller entree pp2 sortie pp1 ;
angle=0. ;
*obtenir angle*flottant ;
 
bell=bell 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 bell) 0.0001 ;
 
 
* On met bout a bout 3 longueur de troncon
 
bell_2 sortie_2 pp1_2 pp2_2 = bell sortie pp1 pp2 'PLUS' (0. 10.) ;
bell_3 sortie_3 pp1_3 pp2_3 = bell_2 sortie_2 pp1_2 pp2_2
 'PLUS' (0. 10.) ;
 
 bell = bell et bell_2 et bell_3;
 pp1  = pp1  et pp1_2  et pp1_3 ;
 pp2  = pp2  et pp2_2  et pp2_3 ;
 sortie = sortie_3 ;
 
 elim (pp1 et pp2 et entree et sortie et bell) 0.0001 ;
 
 Mbell = chan bell QUAF ;
 MbellM= chan bell_2 QUAF ;
 Mentree = chan entree QUAF ;
 Msortie = chan sortie QUAF ;
 Mpp1    = chan pp1    QUAF ;
 Mpp2    = chan pp2    QUAF ;
 
 Elim (Mbell et Mentree et Msortie et Mpp1 et Mpp2 et MbellM) 1.e-5 ;
 
 $bell=MODE Mbell 'NAVIER_STOKES' DISCR ;
 $bellM=MODE MbellM 'NAVIER_STOKES' DISCR ;
 $entree = mode Mentree 'NAVIER_STOKES' DISCR ;
 $sortie = mode Msortie 'NAVIER_STOKES' DISCR ;
 $pp1    = mode Mpp1    'NAVIER_STOKES' DISCR ;
 $pp2    = mode Mpp2    'NAVIER_STOKES' DISCR ;
 PP1 = DOMA $PP1 MAILLAGE ;
 PP2 = DOMA $PP2 MAILLAGE ;
 pp1p = DOMA $PP1 MMAIL ;
 pp2p = DOMA $PP2 MMAIL ;
 ENTREE = DOMA $ENTREE MAILLAGE ;
 SORTIE = DOMA $SORTIE MAILLAGE ;
 SORTIP = DOMA $SORTIE MSOMMET  ;
 
 
  mu=5.E-2 ;
  ro=1 ;
  nu=mu/ro ;
 
  ck=  1.E+5 1.E+2;
  cb= 2. 2. ;
 
  DT = 1. ;
 
  rv=eqex $bell 'DUMP' ALFA 1.   ITMA 10
  OPTI 'EF' 'IMPL'
  ZONE  $BELLM   OPER    FROT  CK CB UN  INCO 'UN'
  ZONE  $BELL    OPER    NS 1. 'UN' NU   INCO 'UN'
  OPTI 'EF' 'IMPL' 'CENTREE'
  ZONE  $BELL   OPER    DFDT 1. 'UNM' DT INCO 'UN'
  CLIM
  UN UIMP (PP1 ET PP2 et entree) 0.
  UN VIMP  entree 1. UN VIMP (PP1 ET PP2) 0.
  ;
 
   RVP= EQEX
    'OPTI' 'EF' KPRES
    ZONE  $bell   OPER  KBBT  (-1.)     INCO 'UN' 'PRES'
   CLIM PRES TIMP sortip 0.
   ;
 
    RV.'PROJ'= RVP ;
 
 
rv.'INCO'=table 'INCO' ;
rv.'INCO'.'UN' = kcht $bell VECT SOMMET (0 1) ;
rv.'INCO'.'UNM' = kcht $bell VECT SOMMET (0 1) ;
rv.inco.'PRES'= kcht $bell scal kpres 0.;
 
 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 ;
 
*dessin (his.'TABD') (his.'1UN') ;
*dessin (his.'TABD') (his.'2UN') ;
 
Si (EGA KPRES 'MSOMMET') ;
mtp= doma $bell'MMAIL' ;
p= (rv.inco.pression);
 
evolp1 = EVOL 'CHPO' p 'PRES' pp1p ;
evolp2 = EVOL 'CHPO' p 'PRES' pp2p ;
 
lpr1=(extr evolp1 'ORDO') ;
lpr2=(extr evolp2 'ORDO') ;
list lpr1 ; list lpr2 ;
 
 Lpr= prog
   110.68       110.06       109.23       108.66       108.08
   107.54       106.83       106.48       105.36       105.84
   103.10       94.036       84.523       75.032       65.533
   56.035       46.538       37.039       27.546       18.038
   8.8027       4.5193       4.5952       4.0910       3.4839
   2.8817       2.2759       1.7034       1.0727       .55958
   0.00000E+00;
 
      ERPR1=SOMM( abs (lPR  - lpr1 ) ) /100. ;
      ERPR2=SOMM( abs (lPR  - lpr2 ) ) /100. ;
 
Si GRAPH ;
dess (evolp1 et evolp2 ) ;
trace p mtp (cont mtp) ;
FINSI ;
 
Sinon;
p=elno $mt (kcht $mt scal kpres
 (exco (rv.inco.pression) 'PRES'))  kpres;
trace p mt TITRE ctit;
finsi ;
 
 
 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 ;
 
Si GRAPH ;
 trace ung1 bell ;
FINSI ;
 
Mess ' ERPR1=' ERPR1 ' ERPR2=' ERPR2 ;
Si (ERPR1 '>' 2.E-2) ; erreur 5 ; Finsi ;
Si (ERPR2 '>' 2.e-2) ; erreur 5 ; Finsi ;
 
FIN ;