* fichier : colline.dgibi
**---------------------------------
** ECOULEMENT AUTOUR D'UNE COLLINE
** G. TURBELIN 14/12/99
**----------------------------------
opti dime 2 elem qua8;
COMPLET=FAUX;
GRAPH='N' ;
Si complet ;
iitma=20000 ;
sinon;
iitma=5 ;
finsi ;
*------------------ PROCEDURE FILTREKE ---------------------------------
DEBP FILTREKE ;
ARGU RX*TABLE ;
* Filtre sur K et Epsilon
* - Echelle de vitesse (K**0.5) inférieure à une fraction (alfk)
* de Uref (vitesse caractéristique) (alfk=1 pour l'instant)
* Uref=max(UN,U0)
* - K > K0
* - Epsilon tel que l'echelle de longueur reste inférieure
* à (L0/a) où L0 = diamètre enceinte et a=f(Re)
* => Nut < Uref*L0/a
rv=rx.'EQEX' ;
rvp=rv.'PRESSION' ;
iarg=rx.'IARG' ;
NASTOK = rv.'NAVISTOK' ;
si( non ( ega iarg 4)) ;
mess 'Procedure FILTREKE : nombre d arguments incorrect ' iarg ;
quitter FILTREKE ;
finsi ;
si ( ega ('TYPE' rx.'ARG1') 'MOT ') ;
U1=rv.'INCO'.(rx.'ARG1') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG1')) 'FLOTTANT') ;
U1=rx.'ARG1' ;
sinon ;
mess 'Procedure FILTREKE : type argument 1 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
si ( ega ('TYPE' rx.'ARG2') 'MOT ') ;
L0=rv.'INCO'.(rx.'ARG2') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG2')) 'FLOTTANT') ;
L0=rx.'ARG2' ;
sinon ;
mess 'Procedure FILTREKE : type argument 2 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
si ( ega ('TYPE' rx.'ARG3') 'MOT ') ;
NU=rv.'INCO'.(rx.'ARG3') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG3')) 'FLOTTANT') ;
NU=rx.'ARG3' ;
sinon ;
mess 'Procedure FILTREKE : type argument 3 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
si ( ega ('TYPE' rx.'ARG4') 'MOT ') ;
UN=rv.'INCO'.(rx.'ARG4') ;
sinon ;
si ( ega ('TYPE' (rx.'ARG4')) 'CHPOINT') ;
UN=rx.'ARG4' ;
sinon ;
mess 'Procedure FILTREKE : type argument 4 invalide ' ;
quitter FILTREKE ;
finsi ;
finsi ;
nic=dime (rx.'LISTINCO') ;
si( non ( ega nic 2)) ;
mess 'Procedure FILTREKE : nombre d inconnues incorrect ' nic ;
quitter FILTREKE ;
finsi ;
nomi1=extr 1 (rx.'LISTINCO');
nomi2=extr 2 (rx.'LISTINCO');
nom1= mot (text (chai nomi1)) ;
nom2= mot (text (chai nomi2)) ;
en=rv.'INCO'.nom2 ;
kn=rv.'INCO'.nom1 ;
Rec=100.;
k0 = 1.e-10 ;
cnu=0.09;
mdu=kops un 'PSCA' un ;
mdu=mdu ** 0.5 ;
Re=kops (kops (kops mdu '*' L0) '/' nu) '+' (Rec / 10.) ;
a= exp (kops Rec '/' Re ) ;
mdu = kops mdu '|<' u1 ;
mdu2= kops mdu '*' mdu ;
kn=kops kn '|<' k0 ;
kn=kops kn '>|' mdu2 ;
E0= kops (kops kn '**' 1.5) '*' (a / L0) ;
en=kops en '|<' E0 ;
rv.'INCO'.nom2=en ;
rv.'INCO'.nom1=kn ;
as2 ama1 = 'KOPS' 'MATRIK' ;
*RESPRO as2 ama1 ;
FINPROC ;
*------------------------------- FIN PROCEDURE FILTREKE ----------------
** CONSTRUCTION DU MAILLAGE
**------------------------
**Définition des longueurs
**-----------------------
**hauteur du domaine (mm)
Hdo = 170.;
**hauteur de la colline (mm)
H_ref = 28.;
**abscisses caractéristiques (mm)
Xmind = -326.; Xmaxd = 500.; X1d=0.;
X2d=9.; X3d= 14.; X4d= 20.;
X5d=30.; X6d=40.; X7d=54.;
**abscisses caractéristiques adim
Xmin = Xmind/H_ref ; Xmax = Xmaxd/H_ref; X1=X1d/H_ref;
X2=X2d/H_ref; X3= X3d/H_ref; X4=X4d/H_ref;
X5=X5d/H_ref; X6=X6d/H_ref; X7=X7d/H_ref;
Hd = Hdo/H_ref ;
**Nb de maille
**Nombre de maille entre x=9 et x=14(colline)
Nbcol = 1;
**Nb de maille entre x=-326 et x =-54 (bas1 & haut4)
Nbbas1 = 25;
****Nb de maille entre x=54 et x =500 (bas2 & haut1)
Nbbas2 = 40;
**Nb de maille entre y=0 et x =170 (sortie et entree)
**
Nbsort = 20 ;
**Densité
DIbasG = -1*(Xmin+X7)/Nbbas1 ;
DIbas1 = 2*DIbasG ; DFbas1 = 0.5*DIbasG;
DIbasD = (Xmax-X7)/Nbbas1 ;
DIbas2 = 0.8*DIbasD ; DFbas2 = 1.5*DIbasD;
Dsort = Hd '/' Nbsort ;
DIsort = 0.3* Dsort;
DFsort = 2.* Dsort;
**-----------------------------------------
** Construction de la colline a partir
** de polynomes du troisième ordre
**-----------------------------------------
**
** Cote GAUCHE
coG60 = 0. (56.39011190988/H_ref);
coG61 = 1. (+2.010520359035);
coG62 = 0. ((1.644919857549e-2)*H_ref) ;
coG63 = 0. (-2.674976141766e-5 *(H_ref**2));
lignG6 = COURBE (3*Nbcol) 'DINI'1. 'DFIN'1.
coG60 coG61 coG62 coG63 'PARAMETRE' (-1*X7) (-1*X6) ;
coG50 = 0. (17.92461334664/H_ref);
coG51 = 1. -8.743920332081e-1 ;
coG52 = 0. ((-5.567361123058e-2)*H_ref) ;
coG53 = 0. ((-6.277731764683e-4)*(H_ref**2));
lignG5 = COURBE (2*Nbcol) 'DINI'1. 'DFIN'1.
coG50 coG51 coG52 coG53 'PARAMETRE' (-1*X6) (-1*X5) ;
coG40 = 0. (40.46435022819/H_ref);
coG41 = 1. +1.379581654948;
coG42 = 0. ((1.945884504128e-2)*H_ref);
coG43 = 0. ((+2.070318932190e-4)*(H_ref**2));
lignG4 = COURBE (2*Nbcol) 'DINI'1. 'DFIN'1.
coG40 coG41 coG42 coG43 'PARAMETRE' (-1*X5) (-1*X4) ;
coG30 = 0. (25.79601052357/H_ref);
coG31 = 1. -8.20669300745e-1;
coG32 = 0. ((-9.055370274339e-2)*H_ref) ;
coG33 = 0. ((-1.626510569859e-3)*(H_ref**2));
lignG3 = COURBE (1*Nbcol) 'DINI'1. 'DFIN'1.
coG30 coG31 coG32 coG33 'PARAMETRE' (-1*X4) (-1*X3) ;
coG20 = 0. ((2.507355893131e+1)/H_ref);
coG21 = 1. -9.754803562315e-1;
coG22 = 0. ((-1.016116352781e-1)*H_ref) ;
coG23 = 0. ((-1.889794677828e-3)*(H_ref**2)) ;
lignG2 = COURBE (1*Nbcol) 'DINI'1. 'DFIN'1.
coG20 coG21 coG22 coG23 'PARAMETRE' (-1*X3) (-1*X2) ;
coG10 = 0. (28/H_ref);
coG11 = 1. 0.;
coG12 = 0. ((6.775070969851e-3)*H_ref) ;
coG13 = 0. ((+2.124527775800e-3)*(H_ref**2)) ;
lignG1 = COURBE (2*Nbcol) 'DINI'1. 'DFIN'1.
coG10 coG11 coG12 coG13 'PARAMETRE' (-1*X2) X1 ;
colG = lignG6 et lignG5 et lignG4 et lignG3
et lignG2 et lignG1 ;
elim colG 1e-5;
*trac colG;
**Cote Droit
coD10 = 0. (28/H_ref);
coD11 = 1. 0.;
coD12 = 0. ((6.775070969851e-3)*H_ref) ;
coD13 = 0. ((-2.124527775800e-3)*(H_ref**2));
lignD1 = COURBE (2*Nbcol) 'DINI'1. 'DFIN'1. coD10
coD11 coD12 coD13 'PARAMETRE' X1 X2 ;
coD20 = 0. ((2.507355893131e+1)/H_ref);
coD21 = 1. 9.754803562315e-1;
coD22 = 0. ((-1.016116352781e-1)*H_ref) ;
coD23 = 0. ((1.889794677828e-3)*(H_ref**2)) ;
lignD2 = COURBE (1*Nbcol) 'DINI'1. 'DFIN'1. coD20 coD21
coD22 coD23 'PARAMETRE' X2 X3 ;
coD30 = 0. ((25.79601052357)/H_ref);
coD31 = 1. 8.20669300745e-1;
coD32 = 0. ((-9.055370274339e-2)*H_ref) ;
coD33 = 0. ((1.626510569859e-3)*(H_ref**2)) ;
lignD3 = COURBE (1*Nbcol) 'DINI'1. 'DFIN'1. coD30 coD31
coD32 coD33 'PARAMETRE' X3 X4 ;
coD40 = 0. (40.46435022819/H_ref);
coD41 = 1. -1.379581654948;
coD42 = 0. ((1.945884504128e-2)*H_ref);
coD43 = 0. ((-2.070318932190e-4)*(H_ref**2));
lignD4 = COURBE (2*Nbcol) 'DINI'1. 'DFIN'1. coD40
coD41 coD42 coD43 'PARAMETRE' X4 X5 ;
coD50 = 0. (17.92461334664/H_ref);
coD51 = 1. 8.743920332081e-1;
coD52 = 0. ((-5.567361123058e-2)*H_ref) ;
coD53 = 0. ((6.277731764683e-4)*(H_ref**2));
lignD5 = COURBE (2*Nbcol) 'DINI'1. 'DFIN'1. coD50 coD51
coD52 coD53 'PARAMETRE' X5 X6 ;
coD60 = 0. (56.39011190988/H_ref);
coD61 = 1. -2.010520359035;
coD62 = 0. ((1.644919857549e-2)*H_ref) ;
coD63 = 0. ((2.674976141766e-5)*(H_ref**2)) ;
lignD6 = COURBE (3*Nbcol) 'DINI'1. 'DFIN'1. coD60
coD61 coD62 coD63 'PARAMETRE' X6 X7 ;
colD = lignD1 et lignD2 et lignD3 et lignD4
et lignD5 et lignD6 ;
elim colD 1e-6;
colline = colG et colD;
elim colline 1e-6;
*trac colline;opti donn 5 ;
**
**Construction des cotes du domaine
**-----------------------------------
bas1 = (Xmin 0.) d (-1*Nbbas1) ((-1*X7) 0.) DINI DIbas1 DFIN DFbas1;
bas2 = (X7 0.) d (-1*Nbbas2) (Xmax 0.) DINI DIbas2 DFIN DFbas2;
bas = bas1 et colline et bas2;
elim bas 1e-6;
** Le maillage de la face supérieure doit correspondre avec le
** maillage de la face inférieure. Pour ne pas avoir de pb dans le
** post-traitement (définition de lignes a X=Cst) , il est utile d'avoir
** un maillage rectiligne.
*Haut Droite et Haut Gauche
hautD = (Xmax Hd) d (-1*Nbbas2) (X7 Hd) DINI DFbas2 DFIN DIbas2;
hautG = ((-1*X7) Hd) d (-1*Nbbas1) (Xmin Hd) DINI DFbas1 DFIN DIbas1;
**cote Droit
haut1D = (X7 Hd) d (3*Nbcol) (X6 Hd);
haut2D = (X6 Hd) d (2*Nbcol) (X5 Hd);
haut3D = (X5 Hd) d (2*Nbcol) (X4 Hd);
haut4D = (X4 Hd) d (1*Nbcol) (X3 Hd);
haut5D = (X3 Hd) d (1*Nbcol) (X2 Hd);
haut6D = (X2 Hd) d (2*Nbcol) (X1 Hd);
**Cote Gauche
haut6G = ((-1*X1) Hd) d (2*Nbcol) ((-1*X2) Hd);
haut5G = ((-1*X2) Hd) d (1*Nbcol) ((-1*X3) Hd);
haut4G = ((-1*X3) Hd) d (1*Nbcol) ((-1*X4) Hd);
haut3G = ((-1*X4) Hd) d (2*Nbcol) ((-1*X5) Hd);
haut2G = ((-1*X5) Hd) d (2*Nbcol) ((-1*X6) Hd);
haut1G = ((-1*X6) Hd) d (3*Nbcol) ((-1*X7) Hd);
haut = hautD et haut1D et haut2D et haut3D et haut4D et
haut5D et haut6D et haut6G et haut5G et haut4G et
haut3G et haut2G et haut1G et hautG;
elim haut 1e-6;
sortie1 =(Xmax 0.) d (-1*(Nbsort/2)) (Xmax (Hd/2))
DINI DIsort DFIN DFsort;
sortie2 =(Xmax (Hd/2)) d (-1*(Nbsort/2)) (Xmax Hd)
DINI DFsort DFIN DIsort;
sortie = (sortie1 et sortie2);
elim sortie 1e-5;
entree1 =(Xmin Hd) d (-1*(Nbsort/2)) (Xmin (Hd/2))
DINI DIsort DFIN DFsort;
entree2 =(Xmin (Hd/2)) d (-1*(Nbsort/2)) (Xmin 0.)
DINI DFsort DFIN DIsort;
entree = (entree1 et entree2);
elim entree 1e-5;
elim (bas et sortie et haut et entree ) 1e-5;
surf2 = dall bas sortie haut entree plan;
ORIENTER surf2;
*trac surf2; opti donn 5;
** Le maillage est sans dimension ...
**-------------------------------
**CALCUL DE L'ERREUR
** ------------------------------
'DEBPROC' ERR;
'ARGUMENT' rvx*'TABLE';
rv= rvx.'EQEX';
DD = rv.PASDETPS.'NUPASDT' ;
NN = DD/5;
L0 = (DD '-' (5*NN)) 'EGA' 0;
'SI' L0;
UN = RV.INCO.'UN' ;
UNM1 = RV.INCO.'UNM1' ;
unx= kcht (rv.'DOMAINE') scal sommet (exco 'UX' un) ;
unm1x = kcht (rv.'DOMAINE') scal sommet (exco 'UX' unm1) ;
uny= kcht (rv.'DOMAINE') scal sommet (exco 'UY' un) ;
unm1y = kcht (rv.'DOMAINE') scal sommet (exco 'UY' unm1) ;
ERRX = KOPS unx '-' unm1x ;
ERRY = KOPS uny '-' unm1y ;
ELIX = MAXI ERRX 'ABS' ;
ELIY = MAXI ERRY 'ABS' ;
ELIX = (LOG (ELIX + 1.0E-20))/(LOG 20.) ;
ELIY = (LOG (ELIY + 1.0E-20))/(LOG 20.) ;
MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' ' ERREUR LINF ' ELIX ELIY
'MAX NUT/NU = ' ((MAXI RV.INCO.'NUT')/NU) ;
RV.INCO.'UNM1'= KCHT (rv.'DOMAINE') vect sommet (RV.INCO.'UN') ;
IT = PROG RV.PASDETPS.'NUPASDT' ;
ERY = PROG ELIY ;
ERX = PROG ELIX ;
RV.INCO.'IT' = (RV.INCO.'IT') ET IT ;
RV.INCO.'ERY' = (RV.INCO.'ERY') ET ERY ;
RV.INCO.'ERX' = (RV.INCO.'ERX') ET ERX ;
'FINSI' ;
'FINPROC' ;
**------------------------------------------------
**
** Définition des domaines
**------------------------
$surf = DOMA surf2 1.e-5 MACRO;
$bas = DOMA bas INCL $surf 1.e-5 MACRO;
$entree = DOMA entree INCL $surf 1.e-5 MACRO;
$haut = DOMA haut INCL $surf 1.e-5 MACRO;
$sortie = DOMA sortie INCL $surf 1.e-5 MACRO;
surf=$surf.maillage;
ORIENTER surf;
entree= $entree.maillage;
haut = $haut.maillage ;
bas = $bas.maillage ;
sortie= $sortie.maillage ;
** CONSTANTES (avec Dim)
**----------------------
**Hauteur du canal (m)
Hc = 170e-3;
**Hauteur de reference (m)
H_ref = 28.e-3 ;
*vitesse de reference (m/s) = Vx(85)
U_ref = 2.147;
* constante de Von Karman
C1 = 0.4;
* constante du modele k-e
CNU = 0.09;
* taux de turbulence
INT = 0.03 ;
*vitesse de frottement (aerodynamique) en m/s
Uf = 0.079;
Lref = (Hc/H_ref) ;
** Proprietes physiques
**-------------------
** nu de l'eau (m**2/s)
nu_eau = 1e-6;
nu = nu_eau;
** NUT_ent/Nu
alf1 = 40.;
**Reynolds de l'écoulement
Re = (U_ref*H_ref)/nu ;
iRe = 1/Re ;
** Coordonnees
**------------
** Ordonnées de tous les points du maillage (sans dim)
Y1 = coor 2 surf;
** Ordonnées de tous les points de l'entree (CHPO et LISTREEL) (sans dim)
Y2 = coor 2 entree;
YC2 = EVOL CHPO Y2 scal (inve entree) ;
** Nombre de noeud en hauteur
NH = DIME (extr YC2 ordo);
*list NH ;
** DISTANCE A LA PAROI (1/2 hauteur de la premiere maille) (sans dim)
**-----------------------------------------------
YP = (EXTR (extr YC2 ordo) 2)/2;
*list Yp;
*opti donn 5;
** Prise en compte de Yp dans les ordonnees (y=y+yp)
Yyp = MANU CHPO entree 1 scal YP;
Yyp = Y2 + Yyp;
** CALCUL DES PROFILS D'ENTREE (un peu complique ...)
**---------------------------------------------------
** (U= a*lnYyp + b) et (K= c*Yyp + d)
** (a,b,c,d) sont obtenus a partir des profils experimentaux.
** Definition du domaine ou les profils seront imposes
M_imp =entree;
** Calcul du profil symetrique
Yc = MANU CHPO entree 1 scal ((85.e-3/H_ref)+YP);
Ysym = Yc -(abs(Yyp - Yc));
Yslog = log Ysym;
** Profil de vitesse (sans dim) : Calcul des coeff a et b
**------------------
** De maniere grossiere
a = 0.873/( (log((85.e-3/H_ref)+YP)) - (log((1.e-3/ H_ref)+YP)) ) ;
*list a;
b= 2.147-( a*(log((85.e-3/H_ref)+YP)) );
*list b ;
b=MANU CHPO entree 1 scal b;
*opti donn 5;
Uent = (Yslog*a)+b ;
Uent = Uent /U_ref;
*List Uent ;
*opti donn 5;
** Transformation du chpo scalaire en
** chpo vecteur + imposition sur le domaine
U1 = Uent NOMC 'UX' ;
U1 = KCHT $surf vect sommet U1;
** Valeurs IMPOSEES
UX_IMP = EXCO (REDU U1 M_imp) 'UX' SCAL;
UY_IMP = EXCO (REDU U1 M_imp) 'UY' SCAL;
** Valeurs INITIALES (sans dim)
U_ini = MANU CHPO surf 1 scal 1.;
U_ini = U_ini NOMC 'UX' ;
U_ini = KCHT $surf vect sommet U_ini;
*opti donn 5;
**ENERGIE CINETIQUE
**-----------------
** Kadim = K/(U_ref**2)
*
** Profil ECT : coeff c et d
**------------------
c = -1.85184878e-1*H_ref;
d = 1.912227164e-2;
d = MANU CHPO entree 1 scal d;
Kent = ((Ysym*c)+ d);
Kent = (Kent/(U_ref**2));
** Fonction de l'intensite de la turbulence
K2 = (INT)**2 ;
K2 = 'MANU' 'CHPO' surf 1 scal K2;
***INITIALE
K_ini = KCHT $surf scal sommet K2;
***IMPOSEE
K_IMP = KCHT $surf scal sommet K2 Kent ;
K_IMP = EXCO (REDU K_IMP M_imp) SCAL;
** DISSIPATION
**--------------
** Fonction de K_ini et de NUT
** K_ini est sans dim
E2 = (CNU*(K_ini**2)*Re)/(alf1);
** En entree, fonction de Kent et d'une longueur caractéristique
** K_ent est sans dim
E3 = ((Kent**1.5)*H_ref)/(Lref);
** INITIALE
E_ini = KCHT $surf scal sommet E2;
** IMPOSEE
E_IMP= KCHT $surf scal sommet E2 E3 ;
E_IMP=EXCO (REDU E_IMP M_imp) SCAL ;
*list E_imp;
*opti donn 5;
** NUT INITIALE
** ------------------
** NUTadim = NUT/(Uref*Href)
** NUadim =1/Re
** Fonction de K et E
*N2=CNU*(K_ini**2)*(E_ini**-1);
** Fonction de nu
N3=KCHT $surf scal sommet (alf1/Re) ;
N_ini = noel $surf N3 ;
** U* INITIALES
**-------------
UET_ini1 = kcht $bas scal centre (Uf/U_ref);
UET_ini2 = kcht $haut scal centre (Uf/U_ref);
*opti donn 5;
** Module de resolution
**-------------------
RV = EQEX $surf ITMA iitma ALFA 0.8
OPTI 'SUPG' 'RNG'
ZONE $surf OPER ERR
ZONE $surf OPER NSKE iRe 'NUT' INCO 'UN' 'KN' 'EN'
ZONE $bas OPER FPU iRe 'UET1' YP INCO 'UN' 'KN' 'EN'
ZONE $haut OPER FPU iRe 'UET2' YP INCO 'UN' 'KN' 'EN'
ZONE $surf OPER FILTREKE 1. Lref iRe 'UN' INCO 'KN' 'EN';
RV = EQEX RV
CLIM 'UN' UIMP M_imp UX_imp
CLIM 'UN' VIMP M_imp UY_imp
CLIM 'KN' TIMP M_imp K_imp
CLIM 'EN' TIMP M_imp E_imp ;
RVP = eqpr $surf
ZONE $surf OPER PRESSION
ZONE $bas OPER VNIMP 0.
ZONE $haut OPER VNIMP 0.;
RV.PRESSION = RVP;
RV.INCO = TABLE INCO;
RV.INCO.'UN' = U_ini;
RV.INCO.'KN' = K_ini;
RV.INCO.'EN' = E_ini;
RV.INCO.'NUT' = N_ini;
RV.INCO.'UET1' = UET_ini1;
RV.INCO.'UET2' = UET_ini2;
RV.INCO.'UNM1'= KCHT (rv.'DOMAINE') vect sommet (1.E-3 1.E-3) ;
RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERY' = PROG 0. ;
RV.INCO.'ERX' = PROG 0.;
rv.inco.'TABFK' = TABLE 'FiltreK' ;
rv.inco.'TABFE' = TABLE 'FiltreE' ;
rv.inco.'TTF' = TABLE 'Temps' ;
rv.inco.'iTabF' = PROG ;
lh= (POIN SURF PROC (-2. 3.)) et (POIN SURF PROC (2. 0.5))
et (POIN SURF PROC (2. 3.))
et (POIN SURF PROC (10. 3.)) et (POIN SURF PROC (17. 3.));
his=khis 'UN' 1 lh
'UN' 2 lh
'KN' lh
'EN' lh;
rv.hist=his;
exec RV;
Fin;
© Cast3M 2003 - All rights reserved.
Disclaimer