Exemples de jeux de données Cast3M

* fichier :  colline_expl.dgibi
**---------------------------------
** ECOULEMENT AUTOUR D'UNE COLLINE
** G. TURBELIN 14/12/99
** modifie le 6/02/2000    
** modifie le 15/06/2014 passage EQPR -> EQEX
**     
** 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' ;
 
DISCR  = 'MACRO' ;
KPRESS = 'CENTRE';
BETAP  = 1.      ;
 
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
*************************
**------------------------
 
 Msurfa  = 'CHAN' surfa  'QUAF' ;
 Mbas    = 'CHAN' bas    'QUAF' ;
 Mentree = 'CHAN' entree 'QUAF' ;
 Mhaut   = 'CHAN' haut   'QUAF' ;
 Msortie = 'CHAN' sortie 'QUAF' ;
 
 'ELIM' (Msurfa et Mbas et Mentree et Mhaut et Msortie)1.e-5;
 
 $surfa  = 'MODE' Msurfa  'NAVIER_STOKES' DISCR ;
 $bas    = 'MODE' Mbas    'NAVIER_STOKES' DISCR ;
 $entree = 'MODE' Mentree 'NAVIER_STOKES' DISCR ;
 $haut   = 'MODE' Mhaut   'NAVIER_STOKES' DISCR ;
 $sortie = 'MODE' Msortie 'NAVIER_STOKES' DISCR ;
 
 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 = EQEX 'OPTI' 'EF' KPRESS
    'ZONE' $surfa OPER  KBBT  -1. betap 'INCO' 'UN' 'PRES'
    'ZONE' $bas   OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
    'ZONE' $haut  OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
    ;
 
 RV.'PROJ' = RVP     ;
 
 RV.INCO = TABLE INCO;
 RV.INCO.'PRES' = kcht $surfa SCAL KPRESS 0.;
 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' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
 
 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' ;