* 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=resu( abs (lPR - lpr1 ) ) /100. ;
ERPR2=resu( 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 ;
© Cast3M 2003 - All rights reserved.
Disclaimer