* fichier : carre.dgibi
*****************************************************
** ECOULEMENT AUTOUR D'UNE CYLINDRE DE SECTION CAREE*
** Gregory Turbelin 29/12/1998 *
*****************************************************
GRAPH='N' ;
COMPLET=FAUX ;
'OPTION' 'DIME' 2 'ELEM' qua8;
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' ;
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 ;
FINPROC ;
** PROCEDURES DE CALCUL
**--------------------
'DEBPROC' VMPM;
'ARGUMENT' rvx*'TABLE';
rv= rvx.'EQEX';
DD = rv.PASDETPS.'NUPASDT' ;
NN = DD/20;
L0 = (DD '-' (20*NN)) 'EGA' 0;
'SI' L0;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
** compteur pour connaitre le nombre de valeurs stockées
cont = 1.;
rv.INCO.'cont'= (rv.INCO.'cont' '+' cont);
** CALCUL DE LA VITESSE MOYENNE LE LONG DE L'AXE
eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
rv.INCO.'UXa_12'= (rv.INCO.'UXa_12' '+' UXa_12 );
rv.INCO.'UXa_34'= (rv.INCO.'UXa_34' '+' UXa_34 );
** CALCUL DE LA PRESSION MOYENNE SUR L'OBSTACLE
ePobst = 'EVOL' chpo PP scal ('INVERSE' obstacle);
Pobst = 'EXTRAIRE' ePobst ordo;
rv.INCO.'Pobst'= (rv.INCO.'Pobst' '+' Pobst);
** CALCUL DE L ENERGIE MOYENNE LE LONG DE L'AXE
eka_12 = EVOL CHPO (rv.inco.'KN') scal axe_12;
ka_12 = 'EXTRAIRE' eKa_12 ORDO;
eKa_34 = EVOL CHPO (rv.inco.'KN') scal axe_34;
ka_34 = 'EXTRAIRE' eka_34 ordo;
rv.INCO.'Ka_12'= (rv.INCO.'Ka_12' '+' Ka_12 );
rv.INCO.'Ka_34'= (rv.INCO.'Ka_34' '+' Ka_34 );
'FINSI' ;
'FINPROC' ;
**------------------------------------------------
'DEBPROC' CLCD;
'ARGUMENT' rvx*'TABLE';
rv = rvx.'EQEX';
DD = rv.PASDETPS.'NUPASDT' ;
NN = DD/20;
L0 = (DD '-' (20*NN)) 'EGA' 0;
'SI' L0;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
** CALCUL DES COEFF AERODYNAMIQUES
PPSF1 = 'REDU' PP face1;
PPCf1 = NOEL $face1 ('KCHT' $face1 SCAL SOMMET PPSF1);
PPSF2 = 'REDU' PP face2;
PPCf2 = NOEL $face2 ('KCHT' $face2 SCAL SOMMET PPSF2);
PPSF3 = 'REDU' PP face3;
PPCf3 = NOEL $face3 ('KCHT' $face3 SCAL SOMMET PPSF3);
PPSF4 = 'REDU' PP face4;
PPCf4 = NOEL $face4 ('KCHT' $face4 SCAL SOMMET PPSF4);
i=1;
CD1 = 0.;
CD3 = 0.;
'REPETER' bloc1 (NBEL $face1.centre) ;
PPf1i = 'EXTRAIRE' PPCf1 SCAL ($face1.centre 'POIN' i) ;
Lf1i = 'ABS' ( (coor 1 (face1 'POIN' (i+1))) '-'
(coor 1 (face1 'POIN' i) ) );
CD1i= PPf1i * Lf1i;
CD1 = CD1 '+'CD1i;
*
PPf3i = 'EXTRAIRE' PPCf3 SCAL ($face3.centre POIN i);
Lf3i = 'ABS' ( (coor 1 (face3 'POIN' i)) '-'
(coor 1 (face3 'POIN' (i+1) )) );
CD3i= PPf3i * Lf3i;
CD3 = CD3 '+' CD3i;
i = i+1;
'FIN' bloc1;
CL = CD1 '-' CD3;
*
j=1;
CD2 = 0.;
CD4 = 0.;
'REPETER' bloc2 (NBEL $face2.centre);
PPf2j = 'EXTRAIRE' PPCf2 SCAL ($face2.centre POIN j);
Lf2j = 'ABS' ( (coor 2 (face2 'POIN' j)) '-'
(coor 2 (face2 'POIN' (j+1) )) );
CD2j= PPf2j * Lf2j;
CD2 = CD2 '+' CD2j;
*
PPf4j = 'EXTRAIRE' PPCf4 SCAL ($face4.centre POIN j);
Lf4j = 'ABS' ( (coor 2 (face4 POIN j)) -
(coor 2 (face4 POIN (j+1) )) );
CD4j= PPf4j * Lf4j;
CD4 = CD4 '+' CD4j;
j = j+1;
'FIN' bloc2;
*
CD = CD4 '-' CD2 ;
rv.inco.'CL' = (rv.inco.'CL') et (Prog CL);
rv.inco.'CD' = (rv.inco.'CD') et (Prog CD);
PT = 'PROG' rv.'PASDETPS'.'TPS' ;
rv.INCO.'TP'= rv.INCO.'TP' ET PT;
'FINSI' ;
'FINPROC' ;
** Procedure pour determiner la valeur moyenne
** d'une liste de relles
'DEBPROC' MOYL LR*'LISTREEL';
LSi =0.;
i=1;
'REPETER' boucle ('DIME' LR);
ielem = 'EXTRAIRE' LR i;
LSi = LSi '+' ielem;
moyenne = LSi/i;
i =i+1;
'FIN' boucle;
'FINPROC' moyenne ;
** Ecart 'TYPE'
'DEBPROC' RMS LR*'LISTREEL';
LMi = 0.;
LEi = 0.;
i=1;
'REPETER' boucle ('DIME' LR);
ielem = 'EXTRAIRE' LR i;
ielem2 = (ielem**2);
LMi = LMi '+' ielem;
LEi = LEi '+' ielem2;
moyenne = LMi/i;
moyenne2 = LEi/i;
i =i+1;
'FIN' boucle;
ecart = (moyenne2 '-' ((moyenne)**2))**(0.5);
'FINPROC' ecart ;
** CONSTRUCTION DU MAILLAGE (ADIM)
**--------------------------------
*
**Longueurs sans dim
*
X1 = 4.5 ; X2 = 1. ; X3 = 14.5;
Y1 = 6.5 ; Y2 = 0.5 ; Y3 = 0.5 ; Y4=6.5 ;
*
**Points caracteristiques
*
P11=0. 0.;
P12=X1 0.;
P13=(X1+X2) 0.;
P14=(X1+X2+X3) 0.;
*
P21=0. Y1;
P22=X1 Y1;
P23=(X1+X2) Y1;
P24=(X1+X2+X3) Y1;
*
P31=0. (Y1+Y2);
P32=X1 (Y1+Y2);
P33=(X1+X2) (Y1+Y2);
P34=(X1+X2+X3) (Y1+Y2);
*
P41=0. (Y1+Y2+Y3);
P42=X1 (Y1+Y2+Y3);
P43=(X1+X2) (Y1+Y2+Y3);
P44=(X1+X2+X3) (Y1+Y2+Y3);
P51=0. (Y1+Y2+Y3+Y4);
P52=X1 (Y1+Y2+Y3+Y4);
P53=(X1+X2) (Y1+Y2+Y3+Y4);
P54=(X1+X2+X3) (Y1+Y2+Y3+Y4);
*
**Nb de mailles horizontales, sur obstacle
si complet ;
Nx = 6;
sinon ;
Nx = 2;
finsi ;
**Nb de mailles verticales 1/2 obstacle
si complet ;
Ny = 3;
sinon ;
Ny = 1;
finsi ;
**Nb de mailles horizontales, avant obstacle
Nx1= 3*Nx;
**Nb de mailles horizontales, sur obstacle
Nx2 = Nx ;
**Nb de mailles horizontales, après obstacle
Nx3 =7*Nx;
*
**Nb de mailles verticales sous obstacle
Ny1 = 8*Ny;
**Nb de mailles verticales 1/2 obstacle inf
Ny2 = Ny;
**Nb de mailles verticales 1/2 obstacle sup
Ny3 = Ny;
**Nb de mailles verticales au dessus obstacle
Ny4 = 8*Ny;
*
*
** Definition des densites
**-------------------------
DX1 = (X1/Nx1);
Dx1_1 = Dx1 ; Dx1_2 = 0.3*Dx1;
*
DX2 = (X2/Nx2);
Dx2_2 = Dx2 ; Dx2_3 = Dx2;
*
DX3 = (X3/Nx3);
Dx3_3 = 0.3*Dx3 ; Dx3_4 = 1.5*Dx3;
*
Dy1 = (Y1/NY1);
Dy1_1 =1.*DY1 ; Dy1_2 = 0.3*DY1;
*
DY2 = (Y2/NY2);
DY2_2 = DY2; DY2_3 = DY2;
*
DY3 = (Y3/NY3);
DY3_3 =0.5*DY3 ; DY3_4 = DY3;
*
DY4 = (Y4/NY4);
DY4_4 =0.3*DY4 ; DY4_5 = 1.*DY4;
*
**Definition des Lignes
**--------------------------------
L1x_12 = P11 droite (-1*Nx1) P12 'DINI' Dx1_1 'DFIN' Dx1_2;
L1x_23 = P12 droite (-1*Nx2) P13 'DINI' Dx2_2 'DFIN' Dx2_3;
L1x_34 = P13 droite (-1*Nx3) P14 'DINI' Dx3_3 'DFIN' Dx3_4;
L2x_12 = P21 droite (-1*Nx1) P22 'DINI' Dx1_1 'DFIN' Dx1_2;
L2x_23 = P22 droite (-1*Nx2) P23 'DINI' Dx2_2 'DFIN' Dx2_3;
L2x_34 = P23 droite (-1*Nx3) P24 'DINI' Dx3_3 'DFIN' Dx3_4;
L3x_12 = P31 droite (-1*Nx1) P32 'DINI' Dx1_1 'DFIN' Dx1_2;
L3x_23 = P32 droite (-1*Nx2) P33 'DINI' Dx2_2 'DFIN' Dx2_3;
L3x_34 = P33 droite (-1*Nx3) P34 'DINI' Dx3_3 'DFIN' Dx3_4;
L4x_12 = P41 droite (-1*Nx1) P42 'DINI' Dx1_1 'DFIN' Dx1_2;
L4x_23 = P42 droite (-1*Nx2) P43 'DINI' Dx2_2 'DFIN' Dx2_3;
L4x_34 = P43 droite (-1*Nx3) P44 'DINI' Dx3_3 'DFIN' Dx3_4;
L5x_12 = P51 droite (-1*Nx1) P52 'DINI' Dx1_1 'DFIN' Dx1_2;
L5x_23 = P52 droite (-1*Nx2) P53 'DINI' Dx2_2 'DFIN' Dx2_3;
L5x_34 = P53 droite (-1*Nx3) P54 'DINI' Dx3_3 'DFIN' Dx3_4;
L1x = L1x_12 'ET' L1x_23 'ET' L1x_34;
L2x = L2x_12 'ET' L2x_23 'ET' L2x_34;
L3x = L3x_12 'ET' L3x_23 'ET' L3x_34;
L4x = L4x_12 'ET' L4x_23 'ET' L4x_34;
L5x = L5x_12 'ET' L5x_23 'ET' L5x_34;
L1y_12 = P11 droite (-1*Ny1) P21 'DINI' Dy1_1 'DFIN' Dy1_2;
L1y_23 = P21 droite (-1*Ny2) P31 'DINI' Dy2_2 'DFIN' Dy2_3;
L1y_34 = P31 droite (-1*Ny3) P41 'DINI' Dy3_3 'DFIN' Dy3_4;
L1y_45 = P41 droite (-1*Ny4) P51 'DINI' Dy4_4 'DFIN' Dy4_5;
L2y_12 = P12 droite (-1*Ny1) P22 'DINI' Dy1_1 'DFIN' Dy1_2;
L2y_23 = P22 droite (-1*Ny2) P32 'DINI' Dy2_2 'DFIN' Dy2_3;
L2y_34 = P32 droite (-1*Ny3) P42 'DINI' Dy3_3 'DFIN' Dy3_4;
L2y_45 = P42 droite (-1*Ny4) P52 'DINI' Dy4_4 'DFIN' Dy4_5;
L3y_12 = P13 droite (-1*Ny1) P23 'DINI' Dy1_1 'DFIN' Dy1_2;
L3y_23 = P23 droite (-1*Ny2) P33 'DINI' Dy2_2 'DFIN' Dy2_3;
L3y_34 = P33 droite (-1*Ny3) P43 'DINI' Dy3_3 'DFIN' Dy3_4;
L3y_45 = P43 droite (-1*Ny4) P53 'DINI' Dy4_4 'DFIN' Dy4_5;
L4y_12 = P14 droite (-1*Ny1) P24 'DINI' Dy1_1 'DFIN' Dy1_2;
L4y_23 = P24 droite (-1*Ny2) P34 'DINI' Dy2_2 'DFIN' Dy2_3;
L4y_34 = P34 droite (-1*Ny3) P44 'DINI' Dy3_3 'DFIN' Dy3_4;
L4y_45 = P44 droite (-1*Ny4) P54 'DINI' Dy4_4 'DFIN' Dy4_5;
L1y = L1y_12 'ET' L1y_23 'ET' L1y_34 'ET' L1y_45;
L2y = L2y_12 'ET' L2y_23 'ET' L2y_34 'ET' L2y_45;
L3y = L3y_12 'ET' L3y_23 'ET' L3y_34 'ET' L3y_45;
L4y = L4y_12 'ET' L4y_23 'ET' L4y_34 'ET' L4y_45;
bas = (L1x) 'COULEUR' vert ; 'ELIMINATION' bas 1.e-5;
sortie = (L4y) 'COULEUR' jaune; 'ELIMINATION' sortie 1.e-5;
haut = (L5x) 'COULEUR' bleu ; 'ELIMINATION' haut 1.e-5;
entree = (L1y) 'COULEUR' rouge; 'ELIMINATION' entree 1.e-5;
axe_12 = L3x_12; axe_34=L3x_34;
axe = (axe_12 'ET' axe_34) 'COULEUR' turq; 'ELIMINATION' axe 1.e-5;
face1 = L2x_23 'COULEUR' vert ;
face2 = (L3y_23 'ET' L3y_34) 'COULEUR' jaune;
face3 = L4x_23 'COULEUR' bleu ;
face4 = (L2y_23 'ET' L2y_34) 'COULEUR' rouge ;
obstacle = FACE1 'ET' FACE2 'ET' FACE3 'ET' FACE4 ;
'ELIMINATION' obstacle 1.e-5; ;
**Creation de sous-zones
ZX1Y1 = 'DALLER' L1x_12 L2y_12 L2x_12 L1y_12 PLAN ;
ZX1Y2 = 'DALLER' L2x_12 L2y_23 L3x_12 L1y_23 PLAN ;
ZX1Y3 = 'DALLER' L3x_12 L2y_34 L4x_12 L1y_34 PLAN ;
ZX1Y4 = 'DALLER' L4x_12 L2y_45 L5x_12 L1y_45 PLAN ;
ZX1 = ZX1Y1 'ET' ZX1Y2 'ET' ZX1Y3 'ET' ZX1Y4;
'ELIMINATION' 1.e-5 ZX1 ;
ZX2Y1 = 'DALLER' L1x_23 L3y_12 L2x_23 L2y_12 PLAN ;
ZX2Y4 = 'DALLER' L4x_23 L3y_45 L5x_23 L2y_45 PLAN ;
ZX2 = ZX2Y1 'ET' ZX2Y4;
'ELIMINATION' ZX2 1.e-5;
orienter Zx2;
ZX3Y1 = 'DALLER' L1x_34 L4y_12 L2x_34 L3y_12 PLAN ;
ZX3Y2 = 'DALLER' L2x_34 L4y_23 L3x_34 L3y_23 PLAN ;
ZX3Y3 = 'DALLER' L3x_34 L4y_34 L4x_34 L3y_34 PLAN ;
ZX3Y4 = 'DALLER' L4x_34 L4y_45 L5x_34 L3y_45 PLAN ;
ZX3 = ZX3Y1 'ET' ZX3Y2 'ET' ZX3Y3 'ET' ZX3Y4;
'ELIMINATION' ZX3 1.e-5;
orienter Zx3;
SURFA = Zx1 'ET' ZX2 'ET' Zx3 ;
'ELIMINATION' surfa 1.e-5;
orienter surfa;
*'TRACER' surfa;
*'OPTION' donn 5;
** Définition des domaines
**------------------------
$surfa = 'DOMA' surfa 1.e-5 MACRO;
$bas = 'DOMA' bas INCL $surfa 1.e-5 MACRO;
$entree = 'DOMA' entree INCL $surfa 1.e-5 MACRO;
$haut = 'DOMA' haut INCL $surfa 1.e-5 MACRO;
$sortie = 'DOMA' sortie INCL $surfa 1.e-5 MACRO;
$face1 = 'DOMA' face1 INCL $surfa 1.e-5 MACRO;
$face2 = 'DOMA' face2 INCL $surfa 1.e-5 MACRO;
$face3 = 'DOMA' face3 INCL $surfa 1.e-5 MACRO;
$face4 = 'DOMA' face4 INCL $surfa 1.e-5 MACRO;
$axe_12 = 'DOMA' axe_12 INCL $surfa 1.e-5 MACRO;
$axe_34 = 'DOMA' axe_34 INCL $surfa 1.e-5 MACRO;
surfa=$surfa.maillage;
'ORIENTER' surfa;
entree = $entree.maillage;
haut = $haut.maillage ;
bas = $bas.maillage ;
sortie = $sortie.maillage ;
face1 = $face1.maillage ;
face2 = $face2.maillage ;
face3 = $face3.maillage ;
face4 = $face4.maillage ;
axe_34 = $axe_34.maillage ;
axe_12 = $axe_12.maillage ;
** Proprietes physiques
**-------------------
nu_eau = 1.e-6;
nu = nu_eau;
** Longeur de reference = Diametre du cylindre en m
D_ref = 0.04;
** Hauteur du domaine (sans dim)
Hda = 14.;
** Longueur du domaine (sans dim)
Lda = 20.;
** Reynolds de l'écoulement
Re = 2.2e4;
iRe = 1/Re ;
** Vitesse de reference = vitesse en entree (m/s)
U_ref = (Re*nu)/D_ref;
** Distances Normales aux parois de l'obstacle
**---------------------------------------------
** Coordonnees des points des faces du cylindre
Xface1 = coor 1 face1;
Yface2 = coor 2 face2;
Xface3 = coor 1 face3;
Yface4 = coor 2 face4;
** abscisse curviligne le long des faces
eXface1 = evol chpo Xface1 face1;
eXface1 = 'EXTRAIRE' eXface1 ABSC;
eYface2 = evol chpo Yface2 face2;
eYface2 = 'EXTRAIRE' eYface2 ABSC;
eXface3 = evol chpo Xface3 face3;
eXface3 = 'EXTRAIRE' eXface3 ABSC;
eYface4 = evol chpo Yface4 face4;
eYface4 = 'EXTRAIRE' eYface4 ABSC;
** Yp = 1/5 hauteur de la premiere maille (sans dim)
YP1 = ('EXTR' eXface1 2)/5;
YP2 = ('EXTR' eYface2 2)/5;
YP3 = ('EXTR' eXface3 2)/5;
YP4 = ('EXTR' eYface4 2)/5;
** vitesse de frottement sans dim
uf1 = (150*nu)/yp1;
uf2 = (150*nu)/yp2;
uf3 = (150*nu)/yp3;
uf4 = (150*nu)/yp4;
**constante de Von Karman
C1 = 0.4;
** constante du modele RNG k-e
CNU = 0.0845;
** taux de turbulence
INT = 0.02 ;
** Definition des profils d'entrees et initiaux
**----------------------------------------------
** Definition du domaine d'imposition
M_imp =entree;
** VITESSES
**-----------
U1 = 1. 0.;
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 0.;
U_ini = U_ini 'NOMC' 'UX' ;
U_ini = 'KCHT' $surfa vect sommet U_ini;
**ENERGIE CINETIQUE
**-----------------
** Kadim = K/(U_ref**2)
** En fonction de l'intensite de la turbulence
K1 = 1.5*(INT)**2 ;
K1 = 'MANU' 'CHPO' surfa 1 scal K1;
***INITIALE
K_ini = 'KCHT' $surfa scal sommet K1;
***IMPOSEE
K_IMP = 'KCHT' $surfa scal sommet K1;
K_IMP= 'EXCO' ('REDU' K_IMP M_imp) SCAL;
**DISSIPATION Eadim = E*Href/Uref**3
**--------------
** Fonction de K et de D
** ! K est sans dim et D = Href
E1 = (K1**1.5);
** INITIALE
E_ini = 'KCHT' $surfa scal sommet E1;
**IMPOSEE
E_IMP = 'KCHT' $surfa scal sommet E1;
E_IMP='EXCO' ('REDU' E_IMP M_imp) SCAL;
**NUT INITIALE
**------------------
*fonction de K et E
N1=CNU*(K_ini**2)*(E_ini**-1);
N1 = 'KCHT' $surfa scal sommet N1;
N_ini = 'NOEL' $surfa N1 ;
** U* INITIALES
**-------------
UET_ini1 = 'KCHT' $face1 scal centre Uf1;
UET_ini2 = 'KCHT' $face2 scal centre Uf2;
UET_ini3 = 'KCHT' $face3 scal centre Uf3;
UET_ini4 = 'KCHT' $face4 scal centre Uf4;
** Module de resolution
**-------------------
RV = 'EQEX' $surfa ITMA 42 ALFA 0.8
OPTI 'RNG'
ZONE $surfa OPER FILTREKE 1. Hda iRe 'UN' INCO 'KN' 'EN'
ZONE $surfa OPER NSKE iRe 'NUT' INCO 'UN' 'KN' 'EN'
ZONE $face4 OPER FPU iRe 'UET4' YP4 INCO 'UN' 'KN' 'EN'
ZONE $face2 OPER FPU iRe 'UET2' YP2 INCO 'UN' 'KN' 'EN'
ZONE $face3 OPER FPU iRe 'UET3' YP3 INCO 'UN' 'KN' 'EN'
ZONE $face1 OPER FPU iRe 'UET1' YP1 INCO 'UN' 'KN' 'EN' ;
RV = 'EQEX' RV
ZONE $surfa OPER CLCD
ZONE $surfa OPER VMPM;
RV = 'EQEX' RV
CLIM 'UN' UIMP M_imp UX_imp
CLIM 'UN' VIMP M_imp UY_imp
CLIM 'UN' VIMP (haut 'ET' bas) 0.
CLIM 'KN' TIMP M_imp K_imp
CLIM 'EN' TIMP M_imp E_imp ;
RVP = 'EQPR' $surfa
ZONE $face4 OPER VNIMP 0.
ZONE $face2 OPER VNIMP 0.
ZONE $face3 OPER VNIMP 0.
ZONE $face1 OPER VNIMP 0.
ZONE $surfa OPER PRESSION ;
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.'UET4' = UET_ini4;
RV.INCO.'UET3' = UET_ini3;
RV.INCO.'UET2' = UET_ini2;
RV.INCO.'UET1' = UET_ini1;
RV.INCO.'CL' = 'PROG' 0.;
RV.INCO.'CD' = 'PROG' 0.;
RV.INCO.'TP' = 'PROG' 0.;
RV.INCO.'cont' = 0.;
RV.INCO.'UXa_12'= 'PROG' ('NBNO' axe_12)*0. ;
RV.INCO.'UXa_34'= 'PROG' ('NBNO' axe_34)*0. ;
RV.INCO.'Pobst' = 'PROG' (('NBNO' obstacle) '+' 1 )*0. ;
rv.INCO.'Ka_12'= 'PROG' ('NBNO' axe_12)*0. ;
rv.INCO.'Ka_34'= 'PROG' ('NBNO' axe_34)*0. ;
lh= (POIN SURFA PROC (2. 7.)) 'ET' (POIN SURFA PROC (6. 7.))
'ET' (POIN SURFA PROC (8. 7.)) 'ET' (POIN SURFA PROC (10. 7.))
'ET' (POIN SURFA PROC (13. 7.)) 'ET' (POIN SURFA PROC (18. 7.));
his= 'KHIS' 'UN' 1 lh
'UN' 2 lh
'KN' lh
'EN' lh;
rv.hist=his;
EXEC RV;
* 'OPTION' donn 5;
* 'FIN' ;
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;
* 'OPTION' donn 5;
*****HISTORIQUES*****
*
TAB1=his.TABD;
TAB1.'TITRE'= TABLE ;
TAB1.'TITRE' . 1 = MOT 'COURBE 1';
TAB1.'TITRE' . 2 = MOT 'COURBE 2';
TAB1.'TITRE' . 3 = MOT 'COURBE 3';
TAB1.'TITRE' . 4 = MOT 'COURBE 4';
TAB1.'TITRE' . 5 = MOT 'COURBE 5';
TAB1.'TITRE' . 6 = MOT 'COURBE 6';
si ('EGA' graph 'O' );
'DESSIN' (his.'1UN') 'TITRE' 'HISTORIQUES POUR UX ' LEGE TAB1;
'DESSIN' (his.'2UN') 'TITRE' 'HISTORIQUES POUR UY ' LEGE TAB1;
'DESSIN' (his.'KN') 'TITRE' 'HISTORIQUES POUR K' LEGE TAB1;
'DESSIN' (his.'EN') 'TITRE' 'HISTORIQUES POUR E' LEGE TAB1;
finsi ;
** 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;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
si ('EGA' graph 'O' );
*'OPTION' isov ligne;
titre 'Vecteur Vitesse';
'TRACER' vv surfa ('CONTOUR' surfa);
titre 'vecteur vitesse composante X';
'TRACER' Vx surfa ('CONTOUR' surfa) ;
titre 'vecteur vitesse composante Y';
'TRACER' VY surfa ('CONTOUR' surfa) ;
titre 'CHAMP DE PRESSION';
'TRACER' PP surfa ('CONTOUR' surfa);
titre 'Energie Cinetique Turbulente';
'TRACER' kk surfa ('CONTOUR' surfa);
titre 'dissipation';
'TRACER' ee surfa ('CONTOUR' surfa);
titre 'NUT';
'TRACER' nuts surfa ('CONTOUR' surfa);
finsi ;
* 'OPTION' donn 5;
** 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
unnh=KCHT $HAUT SCAL SOMMET (un1);
unnh=NOEL $HAUT unnh;
unnh=KCHT $HAUT SCAL CENTRE unnh;
unnb=KCHT $bas SCAL SOMMET (-1*un1);
unnb=NOEL $bas unnb;
unnb=KCHT $bas SCAL CENTRE unnb;
unne=KCHT $entree SCAL SOMMET (un2);
unne=NOEL $entree unne;
unne=KCHT $entree SCAL CENTRE unne;
unns=KCHT $sortie SCAL SOMMET (-1*un2);
unns=NOEL $sortie unns;
unns=KCHT $sortie SCAL CENTRE unns;
rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
ZONE $SURFA OPER LAPN 1. inco 'psi'
ZONE $SURFA OPER FIMP sw inco 'psi'
ZONE $HAUT OPER FIMP unnh inco 'psi'
ZONE $bas OPER FIMP unnb inco 'psi'
ZONE $sortie OPER FIMP unns inco 'psi'
ZONE $entree OPER FIMP unne inco 'psi'
CLIM 'psi' TIMP obstacle 0. ;
rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
exec rk;
si ('EGA' graph 'O' );
psi=rk.inco.'psi';
'OPTION' ISOV LIGNE;
toto = 'PROG' -0.51 PAs 0.02 Npas 50;
'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
finsi ;
** CALCUL DES COEFFS AERODYNAMIQUES
**------------------------------
** Trainee 'ET' portance
cd = 2*(rv.inco.'CD');
cl = 2*(rv.inco.'CL');
tp = 'ENLEVER' (rv.inco.'TP') 1;
cd = 'ENLEVER' cd 1;
cl = 'ENLEVER' cl 1;
Ntp = ('DIME' tp);
PDT1 = 'EXTRAIRE' tp 1;
PDTN = 'EXTRAIRE' tp Ntp;
** duree moyenne de 20 pas de tps
PDTx20 = (PDTN '-'PDT1 )/( Ntp '-' 1);
** Pas de tps moyen
PDTmoy = (PDTx20 '/' 20);
*'LISTE' PDTmoy;
*'OPTION' donn 5;
ecd = 'EVOL' 'MANUEL' tp cd;
ecl = 'EVOL' 'MANUEL' tp cl;
moycd = MOYL cd;
moycl = MOYL cl;
rmscd = RMS cd;
rmscl = RMS cl;
'LISTE' moycd; 'LISTE' rmscd; 'LISTE' moycl; 'LISTE' rmscl;
si ('EGA' graph 'O' );
'DESSIN' ecl 'TITRE' 'cl';
'DESSIN' ecd 'TITRE' 'cd';
DSPCL = DSPR 8 ecl Fmax 0.8;
'DESSIN' DSPCL;
finsi ;
** Evolution du champ de pression autour de l'obstacle
** Presion instantanee
Pref = 'EXTRAIRE' PP Scal P32;
Pref = Pref '-' 0.5;
Pobst = 'REDU' PP obstacle;
CP = 2*(Pobst '-' Pref);
eCP = 'EVOL' chpo CP scal ('INVERSE' obstacle) ;
** On enleve les valeurs aux extremites
Oobst = 'EXTRAIRE' eCP ordo;
Oobst = 'ENLEVER' Oobst 1 ;
Aobst = 'EXTRAIRE' eCP absc;
Aobst = 'ENLEVER' Aobst 1 ;
NbAobst = 'DIME' Aobst ;
Aobst = 'ENLEVER' Aobst NbAobst;
Oobst = 'ENLEVER' Oobst NbAobst;
eCP = 'EVOL' 'MANUEL' Aobst Oobst;
si ('EGA' graph 'O' );
'DESSIN' (eCP) 'TITRE' 'pression instantanee sur obstacle' ;
finsi ;
** Pression moyenne le long de l'obstacle
PM =(rv.INCO.'Pobst') '/' (rv.INCO.'cont');
PM = 2*(PM '-'('PROG'('DIME' PM)*Pref ));
PM = 'ENLEVER' PM 1;
PM = 'ENLEVER' PM ('DIME' Pm);
ePM = 'EVOL' 'MANUEL' Aobst PM;
si ('EGA' graph 'O' );
'DESSIN' (epm) 'TITRE' 'Pression moyenne sur obstacle';
finsi ;
** EVOLUTIONS LE LONG DE L'AXE
**-------------------------------
** ABSC reelles des points de l'axe
Xaxe_12 = 'EVOL' chpo (COOR 1 axe_12) scal axe_12;
Xaxe_34 = 'EVOL' chpo (COOR 1 axe_34) scal axe_34;
Xaxe_12 = 'EXTRAIRE' Xaxe_12 ordo;
Xaxe_34 = 'EXTRAIRE' Xaxe_34 ordo;
** Pour tracer Uy
XUY_34 = ( 'PROG' 5.49) 'ET' Xaxe_34;
** Evolution instantanne le long de l'axe (fct de l'abscisse curviligne)
eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
eUYa_12 = EVOL CHPO (rv.inco.'UN') UY axe_12;
eKa_12 = EVOL CHPO (rv.inco.'KN') scal axe_12;
eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
eUYa_34 = EVOL CHPO (rv.inco.'UN') uy axe_34;
eKa_34 = EVOL CHPO (rv.inco.'KN') scal axe_34;
** valeurs le long de l'axe
UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
UYa_12 = 'EXTRAIRE' eUYa_12 ORDO;
Ka_12 = 'EXTRAIRE' eKa_12 ordo;
UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
UYa_34 = 'EXTRAIRE' eUYa_34 ordo;
UYa_34 = ('PROG' 0.) 'ET' UYa_34;
Ka_34 = 'EXTRAIRE' eKa_34 ordo;
** Evolution le long de l'axe (fct de l'abscisse relle)
eUXa_12 = EVOL VERT 'MANUEL' Xaxe_12 UXa_12;
eUXa_34 = EVOL VERT 'MANUEL' Xaxe_34 UXa_34;
eUXaxe = eUXa_12 'ET' eUXa_34;
eUya_12 = EVOL VERT 'MANUEL' Xaxe_12 Uya_12;
eUya_34 = EVOL VERT 'MANUEL' XUY_34 Uya_34;
eUyaxe = eUya_12 'ET' eUya_34;
eKa_12 = EVOL VERT 'MANUEL' Xaxe_12 Ka_12;
eKa_34 = EVOL VERT 'MANUEL' Xaxe_34 Ka_34;
eKaxe = eKa_12 'ET' eKa_34;
** Ux moyen le long de l'axe
UM_12 = (rv.INCO.'UXa_12') '/' (rv.INCO.'cont');
UM_34 = (rv.INCO.'UXa_34') '/' (rv.INCO.'cont');
eUM_12 = EVOL VERT 'MANUEL' Xaxe_12 UM_12;
eUM_34 = EVOL VERT 'MANUEL' Xaxe_34 UM_34;
eUMaxe = eUM_12 'ET' eUM_34;
** K moyen le long de l'axe
KM_12 = (rv.INCO.'Ka_12') '/' (rv.INCO.'cont');
KM_34 = (rv.INCO.'Ka_34') '/' (rv.INCO.'cont');
eKM_12 = EVOL VERT 'MANUEL' Xaxe_12 KM_12;
eKM_34 = EVOL VERT 'MANUEL' Xaxe_34 KM_34;
eKMaxe = eKM_12 'ET' eKM_34;
** Axes
LX0_12 = 'PROG' 0. 4.5; LX0_34 = 'PROG' 5.5 20.;
LY0 = 'PROG' 0. 0.;
L0_12 = 'EVOL' 'MANUEL' LX0_12 Ly0;
L0_34 = 'EVOL' 'MANUEL' LX0_34 Ly0;
L0 = L0_12 'ET' L0_34;
si ('EGA' graph 'O' );
'DESSIN' (eUMaxe 'ET' L0 ) 'TITRE' 'Ux moyen sur axe';
'DESSIN' (eUXaxe 'ET' L0) titr 'PROFIL Ux sur axe'
tity 'Ux/U0' titx 'X/D';
'DESSIN' (eUya_34 'ET' L0_34) titr 'PROFIL Uy sur axe'
tity 'Uy/U0' titx 'X/D';
'DESSIN' (eKMaxe 'ET' L0 ) 'TITRE' 'K moyen sur axe';
'DESSIN' (eKaxe 'ET' L0) titr 'PROFIL K sur axe'
tity 'K/U0**2' titx 'X/D';
finsi ;
Fin;
** Pas de tps moyen
*PDTmoy = (PDTx20 '/' 20);
*'LISTE' PDTmoy;
*'OPTION' donn 5;
** Evolution reguliere dans le tps
* RegTPS = 'PROG' PDT1 PAS PDTx20 Npas (Nval '-' 1) ;
*'OPTION' donn 5;
*
UYFDT_1= 'EXTRAIRE' his.'2UN' ORDO 1;
UYFDT_2= 'EXTRAIRE' his.'2UN' ORDO 2;
UYFDT_3= 'EXTRAIRE' his.'2UN' ORDO 3;
UYFDT_4= 'EXTRAIRE' his.'2UN' ORDO 4;
UYFDT_5= 'EXTRAIRE' his.'2UN' ORDO 5;
UYFDT_6= 'EXTRAIRE' his.'2UN' ORDO 6;
*
eUYFDT_1 = 'EVOL' 'MANUEL' hisTPS UYFDT_1;
eUYFDT_2 = 'EVOL' 'MANUEL' hisTPS UYFDT_2;
eUYFDT_3 = 'EVOL' 'MANUEL' hisTPS UYFDT_3;
eUYFDT_4 = 'EVOL' 'MANUEL' hisTPS UYFDT_4;
eUYFDT_5 = 'EVOL' 'MANUEL' hisTPS UYFDT_5;
eUYFDT_6 = 'EVOL' 'MANUEL' hisTPS UYFDT_6;
UxFDT_1= 'EXTRAIRE' his.'1UN' ORDO 1;
UxFDT_2= 'EXTRAIRE' his.'1UN' ORDO 2;
UXFDT_3= 'EXTRAIRE' his.'1UN' ORDO 3;
UxFDT_4= 'EXTRAIRE' his.'1UN' ORDO 4;
UxFDT_5= 'EXTRAIRE' his.'1UN' ORDO 5;
UxFDT_6= 'EXTRAIRE' his.'1UN' ORDO 6;
*
eUxFDT_1 = 'EVOL' 'MANUEL' hisTPS UxFDT_1;
eUxFDT_2 = 'EVOL' 'MANUEL' hisTPS UxFDT_2;
eUxFDT_3 = 'EVOL' 'MANUEL' hisTPS UxFDT_3;
eUxFDT_4 = 'EVOL' 'MANUEL' hisTPS UxFDT_4;
eUxFDT_5 = 'EVOL' 'MANUEL' hisTPS UxFDT_5;
eUxFDT_6 = 'EVOL' 'MANUEL' hisTPS UxFDT_6;
*
*TABUY=Table;
* TABUY.1 ='MARQ ETOI';
'DESSIN' eUYFDT_1 'TITRE''EVOLUTION UY_axe en X=1.94'
tity 'UY' titx 'tps';
'DESSIN' eUYFDT_2 'TITRE''EVOLUTION UY_axe en X=5.9'
tity 'UY' titx 'tps' ;
'DESSIN' eUYFDT_3 'TITRE''EVOLUTION UY_axe en X=8'
tity 'UY' titx 'tps'
'DESSIN' eUYFDT_4 'TITRE''EVOLUTION UY_axeen X=10'
tity 'UY' titx 'tps' ;
'DESSIN' eUYFDT_5 'TITRE''EVOLUTION UY_axe en fonction du temps
en X=13' tity 'UY' titx 'tps' ;
'DESSIN' eUYFDT_6 'TITRE''EVOLUTION UY_axe en fonction du temps
en X=18' tity 'UY' titx 'tps' ;
'DESSIN' eUxFDT_1 'TITRE''EVOLUTION Ux_axe en fonction du temps
en X=1.94' tity'Ux' titx 'tps' ;
'DESSIN' eUxFDT_2 'TITRE''EVOLUTION Ux_axe en fonction dutemps
en X=5.9' tity 'Ux' titx 'tps' ;;
'DESSIN' eUxFDT_3 'TITRE''EVOLUTION Ux_axe en fonction du temps
en X=8' tity 'Ux' titx 'tps' ;;
'DESSIN' eUxFDT_4 'TITRE''EVOLUTION Ux_axe en fonction du temps
en X=10' tity 'TABUY Ux' titx 'tps' ;;
'DESSIN' eUxFDT_5 'TITRE''EVOLUTION Ux_axe en fonction du temps
en X=13' tity 'Ux' titx 'tps' ;;
'DESSIN' eUxFDT_6 'TITRE''EVOLUTION Ux_axe en fonction du temps
en X=18' tity 'Ux' titx 'tps' ;;
'OPTION' donn 5;
*********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' ;
kk_ent = EXCO(REDU kk entree) SCAL ;
EE = rv.inco.'EN';
ee_ent = 'EXCO' ('REDU' ee entree) SCAL ;
nuts = 'ELNO' ('KCHT' $surfa scal centre (RV.INCO.'NUT')) $surfa;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
*'OPTION' isov ligne;
* titre 'Vecteur Vitesse';
* 'TRACER' vv surfa ('CONTOUR' surfa);
* titre 'vecteur vitesse composante X';
* 'TRACER' Vx surfa ('CONTOUR' surfa) ;
* titre 'vecteur vitesse composante Y';
* 'TRACER' VY surfa ('CONTOUR' surfa) ;
* titre 'CHAMP DE PRESSION';
* 'TRACER' PP surfa ('CONTOUR' surfa);
* titre 'Energie Cinetique Turbulente';
* 'TRACER' kk surfa ('CONTOUR' surfa);
* titre 'dissipation';
* 'TRACER' ee surfa ('CONTOUR' surfa);
* titre 'NUT';
* 'TRACER' nuts surfa ('CONTOUR' surfa);
* 'OPTION' donn 5;
** FONCTION DE COURANT
** PREMIERE METHODE
**-------------------
** CALCUL DU ROTATIONNEL DU CHAMP DE VITESSE
* sw ='KOPS' (rv.inco.'UN') 'ROT' $surfa;
* DOMLIM = Haut 'ET' BAS 'ET' entree ;
* 'ELIMINATION' domliM 1.e-5 ;
* CHPY = 'COORDONNEE' 2 DOMLIM ;
* PSin = CHPY ;
* CHPTOT = 'KCHT' $surfa scal sommet psin ;
* rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
* ZONE $SURFA OPER LAPN 1. inco 'psi'
* ZONE $SURFA OPER FIMP sw inco 'psi'
* CLIM 'psi' TIMP DOMLIM CHPTOT ;
* rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
* exec rk;
* psi=rk.inco.'psi';
*
* 'OPTION' ISOV LIGNE;
* 'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
* 'OPTION' donn 5 ;
** FONCTION DE COURANT
** SECONDE METHODE
**---------------------
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 !!
unnh=KCHT $HAUT SCAL SOMMET (un1);
unnh=NOEL $HAUT unnh;
unnh=KCHT $HAUT SCAL CENTRE unnh;
unnb=KCHT $bas SCAL SOMMET (-1*un1);
unnb=NOEL $bas unnb;
unnb=KCHT $bas SCAL CENTRE unnb;
unne=KCHT $entree SCAL SOMMET (un2);
unne=NOEL $entree unne;
unne=KCHT $entree SCAL CENTRE unne;
unns=KCHT $sortie SCAL SOMMET (-1*un2);
unns=NOEL $sortie unns;
unns=KCHT $sortie SCAL CENTRE unns;
rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
ZONE $SURFA OPER LAPN 1. inco 'psi'
ZONE $SURFA OPER FIMP sw inco 'psi'
ZONE $HAUT OPER FIMP unnh inco 'psi'
ZONE $bas OPER FIMP unnb inco 'psi'
ZONE $sortie OPER FIMP unns inco 'psi'
ZONE $entree OPER FIMP unne inco 'psi'
CLIM 'psi' TIMP obstacle 0. ;
rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
exec rk;
psi=rk.inco.'psi';
*
'OPTION' ISOV LIGNE;
toto = 'PROG' -0.51 PAs 0.02 Npas 50;
'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
'OPTION' donn 5 ;
* 'FIN' ;
*
** CALCUL DES COEFFS AERODYNAMIQUES
**------------------------------
** Traine 'ET' portance
cd = 2*(rv.inco.'CD');
cl = 2*(rv.inco.'CL');
tp = 'ENLEVER' (rv.inco.'TP') 1;
cd = 'ENLEVER' cd 1;
cl = 'ENLEVER' cl 1;
cd = tronc cd 400 ;
cl = tronc cl 400 ;
tp = tronc tp 400 ;
Ntp = ('DIME' tp);
PDT1 = 'EXTRAIRE' tp 1;
PDTN = 'EXTRAIRE' tp Ntp;
** duree moyenne de 20 pas de tps
PDTx20 = (PDTN '-'PDT1 )/( Ntp '-' 1);
** Pas de tps moyen
PDTmoy = (PDTx20 '/' 20);
'LISTE' PDTmoy;
'OPTION' donn 5;
** Evolution reguliere dans le tps
RegTP = 'PROG' PDT1 PAS PDTx20 Npas (Ntp '-' 1) ;
*'OPTION' donn 5;
ecd = 'EVOL' 'MANUEL' tp cd;
ecl = 'EVOL' 'MANUEL' tp cl;
moycd = MOYL cd;
moycl = MOYL cl;
rmscd = RMS cd;
rmscl = RMS cl;
'LISTE' moycd;
'LISTE' rmscd;
'LISTE' moycl;
'LISTE' rmscl;
'DESSIN' ecl 'TITRE' 'cl';
'DESSIN' ecd 'TITRE' 'cd';
DSPCL = DSPR 8 ecl Fmax 0.8;
'DESSIN' DSPCL;
'OPTION' donn 5;
* Evolution du champ de pression autour de l'obstacle
** Presion instantanee
Pref = 'EXTRAIRE' PP Scal P32;
Pref = Pref '-' 0.5;
Pobst = 'REDU' PP obstacle;
CP = 2*(Pobst '-' Pref);
eCP = 'EVOL' chpo CP scal ('INVERSE' obstacle) ;
* titi = 'PROG' Pref;
'OPTION' isov ligne;
titre 'CHAMP DE PRESSION';
* 'TRACER' PP surfa toto ('CONTOUR' surfa);
* 'TRACER' PP surfa ('CONTOUR' surfa);
* 'OPTION' donn 5;
** On enleve les valeurs aux extremites
Oobst = 'EXTRAIRE' eCP ordo;
Oobst = 'ENLEVER' Oobst 1 ;
Aobst = 'EXTRAIRE' eCP absc;
Aobst = 'ENLEVER' Aobst 1 ;
NbAobst = 'DIME' Aobst ;
Aobst = 'ENLEVER' Aobst NbAobst;
Oobst = 'ENLEVER' Oobst NbAobst;
eCP = 'EVOL' 'MANUEL' Aobst Oobst;
'DESSIN' (eCP) 'TITRE' 'pression instantanee sur obstacle' ;
** Pression moyenne le long de l'obstacle
PM =(rv.INCO.'Pobst') '/' (rv.INCO.'cont');
PM = 2*(PM '-'('PROG'('DIME' PM)*Pref ));
PM = 'ENLEVER' PM 1;
PM = 'ENLEVER' PM ('DIME' Pm);
ePM = 'EVOL' 'MANUEL' Aobst PM;
'DESSIN' (epm) 'TITRE' 'Pression moyenne sur obstacle';
* 'OPTION' donn 5;
**
** EVOLUTIONS LE LONG DE L'AXE
**-------------------------------
** ABSC reelles des points de l'axe
Xaxe_12 = 'EVOL' chpo (COOR 1 axe_12) scal axe_12;
Xaxe_34 = 'EVOL' chpo (COOR 1 axe_34) scal axe_34;
Xaxe_12 = 'EXTRAIRE' Xaxe_12 ordo;
* 'LISTE' Xaxe_12;
Xaxe_34 = 'EXTRAIRE' Xaxe_34 ordo;
** Pour tracer Uy
XUY_34 = ( 'PROG' 5.49) 'ET' Xaxe_34;
* 'LISTE' Xaxe_34;
* 'OPTION' donn 5;
** Evolution instantanne
** le long de l'axe (fct de l'abscisse curviligne)
eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
eUYa_12 = EVOL CHPO (rv.inco.'UN') UY axe_12;
eKa_12 = EVOL CHPO (rv.inco.'KN') scal axe_12;
eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
eUYa_34 = EVOL CHPO (rv.inco.'UN') uy axe_34;
eKa_34 = EVOL CHPO (rv.inco.'KN') scal axe_34;
**valeurs le long de l'axe
UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
UYa_12 = 'EXTRAIRE' eUYa_12 ORDO;
Ka_12 = 'EXTRAIRE' eKa_12 ordo;
UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
UYa_34 = 'EXTRAIRE' eUYa_34 ordo;
UYa_34 = ('PROG' 0.) 'ET' UYa_34;
Ka_34 = 'EXTRAIRE' eKa_34 ordo;
** Evolution le long de l'axe (fct de l'abscisse relle)
eUXa_12 = EVOL VERT 'MANUEL' Xaxe_12 UXa_12;
eUXa_34 = EVOL VERT 'MANUEL' Xaxe_34 UXa_34;
eUXaxe = eUXa_12 'ET' eUXa_34;
eUya_12 = EVOL VERT 'MANUEL' Xaxe_12 Uya_12;
eUya_34 = EVOL VERT 'MANUEL' XUY_34 Uya_34;
eUyaxe = eUya_12 'ET' eUya_34;
eKa_12 = EVOL VERT 'MANUEL' Xaxe_12 Ka_12;
eKa_34 = EVOL VERT 'MANUEL' Xaxe_34 Ka_34;
eKaxe = eKa_12 'ET' eKa_34;
** Ux moyen le long de l'axe
UM_12 = (rv.INCO.'UXa_12') '/' (rv.INCO.'cont');
UM_34 = (rv.INCO.'UXa_34') '/' (rv.INCO.'cont');
eUM_12 = EVOL VERT 'MANUEL' Xaxe_12 UM_12;
eUM_34 = EVOL VERT 'MANUEL' Xaxe_34 UM_34;
eUMaxe = eUM_12 'ET' eUM_34;
** K moyen le long de l'axe
KM_12 = (rv.INCO.'Ka_12') '/' (rv.INCO.'cont');
KM_34 = (rv.INCO.'Ka_34') '/' (rv.INCO.'cont');
eKM_12 = EVOL VERT 'MANUEL' Xaxe_12 KM_12;
eKM_34 = EVOL VERT 'MANUEL' Xaxe_34 KM_34;
eKMaxe = eKM_12 'ET' eKM_34;
* opti donn 5;
TAB2=TABLE;
* TAB2.'TITRE'= TABLE ;
* TAB2.1='MARQ ETOI NOLI';
* TAB2.'TITRE' . 1 = MOT 'Calcul';
* TAB2.2='MARQ ETOI NOLI';
* TAB2.'TITRE' . 2 = MOT 'en entree';
TAB2.3='TIRL';
* TAB2.'TITRE' . 3 = MOT 'en entree';
TAB2.4='TIRL';
** affichage du cube
LX0_12 = 'PROG' 0. 4.5;
LX0_34 = 'PROG' 5.5 20.;
LY0 = 'PROG' 0. 0.;
L0_12 = 'EVOL' 'MANUEL' LX0_12 Ly0;
L0_34 = 'EVOL' 'MANUEL' LX0_34 Ly0;
L0 = L0_12 'ET' L0_34;
CUBEX = 'PROG' 4.5 4.5 5.5 5.5 5.5 4.5 4.5;
CUBEX2 = 'PROG' 4.49 4.49 5.5 5.5 5.5 4.49 4.49;
CUBE1 = 'PROG' 0. 0.05 0.05 0. -0.05 -0.05 0. ;
CUBE2 = 'PROG' 0. 0.02 0.02 0. -0.02 -0.02 0. ;
CUBE3 = 'PROG' 0. 1e-3 1e-3 0. -1e-3 -1e-3 0. ;
CUBEux = 'EVOL' ROUGE 'MANUEL' CUBEx CUBE1;
CUBEuy = 'EVOL' ROUGE 'MANUEL' CUBEx2 CUBE2;
CUBEK = 'EVOL' ROUGE 'MANUEL' CUBEx CUBE3;
** SORTIES GRAPHIQUES
**-------------------
*** VITESSE
**-------------------------------
'DESSIN' (eUMaxe 'ET' L0 'ET' cubeux ) 'TITRE' 'Ux moyen sur axe';
'DESSIN' (eUXaxe 'ET' L0 'ET' cubeux) titr 'PROFIL Ux sur axe'
TAB2 tity 'Ux/U0' titx 'X/D';
* dess (eUyaxe 'ET' L0 'ET' cubeuy) titr 'PROFIL Uy sur axe'
* TAB2 tity 'Uy/U0' titx 'X/D';
dess (eUya_34 'ET' L0_34 'ET' cubeuy) titr 'PROFIL Uy sur axe'
tity 'Uy/U0' titx 'X/D';
'DESSIN' (eKMaxe 'ET' L0 'ET' cubeK ) 'TITRE' 'K moyen sur axe';
dess (eKaxe 'ET' L0 'ET' cubeK ) titr 'PROFIL K sur axe'
TAB2 tity 'K/U0**2' titx 'X/D';
* opti donn 5;
fin;
© Cast3M 2003 - All rights reserved.
Disclaimer