* fichier : colline_expl.dgibi
**---------------------------------
** ECOULEMENT AUTOUR D'UNE COLLINE
** G. TURBELIN 14/12/99
** modifie le 6/02/2000
**
** CAS TEST DE ALMEIDA et Al.
** INFO DOMAINE & RESULTATS EXPERIMENTAUX A L_ADRESSE
** http://vortex.mech.surrey.ac.uk/database/test18/test18.html
**----------------------------------
opti dime 2 elem qua8;
COMPLET=FAUX;
GRAPH='N' ;
Si complet ;
iitma=20000 ;
sinon;
iitma=42 ;
finsi ;
**************************
** PROCEDURE DE CALCUL
** POUR ESTIMER LA CONVERGENCE
****************************************
'DEBPROC' ERR;
'ARGUMENT' rvx*'TABLE';
rv= rvx.'EQEX';
DD = rv.PASDETPS.'NUPASDT' ;
NN = DD/20;
T0 = (DD '-' (20*NN)) 'EGA' 0;
'SI' T0;
UN = RV.INCO.'UN' ;
UNM1 = RV.INCO.'UNM1' ;
unx= kcht $surfa scal sommet (exco 'UX' un) ;
unm1x = kcht $surfa scal sommet (exco 'UX' unm1) ;
uny= kcht $surfa scal sommet (exco 'UY' un) ;
unm1y = kcht $surfa 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.) ;
***! pour afficher NUT/NU, attention aux dimensions
MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' ' ERREUR LINF ' ELIX ELIY
'MAX NUT/NU = ' ((MAXI RV.INCO.'NUT')*Re) ;
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' ;
RV.INCO.'UNM1'= KCHT $surfa vect sommet (RV.INCO.'UN') ;
as2 ama1 = 'KOPS' 'MATRIK' ;
'FINPROC' as2 ama1 ;
**--------------------------------------------------
************************************
** CONSTRUCTION DU MAILLAGE (ADIM)
*************************************
**
**Definition des longueurs
**-----------------------
**hauteur du domaine (mm)
Hdo = 170.;
**hauteur de la colline (mm)
H_ref = 28.;
**abscisses caracteristiques (mm)
Xmind = -326.; Xmaxd = 500.; X1d=0.;
X2d=9.; X3d= 14.; X4d= 20.;
X5d=30.; X6d=40.; X7d=54.;
**abscisses caracteristiques 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)
si complet ;
Nbcol = 2;
sinon ;
Nbcol = 1;
finsi ;
**Nb de maille entre x=-326 et x =-54 (bas1 & haut4)
Nbbas1 = 13*Nbcol;
****Nb de maille entre x=54 et x =500 (bas2 & haut1)
Nbbas2 = 19*Nbcol;
**Nb de maille entre y=0 et y=170 (sortie et entree)
Nbsort = 10*Nbcol ;
**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 troisieme 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 (1*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 (1*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 (1*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 (1*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 (1*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 (1*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 (1*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 (1*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 superieure doit correspondre avec le
** maillage de la face inferieure. Pour ne pas avoir de pb dans le
** post-traitement (definition 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 (1*Nbcol) (X6 Hd);
haut2D = (X6 Hd) d (1*Nbcol) (X5 Hd);
haut3D = (X5 Hd) d (1*Nbcol) (X4 Hd);
haut4D = (X4 Hd) d (1*Nbcol) (X3 Hd);
haut5D = (X3 Hd) d (1*Nbcol) (X2 Hd);
haut6D = (X2 Hd) d (1*Nbcol) (X1 Hd);
**Cote Gauche
haut6G = ((-1*X1) Hd) d (1*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 (1*Nbcol) ((-1*X5) Hd);
haut2G = ((-1*X5) Hd) d (1*Nbcol) ((-1*X6) Hd);
haut1G = ((-1*X6) Hd) d (1*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;
surfa = dall bas sortie haut entree plan;
surfa = 'CHANGER' surfa 'QUAF' ;
bas = 'CHANGER' bas 'QUAF' ;
entree = 'CHANGER' entree 'QUAF' ;
haut = 'CHANGER' haut 'QUAF' ;
sortie = 'CHANGER' sortie 'QUAF' ;
'ELIMINATION' 1.D-6 (bas 'ET' entree 'ET' haut 'ET' sortie
'ET' surfa) ;
** ! Le maillage est sans dimension ...
*************************
** DEFINITION DES MODELES
*************************
**------------------------
$surfa = 'MODE' surfa 'NAVIER_STOKES' MACRO ;
$bas = 'MODE' bas 'NAVIER_STOKES' MACRO ;
$entree = 'MODE' entree 'NAVIER_STOKES' MACRO ;
$haut = 'MODE' haut 'NAVIER_STOKES' MACRO ;
$sortie = 'MODE' sortie 'NAVIER_STOKES' MACRO ;
surfa = 'DOMA' $surfa 'MAILLAGE' ;
'ORIENTER' surfa;
entree = 'DOMA' $entree 'MAILLAGE' ;
haut = 'DOMA' $haut 'MAILLAGE' ;
bas = 'DOMA' $bas 'MAILLAGE' ;
sortie = 'DOMA' $sortie 'MAILLAGE' ;
si ('EGA' graph 'O' );
'TRACER' surfa;
finsi ;
**********************
**VALEURS DE REFERENCE
**(mm)-->(m)
**********************
**Hauteur du canal (m)
Hdo = Hdo*1.e-3;
**1/2 hauteur (axe)
Haxe = Hdo/2.;
**Hauteur de reference (m)
H_ref = H_ref*1.e-3 ;
**vitesse de reference (m/s) = Vx(85)
U_ref = 2.147;
*vitesse de frottement (aerodynamique) en m/s
Uf = 0.079;
*****************************
*** PROPRIETES PHYSIQUES
***************************
** constante de Von Karman
C1 = 0.4;
** constante du modele RNG k-e
CNU = 0.0845;
** taux de turbulence
INT = 0.03 ;
** Longueur de reference pour FILTREKE
Lref = (Hdo/H_ref) ;
** 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 DES POINTS DU MAILLAGE
**************************************
** ORDO des points du maillage (sans dim)
Y1 = coor 2 surfa;
** 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);
** Distances Normales aux parois (1/2 hauteur de la premiere maille) (sans dim)
**-----------------------------------------------
YP = (EXTR (extr YC2 ordo) 2)/2;
** Prise en compte de Yp dans les ordonnees (y=y+yp)
Yyp = MANU CHPO entree 1 scal YP;
Yyp = Y2 + Yyp;
** Coordonnes pour des profils symetriques / a l_axe
Yc = MANU CHPO entree 1 scal ((Haxe/H_ref)+YP);
Ysym = Yc -(abs(Yyp - Yc));
Yslog = log Ysym;
** Definition des profils d'entrees et initiaux
**----------------------------------------------
** Definition du domaine d'imposition
M_imp =entree;
** CALCUL DES PROFILS D'ENTREE
**---------------------------------------------------
** (U= a*lnYyp + b) et (K= c*Yyp + d)
** (a,b,c,d) sont obtenus a partir des profils experimentaux.
a = 0.25756 ;
b = 1.8552 ;
c = -1.85184878e-1*H_ref;
d = 1.912227164e-2;
d = MANU CHPO entree 1 scal d;
** VITESSES
**-----------
Uent = (Yslog*a)+b ;
Uent = Uent /U_ref;
** Transformation du chpo scalaire en
** chpo vecteur + imposition sur le domaine
U1 = Uent NOMC 'UX' ;
U1 = KCHT $surfa 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 surfa 1 scal 1.;
U_ini = U_ini NOMC 'UX' ;
U_ini = KCHT $surfa vect sommet U_ini;
*opti donn 5;
**ENERGIE CINETIQUE
**-----------------
** Kadim = K/(U_ref**2)
*
Kent = ((Ysym*c)+ d);
Kent = (Kent/(U_ref**2));
** Fonction de l'intensite de la turbulence
K2 = 1.5*(INT**2) ;
K2 = 'MANU' 'CHPO' surfa 1 scal K2;
***INITIALE
K_ini = KCHT $surfa scal sommet K2;
***IMPOSEE
K_IMP = KCHT $surfa 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 $surfa scal sommet E2;
** IMPOSEE
E_IMP= KCHT $surfa scal sommet E2 E3 ;
E_IMP=EXCO (REDU E_IMP M_imp) SCAL ;
*list E_imp;
** 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 $surfa scal sommet (alf1/Re) ;
N_ini = noel $surfa 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 $surfa ITMA iitma ALFA 0.8
OPTI 'SUPG' 'RNG'
ZONE $surfa 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 $surfa OPER FILTREKE 1. Lref iRe 'UN' INCO 'KN' 'EN';
RV = EQEX RV 'OPTI' 'CENTREE'
'ZONE' $surfa 'OPER' 'DFDT' 1. 'UN' 'DELTAT' INCO 'UN'
'ZONE' $surfa 'OPER' 'DFDT' 1. 'KN' 'DELTAT' INCO 'KN'
'ZONE' $surfa 'OPER' 'DFDT' 1. 'EN' 'DELTAT' INCO 'EN' ;
RV = 'EQEX' RV
ZONE $surfa OPER ERR;
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' $surfa
ZONE $surfa 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 $surfa vect sommet (1.E-3 1.E-3) ;
RV.INCO.'IT' = PROG 1 ;
RV.INCO.'ERY' = PROG 0. ;
RV.INCO.'ERX' = PROG 0.;
EXEC RV;
rt=rv.pasdetps;
tps1 = rt.tps;
it1 = rt.nupasdt;
dt = rt.'DELTAT-1';
'MESSAGE' 'No ITERATION ' it1;
'MESSAGE' 'DERNIER PAS ' DT;
'MESSAGE' 'TEMPS PHYSIQUE' tps1;
** PROFILS SUR LE DOMAINE TOTAL
**--------------------------------
vv = 'VECTEUR' rv.inco.'UN' .1 ux uy ;
Vx = rv.inco.'UN' exco UX SCAL ;
Vy = rv.inco.'UN' exco UY SCAL ;
kk = rv.inco.'KN' ;
EE = rv.inco.'EN';
nuts = 'ELNO' ('KCHT' $surfa scal centre (RV.INCO.'NUT')) $surfa;
NUTSNU = nuts * Re ;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
si ('EGA' graph 'O' );
*'OPTION' isov ligne;
titre 'Vecteur Vitesse à t=' tps1;
'TRACER' vv surfa ('CONTOUR' surfa);
titre 'vecteur vitesse composante X à t=' tps1;
'TRACER' Vx surfa ('CONTOUR' surfa) ;
titre 'vecteur vitesse composante Y à t=' tps1;
'TRACER' VY surfa ('CONTOUR' surfa) ;
titre 'CHAMP DE PRESSION';
'TRACER' PP surfa ('CONTOUR' surfa);
titre 'Energie Cinetique Turbulente à t=' tps1;
'TRACER' kk surfa ('CONTOUR' surfa);
titre 'dissipation à t=' tps1;
'TRACER' ee surfa ('CONTOUR' surfa);
titre 'NUT/NU à t=' tps1;
'TRACER' nutsnu surfa ('CONTOUR' surfa);
finsi ;
** FONCTION DE COURANT
**---------------------
* FONCTION DE COURANT
**---------------------
un1=EXCO (rv.inco.'UN') ux;
un2=EXCO (rv.inco.'UN') uy;
un1=KCHT $surfa SCAL SOMMET un1;
un2=KCHT $surfa SCAL SOMMET un2;
** CALCUL DU ROTATIONNEL DU CHAMP DE VITESSE
rt2d=KOPS (rv.inco.'UN') 'ROT' $surfa;
sw= (-1*rt2d);
** DEFINITION DE Utangentielle SUR CHAQUE FACE
** !! la direction de Ut depend de la définition
** de la normale a la paroi n !!
unn1=KCHT $BAS SCAL SOMMET (-1*un1);
unn2=KCHT $SORTIE SCAL SOMMET (-1*un2);
unn3=KCHT $HAUT SCAL SOMMET (un1);
unn4=KCHT $ENTREE SCAL SOMMET (un2);
unn1=NOEL $BAS unn1;
unn1=KCHT $BAS SCAL CENTRE unn1;
unn2=NOEL $SORTIE unn2;
unn2=KCHT $SORTIE SCAL CENTRE unn2;
unn3=NOEL $HAUT unn3;
unn3=KCHT $HAUT SCAL CENTRE unn3;
unn4=NOEL $ENTREE unn4;
unn4=KCHT $ENTREE SCAL CENTRE unn4;
** RESOLUTION DE LAPN(PSI)=ROT(UN)
** avec D(PSI)/D(n) = (+/-)Utangentielle AUX LIMITES DU DOMAINE, (ce qui
** se traduit par des conditions de type flux imposés (FIMP) aux frontieres.
** ON PREND PSI(BAS) = 0
** ICI, SEULE LA CONDITION SUR HAUT EST NECESSAIRE
rk=EQEX $surfa 'OPTI' 'EF' 'IMPL'
ZONE $surfa OPER LAPN 1. inco 'psi'
ZONE $surfa OPER FIMP sw inco 'psi'
ZONE $BAS OPER FIMP unn1 inco 'psi'
ZONE $SORTIE OPER FIMP unn2 inco 'psi'
ZONE $HAUT OPER FIMP unn3 inco 'psi'
ZONE $ENTREE OPER FIMP unn4 inco 'psi'
CLIM 'psi' TIMP BAS 0. ;
rk.inco.'psi'=KCHT $surfa SCAL SOMMET 0.;
exec rk;
** Lignes de courant = ISO de la fonction de courant
si ('EGA' graph 'O' );
psi=rk.inco.'psi';
'OPTION' ISOV LIGNE;
'TRACER' psi surfa ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
finsi ;
*opti donn 5 ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer