* fichier : carre_expl.dgibi
*****************************************************
** ECOULEMENT AUTOUR D'UNE CYLINDRE DE SECTION CAREE*
** Gregory Turbelin 29/12/1998
** modifie le 4/02/2000
**
** CAS TEST DE LYN ET RODI
** RESULTATS EXPERIMENTAUX REFERENCES A L_ADRESSE
** http://vortex.mech.surrey.ac.uk/database/test43/test43.html
*****************************************************
GRAPH='N' ;
COMPLET=FAUX ;
'OPTION' 'DIME' 2 'ELEM' qua8;
**************************
** PROCEDURES DE CALCUL
**************************
**--------------------------
***************************************
**Calcul des valeurs moyennes sur l'axe
***************************************
** En chaque point de l'axe, et a chaque iteration,
** on additionne les valeurs de la variable. La moyenne, en un point donne,
** sera la somme des valeurs, divisee par le nbre de valeurs stockees
'DEBPROC' VMPM;
'ARGUMENT' rvx*'TABLE';
rv= rvx.'EQEX';
DD = rv.PASDETPS.'NUPASDT' ;
NN = DD/20;
T0 = (DD '-' (20*NN)) 'EGA' 0;
'SI' T0;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
UN = RV.INCO.'UN' ;
unx= kcht $surfa scal sommet (exco 'UX' un) ;
** compteur pour connaitre le nombre de valeurs stockées
cont = 1.;
rv.INCO.'cont'= (rv.INCO.'cont' '+' cont);
** VITESSE LE LONG DE L'AXE
eUXa_12 = EVOL CHPO UNX axe_12;
UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
eUXa_34 = EVOL CHPO UNX 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 );
** PRESSION SUR L'OBSTACLE
ePobst = 'EVOL' chpo PP scal ('INVERSE' obstacle);
Pobst = 'EXTRAIRE' ePobst ordo;
rv.INCO.'Pobst'= (rv.INCO.'Pobst' '+' Pobst);
** ENERGIE 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' ;
as2 ama1 = 'KOPS' 'MATRIK' ;
'FINPROC' as2 ama1 ;
***************************************
**CALCUL DES COEFF AERODYNAMIQUES
**PAR INTEGRATION DU CHAMP DE PRESSION
**SUR LES PAROIS DE L_OBSTACLE
***************************************
**------------------------------------------------
'DEBPROC' CLCD;
'ARGUMENT' rvx*'TABLE';
rv = rvx.'EQEX';
DD = rv.PASDETPS.'NUPASDT' ;
NN = DD/20;
T0 = (DD '-' (20*NN)) 'EGA' 0;
'SI' T0;
PP = 'ELNO' ('KCHT' $surfa scal centre (RVP.PRESSION)) $surfa;
** PRESSION FACE1
PP1 = 'REDU' PP face1 ;
** PRESSION FACE2
PP2 = 'REDU' PP face2;
**PRESSION FACE3
PP3 = 'REDU' PP face3;
**PRESSION FACE4
PP4 = 'REDU' PP face4;
**COORD. face1
XX1 = REDU CHPX face1 ;
**COORD. face2
YY2 = REDU CHPY face2 ;
**COORD. face3
XX3 = REDU CHPX face3 ;
**COORD. face4
YY4 = REDU CHPY face4 ;
** CHPO ---> EVOL (CHPO)
ePP1 = EVOL CHPO PP1 face1 ;
ePP2 = EVOL CHPO PP2 face2 ;
ePP3 = EVOL CHPO PP3 face3 ;
ePP4 = EVOL CHPO PP4 face4 ;
eXX1 = EVOL CHPO XX1 face1 ;
eXX3 = EVOL CHPO XX3 face3 ;
eYY2 = EVOL CHPO YY2 face2 ;
eYY4 = EVOL CHPO YY4 face4 ;
** EVOL (CHPO) --> LISTREEL
P1 = EXTR ePP1 ORDO;
P2 = EXTR ePP2 ORDO;
P3 = EXTR ePP3 ORDO;
P4 = EXTR ePP4 ORDO;
X1 = EXTR eXX1 ORDO;
X3 = EXTR eXX3 ORDO;
Y2 = EXTR eYY2 ORDO;
Y4 = EXTR eYY4 ORDO;
** CALCUL DE LA PORTANCE ET TRAINEE
*
** LISTREEL --> EVOL (MANU)
evXP1 = EVOL MANU X1 P1;
evXP3 = EVOL MANU X3 P3;
evyP2 = EVOL MANU Y2 P2 ;
evYP4 = EVOL MANU Y4 P4 ;
** CALCUL INTEGRALE (Pdx)
intXP1 = SOMME evXP1 ;
intXP3 = SOMME evXP3 ;
** CALCUL INTEGRALE (Pdy)
intYP2 = SOMME evYP2 ;
intYP4 = SOMME evYP4 ;
CL = ((-1)*intXP3)+ intXP1 ;
CD = (intYP4)+((-1)*intYP2) ;
rv.inco.'CL' = (rv.inco.'CL') et (CL);
rv.inco.'CD' = (rv.inco.'CD') et (CD);
PT = 'PROG' rv.'PASDETPS'.'TPS' ;
rv.INCO.'TP'= rv.INCO.'TP' ET PT;
'FINSI' ;
as2 ama1 = 'KOPS' 'MATRIK' ;
'FINPROC' as2 ama1 ;
****************************************
** PROCEDURE 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 ;
**-------------------------------------------
** Procedures pour estimer la valeur moyenne
** et l_ecart type 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)
*************************************
**--------------------------------
*** Diametre du cylindre (longueur de reference)
*** D_ref = 0.04 m <==> X2 (= Y2+Y3)
**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;
*
********************************************
**LIGNES POUR UN PARAMETRAGE ZONES PAR ZONES
**--------------------------------
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;
surfa = 'CHANGER' surfa 'QUAF' ;
bas = 'CHANGER' bas 'QUAF' ;
entree = 'CHANGER' entree 'QUAF' ;
haut = 'CHANGER' haut 'QUAF' ;
sortie = 'CHANGER' sortie 'QUAF' ;
face1 = 'CHANGER' face1 'QUAF' ;
face2 = 'CHANGER' face2 'QUAF' ;
face3 = 'CHANGER' face3 'QUAF' ;
face4 = 'CHANGER' face4 'QUAF' ;
axe_12 = 'CHANGER' axe_12 'QUAF' ;
axe_34 = 'CHANGER' axe_34 'QUAF' ;
'ELIMINATION' 1.D-6 (bas 'ET' entree 'ET' haut 'ET' sortie
et face1 et face2 et face3 et face4 et axe_12 et axe_34
'ET' surfa) ;
*************************
** 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 ;
$face1 = 'MODE' face1 'NAVIER_STOKES' MACRO ;
$face2 = 'MODE' face2 'NAVIER_STOKES' MACRO ;
$face3 = 'MODE' face3 'NAVIER_STOKES' MACRO ;
$face4 = 'MODE' face4 'NAVIER_STOKES' MACRO ;
$axe_12 = 'MODE' axe_12 'NAVIER_STOKES' MACRO ;
$axe_34 = 'MODE' axe_34 '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' ;
face1 = 'DOMA' $face1 'MAILLAGE' ;
face2 = 'DOMA' $face2 'MAILLAGE' ;
face3 = 'DOMA' $face3 'MAILLAGE' ;
face4 = 'DOMA' $face4 'MAILLAGE' ;
axe_34 = 'DOMA' $axe_34 'MAILLAGE' ;
axe_12 = 'DOMA' $axe_12 'MAILLAGE' ;
si ('EGA' graph 'O' );
'TRACER' surfa;
finsi ;
** ABSC ET ORDO des points du maillage
CHPX CHPY = COOR SURFA;
*****************************
*** PROPRIETES PHYSIQUES
***************************
** Hauteur du domaine (sans dim)
Hda = (Y1 + Y2 + Y3 + Y4);
** Longueur du domaine (sans dim)
Lda = (X1+X2+X3) ;
** Reynolds de l'écoulement
Re = 2.2e4;
iRe = 1/Re ;
** Distances Normales aux parois de l'obstacle
** exprimees en fonction de la taille des mailles
** sur la face_1
Xface1 = coor 1 face1;
eXface1 = evol chpo Xface1 face1;
eXface1 = 'EXTRAIRE' eXface1 ABSC;
** Yp = 1/5 hauteur de la premiere maille (sans dim)
YP = ('EXTR' eXface1 2)/5;
** vitesse de frottement sans dim
ufini = 150./(yp*Re);
**constante de Von Karman
C1 = 0.4;
** constante du modele RNG k-e
CNU = 0.0845;
** taux de turbulence
INT = 0.02 ;
alpha = 1.5 ;
** Energie cinetique
** Kadim = K/(U_ref**2)
Kk = alpha*((INT)**2) ;
** Taille des tourbillons energetiques
** On suppose que k**0.5/Lk = U0/L0
** ou L0 est une dimension caracteristique
** associe a l_ecoulement moyen
L0 = 1. ;
Lk = (alpha**0.5) * INT * L0;
** Dissipation
Ek = (Kk**1.5)'/'(Lk);
**Reynolds associe aux tourbillons energetiques
Rek = alpha *((INT)**2) * (Re*L0) ;
**NUT/NU
NUTSNU = cnu * ((Kk**2)/Ek) * Re ;
** On doit, normalement, avoir NUTSNU = cnu *Rek
**Verification
Liste Lda ;
List Hda ;
'LISTE' Re;
'LISTE' nu;
'LISTE' yp;
'LISTE' ufini;
list lk ;
list kk ;
list Ek ;
list Rek ;
list (CNU*Rek) ;
list NUTSNU;
** 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
**-----------------
K1 = 'MANU' 'CHPO' surfa 1 scal Kk;
***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
**--------------
E1 = 'MANU' 'CHPO' surfa 1 scal Ek;
** 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 Ufini;
UET_ini2 = 'KCHT' $face2 scal centre Ufini;
UET_ini3 = 'KCHT' $face3 scal centre Ufini;
UET_ini4 = 'KCHT' $face4 scal centre Ufini;
** 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' YP INCO 'UN' 'KN' 'EN'
ZONE $face2 OPER FPU iRe 'UET2' YP INCO 'UN' 'KN' 'EN'
ZONE $face3 OPER FPU iRe 'UET3' YP INCO 'UN' 'KN' 'EN'
ZONE $face1 OPER FPU iRe 'UET1' YP INCO 'UN' '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
ZONE $surfa OPER VMPM
ZONE $surfa OPER CLCD ;
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' ;
RV.INCO.'CD' = 'PROG' ;
RV.INCO.'TP' = 'PROG' ;
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. ;
rv.inco.'IT' = PROG 1 ;
RV.INCO.'UNM1'= KCHT $surfa vect sommet (1.E-3 1.E-3) ;
rv.inco.'ERY' = PROG 0. ;
rv.inco.'ERX' = PROG 0. ;
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;
** 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
**---------------------
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 = RV.INCO.'TP' ;
Ntp = ('DIME' tp);
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 TITR 'DSP';
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;
© Cast3M 2003 - All rights reserved.
Disclaimer