Télécharger la plaquette Cast3M


fr en

Télécharger carre.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  carre.dgibi
  2. *****************************************************
  3. ** ECOULEMENT AUTOUR D'UNE CYLINDRE DE SECTION CAREE*
  4. ** Gregory Turbelin 29/12/1998                      *
  5. ** modifie le 15/06/2014 passage EQPR -> EQEX       *
  6. *****************************************************
  7.  
  8. GRAPH='N' ;
  9. COMPLET=FAUX ;
  10. 'OPTION'  'DIME'  2 'ELEM'  qua8;
  11.  
  12. DISCR  = 'MACRO' ;
  13. KPRESS = 'CENTRE';
  14. BETAP  = 1.      ;
  15.  
  16.  
  17. ** PROCEDURES DE CALCUL
  18. **--------------------
  19. 'DEBPROC' VMPM;
  20.  'ARGUMENT' rvx*'TABLE';
  21.   rv= rvx.'EQEX';
  22.  
  23.   DD = rv.PASDETPS.'NUPASDT' ;
  24.   NN = DD/20;
  25.   L0 = (DD '-' (20*NN)) 'EGA' 0;
  26.   'SI' L0;
  27.   PP = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  28. ** compteur pour connaitre le nombre de valeurs stockées
  29.  cont = 1.;
  30.  rv.INCO.'cont'= (rv.INCO.'cont' '+' cont);
  31.  
  32. ** CALCUL DE LA VITESSE MOYENNE LE LONG DE L'AXE
  33.  eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
  34.  UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
  35.  
  36.  eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
  37.  UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
  38.  
  39.  rv.INCO.'UXa_12'= (rv.INCO.'UXa_12' '+' UXa_12 );
  40.  rv.INCO.'UXa_34'= (rv.INCO.'UXa_34' '+' UXa_34 );
  41.  
  42. ** CALCUL DE LA PRESSION MOYENNE SUR L'OBSTACLE
  43.  
  44.  ePobst = 'EVOL' chpo PP scal ('INVERSE' obstacle);
  45.  Pobst  = 'EXTRAIRE' ePobst ordo;
  46.  rv.INCO.'Pobst'= (rv.INCO.'Pobst' '+' Pobst);
  47.  
  48. ** CALCUL DE L ENERGIE MOYENNE LE LONG DE L'AXE
  49.  eka_12 = EVOL CHPO (rv.inco.'KN') scal axe_12;
  50.  ka_12 = 'EXTRAIRE' eKa_12 ORDO;
  51.  
  52.  eKa_34 = EVOL CHPO (rv.inco.'KN') scal axe_34;
  53.  ka_34 = 'EXTRAIRE' eka_34 ordo;
  54.  
  55.  rv.INCO.'Ka_12'= (rv.INCO.'Ka_12' '+' Ka_12 );                                      
  56.  rv.INCO.'Ka_34'= (rv.INCO.'Ka_34' '+' Ka_34 );
  57.  
  58.  'FINSI' ;
  59.  as2 ama1  = 'KOPS' 'MATRIK' ;
  60.  'FINPROC' as2 ama1 ;
  61.  
  62. **------------------------------------------------ 
  63.  'DEBPROC' CLCD;
  64.  'ARGUMENT' rvx*'TABLE';
  65.   rv = rvx.'EQEX';
  66.  
  67.   DD = rv.PASDETPS.'NUPASDT' ;
  68.   NN = DD/20;
  69.   L0 = (DD '-' (20*NN)) 'EGA' 0;
  70.   'SI' L0;
  71.  
  72.   PP = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  73.  
  74. ** CALCUL DES COEFF AERODYNAMIQUES
  75.  
  76.   PPSF1 = 'REDU' PP face1; 
  77.   PPCf1 = NOEL $face1 ('KCHT' $face1 SCAL SOMMET PPSF1);
  78.   PPSF2 = 'REDU' PP face2; 
  79.   PPCf2 = NOEL $face2 ('KCHT' $face2 SCAL SOMMET PPSF2);
  80.   PPSF3 = 'REDU' PP face3;
  81.   PPCf3 = NOEL $face3 ('KCHT' $face3 SCAL SOMMET PPSF3);
  82.   PPSF4 = 'REDU' PP face4;
  83.   PPCf4 = NOEL $face4 ('KCHT' $face4 SCAL SOMMET PPSF4);
  84.  
  85.     i=1;
  86.     CD1 = 0.;
  87.     CD3 = 0.;
  88.  
  89.     'REPETER' bloc1 (NBEL (doma $face1 centre)) ;
  90.  
  91.     PPf1i = 'EXTRAIRE' PPCf1 SCAL (doma $face1 centre 'POIN' i) ;
  92.     Lf1i = 'ABS' ( (coor 1 (face1 'POIN' (i+1))) '-'
  93.     (coor 1 (face1 'POIN'  i) ) );
  94.     CD1i= PPf1i * Lf1i;
  95.     CD1 = CD1 '+'CD1i;
  96. *
  97.     PPf3i = 'EXTRAIRE' PPCf3 SCAL (doma $face3 centre POIN i);
  98.     Lf3i = 'ABS' ( (coor 1 (face3 'POIN' i)) '-'
  99.     (coor 1 (face3 'POIN' (i+1) )) );
  100.     CD3i= PPf3i * Lf3i;
  101.     CD3 = CD3 '+' CD3i;
  102.     i = i+1;
  103.    'FIN' bloc1;
  104.  
  105.    CL = CD1 '-' CD3; 
  106. *
  107.     j=1;
  108.     CD2 = 0.;
  109.     CD4 = 0.;
  110.    'REPETER' bloc2 (NBEL (doma $face2 centre));
  111.     PPf2j = 'EXTRAIRE' PPCf2 SCAL (doma $face2 centre POIN j);
  112.     Lf2j = 'ABS' ( (coor 2 (face2 'POIN' j)) '-'  
  113.            (coor 2 (face2 'POIN' (j+1) )) );
  114.     CD2j= PPf2j *  Lf2j;
  115.     CD2 = CD2 '+' CD2j;
  116. *
  117.     PPf4j = 'EXTRAIRE' PPCf4  SCAL (doma $face4 centre POIN j);
  118.     Lf4j = 'ABS' ( (coor 2 (face4 POIN j)) - 
  119.            (coor 2 (face4 POIN (j+1) )) );
  120.     CD4j= PPf4j * Lf4j;
  121.     CD4 = CD4 '+' CD4j;
  122.     j = j+1;
  123.    'FIN' bloc2;
  124. *
  125.    CD = CD4 '-' CD2 ; 
  126.  
  127.    rv.inco.'CL'  = (rv.inco.'CL') et (Prog CL);
  128.    rv.inco.'CD'  = (rv.inco.'CD') et (Prog CD);
  129.  
  130.    PT = 'PROG' rv.'PASDETPS'.'TPS' ;
  131.    rv.INCO.'TP'= rv.INCO.'TP' ET PT;
  132.  
  133.   'FINSI' ;
  134.  as2 ama1  = 'KOPS' 'MATRIK' ;
  135.  'FINPROC' as2 ama1 ;
  136.  
  137. ** Procedure pour determiner la valeur moyenne
  138. ** d'une liste de relles
  139.  'DEBPROC' MOYL LR*'LISTREEL';
  140.    LSi =0.;
  141.    i=1;
  142.   'REPETER' boucle  ('DIME' LR);
  143.      ielem = 'EXTRAIRE' LR i;
  144.      LSi = LSi '+' ielem;
  145.      moyenne  =  LSi/i;
  146.    i =i+1;
  147.   'FIN'  boucle;
  148.  'FINPROC' moyenne ;
  149. ** Ecart 'TYPE'
  150.  'DEBPROC' RMS LR*'LISTREEL';
  151.    LMi = 0.;
  152.    LEi = 0.;
  153.    i=1;
  154.   'REPETER' boucle  ('DIME' LR);
  155.      ielem = 'EXTRAIRE' LR i;
  156.      ielem2 = (ielem**2);
  157.  
  158.      LMi = LMi '+' ielem;
  159.      LEi = LEi '+' ielem2;
  160.  
  161.      moyenne  =  LMi/i;
  162.      moyenne2 =  LEi/i;  
  163.    i =i+1;
  164.   'FIN'  boucle;
  165.   ecart = (moyenne2 '-' ((moyenne)**2))**(0.5); 
  166.  'FINPROC' ecart ;
  167.  
  168. ** CONSTRUCTION DU MAILLAGE (ADIM)
  169. **--------------------------------
  170. *
  171. **Longueurs sans dim
  172. *
  173. X1 = 4.5 ; X2 = 1.  ; X3 = 14.5;
  174. Y1 = 6.5 ; Y2 = 0.5 ; Y3 = 0.5 ; Y4=6.5 ;
  175. *
  176. **Points caracteristiques
  177. *
  178. P11=0.  0.;
  179. P12=X1  0.;
  180. P13=(X1+X2)  0.;
  181. P14=(X1+X2+X3)  0.;
  182. *
  183. P21=0.   Y1;
  184. P22=X1  Y1;
  185. P23=(X1+X2)  Y1;
  186. P24=(X1+X2+X3) Y1;
  187. *
  188. P31=0.  (Y1+Y2);
  189. P32=X1  (Y1+Y2);
  190. P33=(X1+X2) (Y1+Y2);
  191. P34=(X1+X2+X3)  (Y1+Y2);
  192. *
  193. P41=0.  (Y1+Y2+Y3);
  194. P42=X1  (Y1+Y2+Y3);
  195. P43=(X1+X2)  (Y1+Y2+Y3);
  196. P44=(X1+X2+X3) (Y1+Y2+Y3);
  197.  
  198. P51=0.  (Y1+Y2+Y3+Y4);
  199. P52=X1  (Y1+Y2+Y3+Y4);
  200. P53=(X1+X2)  (Y1+Y2+Y3+Y4);
  201. P54=(X1+X2+X3) (Y1+Y2+Y3+Y4);
  202. *
  203. **Nb de mailles horizontales, sur obstacle
  204. si complet ;
  205. Nx = 6;
  206. sinon ;
  207. Nx = 2;
  208. finsi ;
  209. **Nb de mailles verticales 1/2 obstacle
  210. si complet ;
  211. Ny = 3; 
  212. sinon ;
  213. Ny = 1;
  214. finsi ;
  215. **Nb de mailles horizontales, avant obstacle
  216. Nx1= 3*Nx;
  217. **Nb de mailles horizontales, sur obstacle
  218. Nx2 = Nx ;
  219. **Nb de mailles horizontales, après obstacle
  220. Nx3 =7*Nx;
  221. *
  222. **Nb de mailles verticales sous obstacle 
  223. Ny1 = 8*Ny;
  224. **Nb de mailles verticales 1/2 obstacle inf
  225. Ny2 = Ny;
  226. **Nb de mailles verticales 1/2 obstacle sup
  227. Ny3 = Ny;
  228. **Nb de mailles verticales au dessus obstacle 
  229. Ny4 = 8*Ny;
  230. *
  231. *
  232. ** Definition des densites
  233. **-------------------------
  234. DX1 = (X1/Nx1);
  235. Dx1_1 = Dx1 ; Dx1_2 = 0.3*Dx1;
  236. *
  237. DX2 = (X2/Nx2);
  238. Dx2_2 = Dx2 ; Dx2_3 = Dx2;
  239. *
  240. DX3 = (X3/Nx3);
  241. Dx3_3 = 0.3*Dx3 ; Dx3_4 = 1.5*Dx3;
  242. *
  243. Dy1 = (Y1/NY1);
  244. Dy1_1 =1.*DY1 ; Dy1_2 = 0.3*DY1;
  245. *
  246. DY2 = (Y2/NY2);
  247. DY2_2 = DY2; DY2_3 = DY2;
  248. *
  249. DY3 = (Y3/NY3);
  250. DY3_3 =0.5*DY3 ;  DY3_4 = DY3;
  251. *
  252. DY4 = (Y4/NY4);
  253. DY4_4 =0.3*DY4 ;  DY4_5 = 1.*DY4;
  254. *
  255. **Definition des Lignes
  256. **--------------------------------
  257. L1x_12 = P11 droite (-1*Nx1) P12 'DINI' Dx1_1 'DFIN' Dx1_2; 
  258. L1x_23 = P12 droite (-1*Nx2) P13 'DINI' Dx2_2 'DFIN' Dx2_3; 
  259. L1x_34 = P13 droite (-1*Nx3) P14 'DINI' Dx3_3 'DFIN' Dx3_4;
  260.  
  261. L2x_12 = P21 droite (-1*Nx1) P22 'DINI' Dx1_1 'DFIN' Dx1_2; 
  262. L2x_23 = P22 droite (-1*Nx2) P23 'DINI' Dx2_2 'DFIN' Dx2_3; 
  263. L2x_34 = P23 droite (-1*Nx3) P24 'DINI' Dx3_3 'DFIN' Dx3_4;
  264.  
  265. L3x_12 = P31 droite (-1*Nx1) P32 'DINI' Dx1_1 'DFIN' Dx1_2; 
  266. L3x_23 = P32 droite (-1*Nx2) P33 'DINI' Dx2_2 'DFIN' Dx2_3; 
  267. L3x_34 = P33 droite (-1*Nx3) P34 'DINI' Dx3_3 'DFIN' Dx3_4;
  268.  
  269. L4x_12 = P41 droite (-1*Nx1) P42 'DINI' Dx1_1 'DFIN' Dx1_2; 
  270. L4x_23 = P42 droite (-1*Nx2) P43 'DINI' Dx2_2 'DFIN' Dx2_3; 
  271. L4x_34 = P43 droite (-1*Nx3) P44 'DINI' Dx3_3 'DFIN' Dx3_4;
  272.  
  273. L5x_12 = P51 droite (-1*Nx1) P52 'DINI' Dx1_1 'DFIN' Dx1_2; 
  274. L5x_23 = P52 droite (-1*Nx2) P53 'DINI' Dx2_2 'DFIN' Dx2_3; 
  275. L5x_34 = P53 droite (-1*Nx3) P54 'DINI' Dx3_3 'DFIN' Dx3_4;
  276.  
  277. L1x =  L1x_12 'ET'  L1x_23 'ET'  L1x_34;
  278. L2x =  L2x_12 'ET'  L2x_23 'ET'  L2x_34;
  279. L3x =  L3x_12 'ET'  L3x_23 'ET'  L3x_34;
  280. L4x =  L4x_12 'ET'  L4x_23 'ET'  L4x_34;
  281. L5x =  L5x_12 'ET'  L5x_23 'ET'  L5x_34; 
  282.  
  283. L1y_12 = P11 droite (-1*Ny1) P21 'DINI' Dy1_1 'DFIN' Dy1_2; 
  284. L1y_23 = P21 droite (-1*Ny2) P31 'DINI' Dy2_2 'DFIN' Dy2_3; 
  285. L1y_34 = P31 droite (-1*Ny3) P41 'DINI' Dy3_3 'DFIN' Dy3_4;
  286. L1y_45 = P41 droite (-1*Ny4) P51 'DINI' Dy4_4 'DFIN' Dy4_5;
  287.  
  288. L2y_12 = P12 droite (-1*Ny1) P22 'DINI' Dy1_1 'DFIN' Dy1_2; 
  289. L2y_23 = P22 droite (-1*Ny2) P32 'DINI' Dy2_2 'DFIN' Dy2_3; 
  290. L2y_34 = P32 droite (-1*Ny3) P42 'DINI' Dy3_3 'DFIN' Dy3_4;
  291. L2y_45 = P42 droite (-1*Ny4) P52 'DINI' Dy4_4 'DFIN' Dy4_5;
  292.  
  293. L3y_12 = P13 droite (-1*Ny1) P23 'DINI' Dy1_1 'DFIN' Dy1_2; 
  294. L3y_23 = P23 droite (-1*Ny2) P33 'DINI' Dy2_2 'DFIN' Dy2_3; 
  295. L3y_34 = P33 droite (-1*Ny3) P43 'DINI' Dy3_3 'DFIN' Dy3_4;
  296. L3y_45 = P43 droite (-1*Ny4) P53 'DINI' Dy4_4 'DFIN' Dy4_5;
  297.  
  298. L4y_12 = P14 droite (-1*Ny1) P24 'DINI' Dy1_1 'DFIN' Dy1_2; 
  299. L4y_23 = P24 droite (-1*Ny2) P34 'DINI' Dy2_2 'DFIN' Dy2_3; 
  300. L4y_34 = P34 droite (-1*Ny3) P44 'DINI' Dy3_3 'DFIN' Dy3_4;
  301. L4y_45 = P44 droite (-1*Ny4) P54 'DINI' Dy4_4 'DFIN' Dy4_5;
  302.  
  303. L1y =  L1y_12 'ET'  L1y_23 'ET'  L1y_34 'ET'  L1y_45;
  304. L2y =  L2y_12 'ET'  L2y_23 'ET'  L2y_34 'ET'  L2y_45;
  305. L3y =  L3y_12 'ET'  L3y_23 'ET'  L3y_34 'ET'  L3y_45; 
  306. L4y =  L4y_12 'ET'  L4y_23 'ET'  L4y_34 'ET'  L4y_45;  
  307.  
  308. bas    = (L1x) 'COULEUR' vert ; 'ELIMINATION' bas 1.e-5;
  309. sortie = (L4y) 'COULEUR' jaune; 'ELIMINATION' sortie 1.e-5;
  310. haut   = (L5x) 'COULEUR' bleu ; 'ELIMINATION' haut 1.e-5;
  311. entree = (L1y) 'COULEUR' rouge; 'ELIMINATION' entree 1.e-5;
  312. axe_12 = L3x_12; axe_34=L3x_34;
  313. axe    = (axe_12 'ET' axe_34) 'COULEUR' turq; 'ELIMINATION' axe 1.e-5;
  314.  
  315. face1 = L2x_23 'COULEUR' vert ;
  316. face2 = (L3y_23 'ET' L3y_34) 'COULEUR' jaune;
  317. face3 = L4x_23 'COULEUR' bleu ;
  318. face4 = (L2y_23 'ET' L2y_34) 'COULEUR' rouge ;
  319.  
  320. obstacle = FACE1 'ET' FACE2 'ET' FACE3 'ET' FACE4 ;
  321. 'ELIMINATION' obstacle 1.e-5; ;
  322.  
  323. **Creation de sous-zones
  324. ZX1Y1 = 'DALLER' L1x_12 L2y_12 L2x_12  L1y_12  PLAN ;
  325. ZX1Y2 = 'DALLER' L2x_12 L2y_23 L3x_12  L1y_23  PLAN ;
  326. ZX1Y3 = 'DALLER' L3x_12 L2y_34 L4x_12  L1y_34  PLAN ;
  327. ZX1Y4 = 'DALLER' L4x_12 L2y_45 L5x_12  L1y_45  PLAN ;
  328.  
  329. ZX1 = ZX1Y1 'ET' ZX1Y2 'ET' ZX1Y3 'ET' ZX1Y4;
  330. 'ELIMINATION' 1.e-5  ZX1 ;
  331.  
  332. ZX2Y1 = 'DALLER' L1x_23 L3y_12 L2x_23  L2y_12  PLAN ;
  333. ZX2Y4 = 'DALLER'  L4x_23 L3y_45 L5x_23  L2y_45  PLAN ;
  334.  
  335. ZX2 = ZX2Y1 'ET' ZX2Y4;
  336. 'ELIMINATION'  ZX2 1.e-5;
  337. orienter Zx2;
  338.  
  339. ZX3Y1 = 'DALLER' L1x_34 L4y_12 L2x_34  L3y_12  PLAN ;
  340. ZX3Y2 = 'DALLER' L2x_34 L4y_23 L3x_34  L3y_23  PLAN ;
  341. ZX3Y3 = 'DALLER' L3x_34 L4y_34 L4x_34  L3y_34  PLAN ;
  342. ZX3Y4 = 'DALLER' L4x_34 L4y_45 L5x_34  L3y_45  PLAN ;
  343.  
  344. ZX3 = ZX3Y1 'ET' ZX3Y2 'ET' ZX3Y3 'ET' ZX3Y4;
  345. 'ELIMINATION'  ZX3 1.e-5;
  346. orienter Zx3;
  347.  
  348. SURFA = Zx1 'ET' ZX2 'ET' Zx3 ;
  349. 'ELIMINATION' surfa 1.e-5;
  350. orienter surfa;
  351.  
  352. *'TRACER' surfa;
  353. *'OPTION' donn 5;
  354.  
  355. ** Définition des domaines
  356. **------------------------
  357.  
  358.  Msurfa  = 'CHAN' surfa  'QUAF' ;
  359.  Mbas    = 'CHAN' bas    'QUAF' ;
  360.  Mentree = 'CHAN' entree 'QUAF' ;
  361.  Mhaut   = 'CHAN' haut   'QUAF' ;
  362.  Msortie = 'CHAN' sortie 'QUAF' ;
  363.  Mface1  = 'CHAN' face1  'QUAF' ;
  364.  Mface2  = 'CHAN' face2  'QUAF' ;
  365.  Mface3  = 'CHAN' face3  'QUAF' ;
  366.  Mface4  = 'CHAN' face4  'QUAF' ;
  367.  Maxe_12 = 'CHAN' axe_12 'QUAF' ;                                
  368.  Maxe_34 = 'CHAN' axe_34 'QUAF' ;
  369.  
  370.  'ELIM' (Msurfa et Mbas et Mentree et Mhaut et Msortie et
  371.  Mface1 et Mface2 et Mface3 et Mface4 et Maxe_12 et Maxe_34)1.e-5;
  372.  
  373.  
  374.  $surfa  = 'MODE' Msurfa  'NAVIER_STOKES' DISCR ;
  375.  $bas    = 'MODE' Mbas    'NAVIER_STOKES' DISCR ;
  376.  $entree = 'MODE' Mentree 'NAVIER_STOKES' DISCR ;
  377.  $haut   = 'MODE' Mhaut   'NAVIER_STOKES' DISCR ;
  378.  $sortie = 'MODE' Msortie 'NAVIER_STOKES' DISCR ;
  379.  $face1  = 'MODE' Mface1  'NAVIER_STOKES' DISCR ;
  380.  $face2  = 'MODE' Mface2  'NAVIER_STOKES' DISCR ;
  381.  $face3  = 'MODE' Mface3  'NAVIER_STOKES' DISCR ;
  382.  $face4  = 'MODE' Mface4  'NAVIER_STOKES' DISCR ;
  383.  $axe_12 = 'MODE' Maxe_12 'NAVIER_STOKES' DISCR ;
  384.  $axe_34 = 'MODE' Maxe_34 'NAVIER_STOKES' DISCR ;
  385.  
  386.  surfa   = 'DOMA' $surfa   maillage ;
  387.  'ORIENTER' surfa;
  388.  
  389.  entree  = 'DOMA'  $entree maillage ;
  390.  haut    = 'DOMA'  $haut   maillage ;
  391.  bas     = 'DOMA'  $bas    maillage ;
  392.  sortie  = 'DOMA'  $sortie maillage ;
  393.  face1   = 'DOMA'  $face1  maillage ;
  394.  face2   = 'DOMA'  $face2  maillage ;
  395.  face3   = 'DOMA'  $face3  maillage ;
  396.  face4   = 'DOMA'  $face4  maillage ;
  397.  axe_34  = 'DOMA'  $axe_34 maillage ;
  398.  axe_12  = 'DOMA'  $axe_12 maillage ;
  399.  
  400. ** Proprietes physiques
  401. **-------------------
  402.   nu_eau = 1.e-6;
  403.   nu = nu_eau;
  404.  
  405. ** Longeur de reference = Diametre du cylindre en m
  406.  D_ref = 0.04;
  407. ** Hauteur du domaine (sans dim)
  408.  Hda = 14.;
  409. ** Longueur du domaine (sans dim)  
  410.  Lda = 20.;
  411. ** Reynolds de l'écoulement
  412.  Re = 2.2e4; 
  413.  iRe = 1/Re ;
  414. ** Vitesse de reference = vitesse en entree (m/s)
  415.  U_ref = (Re*nu)/D_ref;  
  416.  
  417. ** Distances Normales aux  parois de l'obstacle
  418. **---------------------------------------------
  419. ** Coordonnees des points des faces  du cylindre
  420.  Xface1 = coor 1  face1;
  421.  Yface2 = coor 2  face2;
  422.  Xface3 = coor 1  face3;
  423.  Yface4 = coor 2  face4;
  424.  
  425. ** abscisse curviligne le long des faces
  426.  eXface1 = evol chpo Xface1 face1;
  427.  eXface1 = 'EXTRAIRE' eXface1 ABSC;
  428.  eYface2 = evol chpo Yface2 face2;
  429.  eYface2 = 'EXTRAIRE' eYface2 ABSC;
  430.  eXface3 = evol chpo Xface3 face3;
  431.  eXface3 = 'EXTRAIRE' eXface3 ABSC;
  432.  eYface4 = evol chpo Yface4 face4;
  433.  eYface4 = 'EXTRAIRE'  eYface4 ABSC;
  434.  
  435. ** Yp = 1/5 hauteur de la premiere maille  (sans dim)
  436.  YP1 = ('EXTR' eXface1 2)/5;
  437.  YP2 = ('EXTR' eYface2 2)/5;
  438.  YP3 = ('EXTR' eXface3 2)/5;
  439.  YP4 = ('EXTR' eYface4 2)/5;
  440. ** vitesse de frottement sans dim
  441.  uf1 = (150*nu)/yp1;
  442.  uf2 = (150*nu)/yp2;
  443.  uf3 = (150*nu)/yp3;
  444.  uf4 = (150*nu)/yp4;
  445. **constante de Von Karman
  446.  C1  = 0.4;
  447. ** constante du modele RNG  k-e
  448.  CNU = 0.0845;
  449. ** taux de turbulence 
  450.  INT =  0.02 ;
  451.  
  452. ** Definition des profils d'entrees et initiaux
  453. **----------------------------------------------
  454. ** Definition du domaine d'imposition
  455.  M_imp =entree;
  456.  
  457. ** VITESSES
  458. **-----------
  459.  U1 = 1. 0.;
  460.  U1 = 'KCHT' $surfa vect sommet  U1;
  461.  
  462. ** Valeurs IMPOSEES
  463.  UX_IMP = 'EXCO' ('REDU' U1 M_imp) 'UX' SCAL;
  464.  UY_IMP = 'EXCO' ('REDU' U1 M_imp) 'UY' SCAL;
  465.  
  466. ** Valeurs INITIALES (sans dim)
  467.  U_ini = 'MANU' 'CHPO' surfa 1 scal 0.;
  468.  U_ini = U_ini 'NOMC' 'UX' ;
  469.  U_ini = 'KCHT' $surfa vect sommet U_ini;
  470.  
  471. **ENERGIE CINETIQUE
  472. **-----------------
  473. ** Kadim = K/(U_ref**2) 
  474. ** En fonction de l'intensite de la turbulence
  475.  
  476.  K1 = 1.5*(INT)**2 ;
  477.  K1 = 'MANU' 'CHPO' surfa 1 scal K1;
  478. ***INITIALE
  479.  K_ini = 'KCHT' $surfa scal sommet K1;
  480. ***IMPOSEE
  481.  K_IMP = 'KCHT' $surfa scal sommet K1;
  482.  K_IMP= 'EXCO' ('REDU' K_IMP M_imp) SCAL;
  483.  
  484. **DISSIPATION  Eadim = E*Href/Uref**3 
  485. **--------------
  486. ** Fonction de K et de D
  487. ** ! K  est sans dim et D = Href
  488.  E1  = (K1**1.5);
  489. ** INITIALE
  490.  E_ini = 'KCHT' $surfa scal sommet E1;
  491. **IMPOSEE
  492.  E_IMP = 'KCHT' $surfa scal sommet E1;
  493.  E_IMP='EXCO' ('REDU' E_IMP M_imp) SCAL;
  494.  
  495. **NUT INITIALE
  496. **------------------
  497. *fonction de K et E
  498.  N1=CNU*(K_ini**2)*(E_ini**-1);
  499.  N1 = 'KCHT' $surfa scal sommet N1;
  500.  N_ini = 'NOEL'  $surfa N1 ;
  501.  
  502. ** U* INITIALES 
  503. **-------------
  504.  UET_ini1 = 'KCHT' $face1 scal centre Uf1;
  505.  UET_ini2 = 'KCHT' $face2 scal centre Uf2;
  506.  UET_ini3 = 'KCHT' $face3 scal centre Uf3;
  507.  UET_ini4 = 'KCHT' $face4 scal centre Uf4;
  508.  
  509. ** Module de resolution
  510. **-------------------
  511.  
  512.  RV = 'EQEX' $surfa ITMA 42 ALFA 0.8
  513.      OPTI 'RNG'
  514.      ZONE $surfa OPER FILTREKE 1. Hda iRe 'UN' INCO 'KN' 'EN'
  515.      ZONE $surfa OPER NSKE iRe 'NUT' INCO 'UN' 'KN' 'EN'
  516.      ZONE $face4 OPER FPU iRe 'UET4' YP4 INCO 'UN' 'KN' 'EN'
  517.      ZONE $face2 OPER FPU iRe 'UET2' YP2 INCO 'UN' 'KN' 'EN'
  518.      ZONE $face3 OPER FPU iRe 'UET3' YP3 INCO 'UN' 'KN' 'EN'
  519.      ZONE $face1 OPER FPU iRe 'UET1' YP1 INCO 'UN' 'KN' 'EN' ;
  520.  
  521.   RV=EQEX RV OPTI EFM1 'CENTREE'
  522.   'ZONE' $surfa   'OPER' 'DFDT' 1. 'UN' 'DELTAT' 'INCO' 'UN'
  523.   'ZONE' $surfa   'OPER' 'DFDT' 1. 'KN' 'DELTAT' 'INCO' 'KN'
  524.   'ZONE' $surfa   'OPER' 'DFDT' 1. 'EN' 'DELTAT' 'INCO' 'EN'
  525.   ;
  526.  
  527.  RV =  'EQEX' RV  
  528.     ZONE $surfa OPER CLCD
  529.     ZONE $surfa OPER VMPM;
  530.  
  531.  RV =  'EQEX' RV
  532.      CLIM 'UN' UIMP M_imp UX_imp
  533.      CLIM 'UN' VIMP M_imp UY_imp
  534.      CLIM 'UN' VIMP (haut 'ET' bas) 0.
  535.      CLIM 'KN' TIMP M_imp K_imp
  536.      CLIM 'EN' TIMP M_imp E_imp ;
  537.  
  538. RVP = EQEX 'OPTI' 'EF' KPRESS
  539.     'ZONE' $surfa OPER  KBBT  -1. betap 'INCO' 'UN' 'PRES'
  540.     'ZONE' $face4 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  541.     'ZONE' $face2 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  542.     'ZONE' $face3 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  543.     'ZONE' $face1 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  544.     ;
  545.  
  546.  RV.'PROJ' = RVP     ;
  547.  
  548.  RV.INCO = TABLE INCO;
  549.  RV.INCO.'PRES' = kcht $surfa SCAL KPRESS 0.;
  550.  RV.INCO.'UN'  = U_ini;
  551.  RV.INCO.'KN'  = K_ini;
  552.  RV.INCO.'EN'  = E_ini;
  553.  RV.INCO.'NUT' = N_ini;
  554.  RV.INCO.'UET4' = UET_ini4;
  555.  RV.INCO.'UET3' = UET_ini3;
  556.  RV.INCO.'UET2' = UET_ini2;
  557.  RV.INCO.'UET1' = UET_ini1;
  558.  RV.INCO.'CL'  = 'PROG' 0.;                                             
  559.  RV.INCO.'CD'  = 'PROG' 0.;                                             
  560.  RV.INCO.'TP' = 'PROG' 0.;
  561.  RV.INCO.'cont'  = 0.;                                                  
  562.  RV.INCO.'UXa_12'= 'PROG' ('NBNO' axe_12)*0. ;                          
  563.  RV.INCO.'UXa_34'= 'PROG' ('NBNO' axe_34)*0. ;
  564.  RV.INCO.'Pobst' = 'PROG' (('NBNO' obstacle) '+' 1 )*0. ;       
  565.  rv.INCO.'Ka_12'= 'PROG' ('NBNO' axe_12)*0. ;
  566.  rv.INCO.'Ka_34'= 'PROG' ('NBNO' axe_34)*0. ;
  567.  
  568.     lh=  (POIN SURFA PROC  (2. 7.)) 'ET'  (POIN SURFA PROC (6. 7.))      
  569.    'ET'  (POIN SURFA PROC (8.  7.)) 'ET'  (POIN SURFA PROC (10. 7.))     
  570.    'ET'  (POIN SURFA PROC (13. 7.)) 'ET'  (POIN SURFA PROC (18. 7.));    
  571.  
  572.    his= 'KHIS' 'UN' 1 lh
  573.                'UN' 2 lh
  574.                'KN'   lh
  575.                'EN'   lh;
  576.    rv.hist=his; 
  577.  
  578.  EXEC  RV;
  579.  
  580. * 'OPTION' donn 5;
  581. * 'FIN' ;
  582.  
  583.  rt=rv.pasdetps; 
  584.  tps1 = rt.tps;
  585.  it1  = rt.nupasdt;
  586.  dt   = rt.'DELTAT-1';
  587.  
  588.  'MESSAGE'  'No ITERATION '  it1;
  589.  'MESSAGE'  'DERNIER PAS '   DT;
  590.  'MESSAGE'  'TEMPS PHYSIQUE' tps1;
  591.  
  592. * 'OPTION'  donn 5;
  593.  
  594. *****HISTORIQUES*****
  595. *
  596. TAB1=his.TABD;                                                             
  597. TAB1.'TITRE'= TABLE ;                                                   
  598. TAB1.'TITRE' . 1  = MOT 'COURBE 1';                                  
  599. TAB1.'TITRE' . 2  = MOT 'COURBE 2';
  600. TAB1.'TITRE' . 3  = MOT 'COURBE 3';
  601. TAB1.'TITRE' . 4  = MOT 'COURBE 4';
  602. TAB1.'TITRE' . 5  = MOT 'COURBE 5';
  603. TAB1.'TITRE' . 6  = MOT 'COURBE 6';
  604.  
  605.  si ('EGA' graph 'O' );
  606.  'DESSIN' (his.'1UN') 'TITRE' 'HISTORIQUES POUR UX ' LEGE TAB1;
  607.  'DESSIN' (his.'2UN') 'TITRE' 'HISTORIQUES POUR UY ' LEGE TAB1;
  608.  'DESSIN' (his.'KN') 'TITRE' 'HISTORIQUES POUR K' LEGE TAB1;
  609.  'DESSIN' (his.'EN') 'TITRE' 'HISTORIQUES POUR E' LEGE TAB1;
  610.  finsi ;
  611.  
  612. ** PROFILS SUR LE DOMAINE TOTAL
  613. **--------------------------------
  614. vv = 'VECTEUR'  rv.inco.'UN' .1 ux uy  ;
  615. Vx = rv.inco.'UN' exco UX SCAL ;
  616. Vy = rv.inco.'UN' exco UY SCAL ;
  617. kk =  rv.inco.'KN' ;
  618. EE =  rv.inco.'EN';
  619. nuts = 'ELNO' ('KCHT' $surfa scal centre (RV.INCO.'NUT')) $surfa;
  620. PP   = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  621.  
  622.  si ('EGA' graph 'O' );
  623. *'OPTION' isov ligne;
  624.  titre 'Vecteur Vitesse';
  625.  'TRACER' vv surfa ('CONTOUR'  surfa); 
  626.  titre 'vecteur vitesse composante X';
  627.  'TRACER'  Vx surfa ('CONTOUR'  surfa) ;
  628.  titre 'vecteur vitesse composante Y';
  629.  'TRACER' VY surfa ('CONTOUR'  surfa) ;
  630.  titre 'CHAMP DE PRESSION';
  631.  'TRACER' PP surfa ('CONTOUR'  surfa); 
  632.  titre 'Energie Cinetique Turbulente';
  633.  'TRACER'  kk surfa ('CONTOUR'  surfa);
  634.  titre 'dissipation';
  635.  'TRACER' ee surfa ('CONTOUR' surfa); 
  636.  titre 'NUT';
  637.  'TRACER' nuts surfa ('CONTOUR' surfa); 
  638.  finsi ;
  639.  
  640. * 'OPTION' donn 5;
  641.  
  642. ** FONCTION DE COURANT 
  643. **---------------------
  644.  
  645.  un1='EXCO' (rv.inco.'UN') ux;
  646.  un2='EXCO' (rv.inco.'UN') uy;
  647.  
  648.  un1='KCHT' $surfa SCAL SOMMET un1;
  649.  un2='KCHT' $surfa SCAL SOMMET un2;
  650.  
  651. ** CALCUL DU ROTATIONNEL DU CHAMP DE VITESSE
  652.  
  653.  rt2d='KOPS' (rv.inco.'UN') 'ROT' $surfa;
  654.  sw= (-1*rt2d);
  655.  
  656. ** DEFINITION DE Utangentielle SUR CHAQUE FACE
  657. ** !! la direction de Ut depend de la définition 
  658.  
  659.  unnh=KCHT $HAUT  SCAL SOMMET (un1);
  660.  unnh=NOEL $HAUT unnh;
  661.  unnh=KCHT $HAUT SCAL CENTRE unnh;
  662.  
  663.  unnb=KCHT $bas  SCAL SOMMET (-1*un1);
  664.  unnb=NOEL $bas unnb;
  665.  unnb=KCHT $bas SCAL CENTRE unnb;
  666.  
  667.  unne=KCHT $entree  SCAL SOMMET (un2);
  668.  unne=NOEL $entree unne;
  669.  unne=KCHT $entree SCAL CENTRE unne;
  670.  
  671.  unns=KCHT $sortie SCAL SOMMET (-1*un2);
  672.  unns=NOEL $sortie unns;
  673.  unns=KCHT $sortie SCAL CENTRE unns;
  674.  
  675.    rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
  676.    ZONE $SURFA  OPER LAPN 1. inco 'psi'
  677.    ZONE $SURFA  OPER FIMP sw inco 'psi'
  678.    ZONE $HAUT   OPER FIMP unnh inco 'psi'
  679.    ZONE $bas    OPER FIMP unnb inco 'psi'
  680.    ZONE $sortie OPER FIMP unns inco 'psi'
  681.    ZONE $entree OPER FIMP unne inco 'psi'
  682.  
  683.    CLIM 'psi' TIMP obstacle 0. ;
  684.  
  685.    rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
  686.    exec rk;
  687.  
  688.  si ('EGA' graph 'O' );
  689.    psi=rk.inco.'psi'; 
  690.  'OPTION' ISOV LIGNE;
  691.   toto = 'PROG' -0.51 PAs 0.02   Npas 50;  
  692.  'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
  693.  finsi ;
  694.  
  695. ** CALCUL DES COEFFS AERODYNAMIQUES
  696. **------------------------------
  697. ** Trainee 'ET' portance
  698.  cd = 2*(rv.inco.'CD');
  699.  cl = 2*(rv.inco.'CL');
  700.  tp = 'ENLEVER' (rv.inco.'TP') 1;
  701.  cd = 'ENLEVER' cd  1; 
  702.  cl = 'ENLEVER' cl  1;
  703.  
  704. Ntp = ('DIME' tp);
  705.  
  706. PDT1 =  'EXTRAIRE' tp 1;
  707. PDTN =  'EXTRAIRE' tp Ntp;
  708.  
  709. ** duree moyenne de  20  pas de tps
  710. PDTx20 = (PDTN '-'PDT1 )/( Ntp '-' 1);
  711.  
  712. ** Pas de tps  moyen
  713. PDTmoy = (PDTx20 '/' 20);
  714.  
  715. *'LISTE' PDTmoy;
  716. *'OPTION' donn 5;
  717.  
  718.  ecd = 'EVOL' 'MANUEL' tp cd;
  719.  ecl = 'EVOL' 'MANUEL' tp cl; 
  720.  moycd = MOYL cd;
  721.  moycl = MOYL cl;
  722.  rmscd = RMS cd;
  723.  rmscl = RMS cl;
  724.  
  725.  'LISTE' moycd; 'LISTE' rmscd; 'LISTE' moycl; 'LISTE' rmscl; 
  726.  
  727.  si ('EGA' graph 'O' );
  728.  'DESSIN' ecl 'TITRE' 'cl';
  729.  'DESSIN' ecd 'TITRE' 'cd';
  730.  
  731.   DSPCL = DSPR 8 ecl Fmax 0.8;
  732.  'DESSIN' DSPCL; 
  733.  finsi ;
  734.  
  735. ** Evolution du champ de pression autour de l'obstacle
  736. ** Presion instantanee
  737.  Pref  = 'EXTRAIRE' PP Scal  P32;
  738.  Pref = Pref '-' 0.5;
  739.  
  740.  Pobst = 'REDU' PP obstacle;
  741.  CP    = 2*(Pobst '-' Pref);
  742.  eCP   = 'EVOL' chpo CP scal ('INVERSE' obstacle) ;
  743.  
  744. ** On enleve les valeurs aux extremites
  745.  Oobst = 'EXTRAIRE' eCP ordo;
  746.  Oobst = 'ENLEVER' Oobst 1 ;
  747.  Aobst = 'EXTRAIRE' eCP absc;
  748.  Aobst  = 'ENLEVER' Aobst 1 ;
  749.  NbAobst = 'DIME' Aobst ;
  750.  Aobst = 'ENLEVER' Aobst NbAobst;
  751.  Oobst = 'ENLEVER' Oobst NbAobst;
  752.  eCP = 'EVOL' 'MANUEL' Aobst Oobst;
  753.  
  754. si ('EGA' graph 'O' );
  755.  'DESSIN' (eCP)  'TITRE' 'pression instantanee  sur obstacle' ;
  756. finsi ;
  757.  
  758. ** Pression moyenne le long de l'obstacle
  759.  PM =(rv.INCO.'Pobst') '/' (rv.INCO.'cont');
  760.  PM = 2*(PM '-'('PROG'('DIME' PM)*Pref ));
  761.  PM = 'ENLEVER' PM 1;
  762.  PM =  'ENLEVER' PM ('DIME' Pm);
  763.  ePM = 'EVOL' 'MANUEL' Aobst PM;
  764. si ('EGA' graph 'O' );
  765.  'DESSIN' (epm) 'TITRE' 'Pression moyenne sur obstacle';
  766. finsi ;
  767.  
  768. ** EVOLUTIONS LE LONG DE L'AXE
  769. **-------------------------------
  770. ** ABSC reelles des points de l'axe
  771.  Xaxe_12 = 'EVOL' chpo (COOR 1 axe_12) scal axe_12; 
  772.  Xaxe_34 = 'EVOL' chpo (COOR 1 axe_34) scal axe_34;
  773.  Xaxe_12 = 'EXTRAIRE' Xaxe_12 ordo;
  774.  Xaxe_34 = 'EXTRAIRE' Xaxe_34 ordo;
  775.  
  776. ** Pour tracer Uy
  777.   XUY_34 = ( 'PROG' 5.49) 'ET' Xaxe_34;
  778.  
  779. ** Evolution instantanne le long de l'axe (fct de l'abscisse curviligne)
  780.  
  781.  eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
  782.  eUYa_12 = EVOL CHPO (rv.inco.'UN') UY axe_12;
  783.  eKa_12  = EVOL CHPO (rv.inco.'KN') scal axe_12;
  784.  
  785.  eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
  786.  eUYa_34 = EVOL CHPO (rv.inco.'UN') uy axe_34;
  787.  eKa_34  = EVOL CHPO (rv.inco.'KN') scal axe_34;
  788.  
  789. ** valeurs le long de l'axe  
  790.  UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
  791.  UYa_12 = 'EXTRAIRE' eUYa_12 ORDO;
  792.  Ka_12  = 'EXTRAIRE' eKa_12 ordo;
  793.  
  794.  UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
  795.  UYa_34 = 'EXTRAIRE' eUYa_34 ordo;
  796.  UYa_34 = ('PROG' 0.) 'ET' UYa_34;
  797.  Ka_34  = 'EXTRAIRE' eKa_34  ordo;
  798.  
  799. ** Evolution le long de l'axe (fct de l'abscisse relle)
  800.  
  801.  eUXa_12 = EVOL VERT 'MANUEL'  Xaxe_12 UXa_12;
  802.  eUXa_34 = EVOL VERT 'MANUEL'  Xaxe_34 UXa_34;
  803.  eUXaxe = eUXa_12 'ET' eUXa_34;
  804.  
  805.  eUya_12 = EVOL VERT 'MANUEL'  Xaxe_12 Uya_12;
  806.  eUya_34 = EVOL VERT 'MANUEL'  XUY_34 Uya_34;
  807.  eUyaxe = eUya_12 'ET' eUya_34;
  808.  
  809.  eKa_12 = EVOL VERT 'MANUEL'  Xaxe_12 Ka_12;
  810.  eKa_34 = EVOL VERT 'MANUEL'  Xaxe_34 Ka_34;
  811.  eKaxe = eKa_12 'ET' eKa_34;
  812.  
  813. ** Ux moyen le long de l'axe
  814.  
  815.  UM_12 =  (rv.INCO.'UXa_12') '/' (rv.INCO.'cont');
  816.  UM_34 =  (rv.INCO.'UXa_34') '/' (rv.INCO.'cont');
  817.  eUM_12 = EVOL VERT 'MANUEL'  Xaxe_12 UM_12;
  818.  eUM_34 = EVOL VERT 'MANUEL'  Xaxe_34 UM_34;
  819.  eUMaxe = eUM_12 'ET' eUM_34;
  820.  
  821. ** K moyen le long de l'axe
  822.  KM_12 =  (rv.INCO.'Ka_12') '/' (rv.INCO.'cont');
  823.  KM_34 =  (rv.INCO.'Ka_34') '/' (rv.INCO.'cont');
  824.  eKM_12 = EVOL VERT 'MANUEL'  Xaxe_12 KM_12;
  825.  eKM_34 = EVOL VERT 'MANUEL'  Xaxe_34 KM_34;
  826.  eKMaxe = eKM_12 'ET' eKM_34;
  827.  
  828. ** Axes  
  829.  LX0_12 = 'PROG' 0. 4.5; LX0_34 = 'PROG' 5.5 20.;
  830.  LY0 = 'PROG' 0. 0.;
  831.  L0_12 = 'EVOL'  'MANUEL' LX0_12 Ly0;
  832.  L0_34 = 'EVOL'  'MANUEL' LX0_34 Ly0;
  833.  L0 = L0_12 'ET' L0_34; 
  834.  
  835.  si ('EGA' graph 'O' );
  836. 'DESSIN' (eUMaxe 'ET' L0 ) 'TITRE' 'Ux moyen sur axe';
  837.  'DESSIN'  (eUXaxe  'ET' L0) titr 'PROFIL Ux sur axe'
  838.   tity 'Ux/U0' titx 'X/D';
  839.  'DESSIN' (eUya_34  'ET' L0_34) titr 'PROFIL Uy sur axe'
  840.   tity 'Uy/U0' titx 'X/D';
  841.  'DESSIN' (eKMaxe 'ET' L0 ) 'TITRE' 'K  moyen sur axe';
  842.  'DESSIN' (eKaxe 'ET' L0) titr 'PROFIL K sur axe'
  843.   tity 'K/U0**2' titx 'X/D';
  844.  finsi ;
  845.  
  846. Fin;
  847.  
  848.  
  849.  
  850.  
  851. ** Pas de tps  moyen 
  852. *PDTmoy = (PDTx20 '/' 20);
  853.  
  854.  
  855. *'LISTE' PDTmoy;
  856. *'OPTION' donn 5;
  857.  
  858. ** Evolution reguliere dans le tps
  859. * RegTPS = 'PROG' PDT1 PAS PDTx20 Npas (Nval '-' 1) ;
  860.  
  861. *'OPTION'  donn 5;
  862. *
  863. UYFDT_1= 'EXTRAIRE'  his.'2UN' ORDO 1;
  864. UYFDT_2= 'EXTRAIRE'  his.'2UN' ORDO 2;
  865. UYFDT_3= 'EXTRAIRE'  his.'2UN' ORDO 3;
  866. UYFDT_4= 'EXTRAIRE'  his.'2UN' ORDO 4;
  867. UYFDT_5= 'EXTRAIRE'  his.'2UN' ORDO 5;
  868. UYFDT_6= 'EXTRAIRE'  his.'2UN' ORDO 6;
  869. *
  870. eUYFDT_1 = 'EVOL' 'MANUEL' hisTPS UYFDT_1;
  871. eUYFDT_2 = 'EVOL' 'MANUEL' hisTPS UYFDT_2;
  872. eUYFDT_3 = 'EVOL' 'MANUEL' hisTPS UYFDT_3;
  873. eUYFDT_4 = 'EVOL' 'MANUEL' hisTPS UYFDT_4;
  874. eUYFDT_5 = 'EVOL' 'MANUEL' hisTPS UYFDT_5;
  875. eUYFDT_6 = 'EVOL' 'MANUEL' hisTPS UYFDT_6;
  876.  
  877. UxFDT_1= 'EXTRAIRE'  his.'1UN' ORDO 1;
  878. UxFDT_2= 'EXTRAIRE'  his.'1UN' ORDO 2;
  879. UXFDT_3= 'EXTRAIRE'  his.'1UN' ORDO 3;
  880. UxFDT_4= 'EXTRAIRE'  his.'1UN' ORDO 4;
  881. UxFDT_5= 'EXTRAIRE'  his.'1UN' ORDO 5;
  882. UxFDT_6= 'EXTRAIRE'  his.'1UN' ORDO 6;
  883. *
  884. eUxFDT_1 = 'EVOL' 'MANUEL' hisTPS UxFDT_1;
  885. eUxFDT_2 = 'EVOL' 'MANUEL' hisTPS UxFDT_2;
  886. eUxFDT_3 = 'EVOL' 'MANUEL' hisTPS UxFDT_3;
  887. eUxFDT_4 = 'EVOL' 'MANUEL' hisTPS UxFDT_4;
  888. eUxFDT_5 = 'EVOL' 'MANUEL' hisTPS UxFDT_5;
  889. eUxFDT_6 = 'EVOL' 'MANUEL' hisTPS UxFDT_6;
  890. *
  891. *TABUY=Table;
  892. * TABUY.1 ='MARQ ETOI';
  893.  
  894.  'DESSIN' eUYFDT_1 'TITRE''EVOLUTION UY_axe en X=1.94'
  895.   tity 'UY' titx 'tps';
  896.  'DESSIN' eUYFDT_2 'TITRE''EVOLUTION UY_axe en X=5.9'
  897.   tity 'UY' titx 'tps' ;
  898.  'DESSIN' eUYFDT_3  'TITRE''EVOLUTION UY_axe en X=8'
  899.   tity 'UY' titx 'tps'
  900.   'DESSIN' eUYFDT_4 'TITRE''EVOLUTION UY_axeen X=10'
  901.   tity 'UY' titx 'tps' ;
  902.  'DESSIN' eUYFDT_5 'TITRE''EVOLUTION UY_axe en fonction du temps
  903.   en X=13' tity 'UY' titx 'tps' ;
  904.  'DESSIN' eUYFDT_6 'TITRE''EVOLUTION UY_axe en fonction du temps
  905.   en X=18'  tity 'UY' titx 'tps' ;
  906.  
  907.  'DESSIN' eUxFDT_1 'TITRE''EVOLUTION Ux_axe en fonction du temps
  908.   en X=1.94'   tity'Ux' titx 'tps' ;
  909.  'DESSIN' eUxFDT_2 'TITRE''EVOLUTION Ux_axe en fonction dutemps
  910.   en X=5.9'  tity 'Ux' titx 'tps' ;;
  911.  'DESSIN' eUxFDT_3 'TITRE''EVOLUTION Ux_axe en fonction du temps
  912.   en X=8' tity 'Ux' titx 'tps' ;;
  913.  'DESSIN' eUxFDT_4 'TITRE''EVOLUTION Ux_axe en fonction du temps
  914.   en X=10' tity 'TABUY Ux' titx 'tps' ;;
  915.  'DESSIN' eUxFDT_5 'TITRE''EVOLUTION Ux_axe en fonction du temps
  916.   en X=13' tity 'Ux' titx 'tps' ;;
  917.  'DESSIN' eUxFDT_6 'TITRE''EVOLUTION Ux_axe en fonction du temps
  918.   en X=18'  tity 'Ux' titx 'tps' ;;
  919.  
  920.  
  921. 'OPTION'  donn 5;
  922.  
  923. *********PROFILS SUR LE DOMAINE TOTAL********
  924.  
  925. vv = 'VECTEUR'  rv.inco.'UN' .1 ux uy  ;
  926.  
  927. Vx = rv.inco.'UN' exco UX SCAL ;
  928.  
  929. Vy = rv.inco.'UN' exco UY SCAL ;
  930.  
  931. kk =  rv.inco.'KN' ;
  932.  
  933. kk_ent = EXCO(REDU kk entree) SCAL ;
  934.  
  935. EE =  rv.inco.'EN';
  936.  
  937. ee_ent = 'EXCO' ('REDU' ee entree) SCAL ;
  938.  
  939. nuts = 'ELNO' ('KCHT' $surfa scal centre (RV.INCO.'NUT')) $surfa;
  940.  
  941. PP   = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  942.  
  943. *'OPTION' isov ligne;
  944.  
  945. * titre 'Vecteur Vitesse';
  946. * 'TRACER' vv surfa ('CONTOUR'  surfa); 
  947.  
  948. *  titre 'vecteur vitesse composante X';
  949. * 'TRACER'  Vx surfa ('CONTOUR'  surfa) ;
  950.  
  951. * titre 'vecteur vitesse composante Y';
  952. * 'TRACER' VY surfa ('CONTOUR'  surfa) ;
  953.  
  954. * titre 'CHAMP DE PRESSION';
  955. * 'TRACER' PP surfa ('CONTOUR'  surfa); 
  956.  
  957. * titre 'Energie Cinetique Turbulente';
  958. * 'TRACER'  kk surfa ('CONTOUR'  surfa);
  959.  
  960. * titre 'dissipation';
  961. * 'TRACER' ee surfa ('CONTOUR' surfa); 
  962.  
  963. * titre 'NUT';
  964. * 'TRACER' nuts surfa ('CONTOUR' surfa); 
  965.  
  966. * 'OPTION' donn 5;
  967.  
  968. ** FONCTION DE COURANT
  969. ** PREMIERE METHODE
  970. **-------------------
  971. ** CALCUL DU ROTATIONNEL DU CHAMP DE VITESSE
  972.  
  973. *  sw ='KOPS' (rv.inco.'UN') 'ROT' $surfa;
  974.  
  975.  
  976. *  DOMLIM = Haut 'ET' BAS 'ET' entree ;
  977.  
  978. *  'ELIMINATION' domliM 1.e-5 ;
  979.  
  980. *  CHPY = 'COORDONNEE' 2 DOMLIM ;
  981.  
  982. *  PSin = CHPY ;
  983.  
  984. *  CHPTOT = 'KCHT' $surfa scal sommet psin ;
  985.  
  986. *   rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
  987. *   ZONE $SURFA  OPER LAPN 1. inco 'psi'
  988. *   ZONE $SURFA  OPER FIMP sw inco 'psi' 
  989.  
  990. *   CLIM 'psi' TIMP DOMLIM  CHPTOT ;
  991.  
  992. *   rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
  993. *   exec rk;
  994.  
  995. *   psi=rk.inco.'psi';
  996. * 
  997. *  'OPTION' ISOV LIGNE;
  998.  
  999. *   'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
  1000.  
  1001. *  'OPTION' donn 5 ; 
  1002.  
  1003.  
  1004. ** FONCTION DE COURANT 
  1005. ** SECONDE METHODE
  1006. **---------------------
  1007.  
  1008.  un1='EXCO' (rv.inco.'UN') ux;
  1009.  un2='EXCO' (rv.inco.'UN') uy;
  1010.  
  1011.  un1='KCHT' $surfa SCAL SOMMET un1;
  1012.  un2='KCHT' $surfa SCAL SOMMET un2;
  1013.  
  1014. ** CALCUL DU ROTATIONNEL DU CHAMP DE VITESSE
  1015.  
  1016.  rt2d='KOPS' (rv.inco.'UN') 'ROT' $surfa;
  1017.  sw= (-1*rt2d);
  1018.  
  1019. ** DEFINITION DE Utangentielle SUR CHAQUE FACE
  1020. ** !! la direction de Ut depend de la définition 
  1021. ** de la normale a la paroi n  !!
  1022.  
  1023.  
  1024.  unnh=KCHT $HAUT  SCAL SOMMET (un1);
  1025.  unnh=NOEL $HAUT unnh;
  1026.  unnh=KCHT $HAUT SCAL CENTRE unnh;
  1027.  
  1028.  unnb=KCHT $bas  SCAL SOMMET (-1*un1);
  1029.  unnb=NOEL $bas unnb;
  1030.  unnb=KCHT $bas SCAL CENTRE unnb;
  1031.  
  1032.  unne=KCHT $entree  SCAL SOMMET (un2);
  1033.  unne=NOEL $entree unne;
  1034.  unne=KCHT $entree SCAL CENTRE unne;
  1035.  
  1036.  unns=KCHT $sortie SCAL SOMMET (-1*un2);
  1037.  unns=NOEL $sortie unns;
  1038.  unns=KCHT $sortie SCAL CENTRE unns;
  1039.  
  1040.  
  1041.  
  1042.    rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
  1043.    ZONE $SURFA  OPER LAPN 1. inco 'psi'
  1044.    ZONE $SURFA  OPER FIMP sw inco 'psi'
  1045.    ZONE $HAUT   OPER FIMP unnh inco 'psi'
  1046.    ZONE $bas    OPER FIMP unnb inco 'psi'
  1047.    ZONE $sortie OPER FIMP unns inco 'psi'
  1048.    ZONE $entree OPER FIMP unne inco 'psi'
  1049.  
  1050.    CLIM 'psi' TIMP obstacle 0. ;
  1051.  
  1052.    rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
  1053.    exec rk;
  1054.  
  1055.    psi=rk.inco.'psi';
  1056. * 
  1057.   'OPTION' ISOV LIGNE;
  1058.     toto = 'PROG' -0.51 PAs 0.02   Npas 50;  
  1059.  
  1060.    'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
  1061.  
  1062.  
  1063.  
  1064.   'OPTION' donn 5 ; 
  1065. * 'FIN' ;
  1066. *
  1067.  
  1068.  
  1069. ** CALCUL DES COEFFS AERODYNAMIQUES
  1070. **------------------------------
  1071. ** Traine 'ET' portance
  1072.  cd = 2*(rv.inco.'CD');
  1073.  cl = 2*(rv.inco.'CL');
  1074.  tp = 'ENLEVER' (rv.inco.'TP') 1;
  1075.  cd = 'ENLEVER' cd  1; 
  1076.  cl = 'ENLEVER' cl  1;
  1077.  
  1078.  cd = tronc cd 400 ;
  1079.  cl = tronc cl 400 ;
  1080.  tp = tronc tp 400 ;
  1081.  
  1082. Ntp = ('DIME' tp);
  1083.  
  1084. PDT1 =  'EXTRAIRE' tp 1;
  1085. PDTN =  'EXTRAIRE' tp Ntp;
  1086.  
  1087. ** duree moyenne de  20  pas de tps
  1088. PDTx20 = (PDTN '-'PDT1 )/( Ntp '-' 1);
  1089.  
  1090. ** Pas de tps  moyen
  1091. PDTmoy = (PDTx20 '/' 20);
  1092.  
  1093. 'LISTE' PDTmoy;
  1094. 'OPTION' donn 5;
  1095.  
  1096. ** Evolution reguliere dans le tps
  1097.  RegTP = 'PROG' PDT1 PAS PDTx20 Npas (Ntp '-' 1) ;
  1098. *'OPTION'  donn 5;
  1099.  
  1100.  ecd = 'EVOL' 'MANUEL' tp cd;
  1101.  ecl = 'EVOL' 'MANUEL' tp cl; 
  1102.  
  1103.  moycd = MOYL cd;
  1104.  moycl = MOYL cl;
  1105.  
  1106.  rmscd = RMS cd;
  1107.  rmscl = RMS cl;
  1108.  
  1109.  'LISTE' moycd;
  1110.  'LISTE' rmscd;
  1111.  'LISTE' moycl;
  1112.  'LISTE' rmscl; 
  1113.  
  1114.  'DESSIN' ecl 'TITRE' 'cl';
  1115.  'DESSIN' ecd 'TITRE' 'cd';
  1116.  
  1117.  DSPCL = DSPR 8 ecl Fmax 0.8;
  1118.  'DESSIN' DSPCL; 
  1119.  
  1120.  'OPTION' donn 5;
  1121.  
  1122. * Evolution du champ de pression autour de l'obstacle
  1123. ** Presion instantanee
  1124.  Pref  = 'EXTRAIRE' PP Scal  P32;
  1125.  Pref = Pref '-' 0.5;
  1126.  
  1127.  Pobst = 'REDU' PP obstacle;
  1128.  CP    = 2*(Pobst '-' Pref);
  1129.  eCP   = 'EVOL' chpo CP scal ('INVERSE' obstacle) ;
  1130.  
  1131.  
  1132. * titi = 'PROG' Pref;
  1133.  
  1134.  'OPTION' isov ligne;
  1135.  titre 'CHAMP DE PRESSION';
  1136. * 'TRACER' PP surfa toto ('CONTOUR'  surfa);
  1137. * 'TRACER' PP surfa ('CONTOUR'  surfa); 
  1138.  
  1139.  
  1140. * 'OPTION' donn 5;
  1141.  
  1142. ** On enleve les valeurs aux extremites
  1143.  
  1144.  Oobst = 'EXTRAIRE' eCP ordo;
  1145.  Oobst = 'ENLEVER' Oobst 1 ;
  1146.  
  1147.  Aobst = 'EXTRAIRE' eCP absc;
  1148.  Aobst  = 'ENLEVER' Aobst 1 ;
  1149.  
  1150.  NbAobst = 'DIME' Aobst ;
  1151.  
  1152.  Aobst = 'ENLEVER' Aobst NbAobst;
  1153.  Oobst = 'ENLEVER' Oobst NbAobst;
  1154.  
  1155.  eCP = 'EVOL' 'MANUEL' Aobst Oobst;
  1156.  
  1157.  'DESSIN' (eCP)  'TITRE' 'pression instantanee  sur obstacle' ;
  1158.  
  1159. ** Pression moyenne le long de l'obstacle
  1160.  PM =(rv.INCO.'Pobst') '/' (rv.INCO.'cont');
  1161.  PM = 2*(PM '-'('PROG'('DIME' PM)*Pref ));
  1162.  PM = 'ENLEVER' PM 1;
  1163.  PM =  'ENLEVER' PM ('DIME' Pm);
  1164.  
  1165.  ePM = 'EVOL' 'MANUEL' Aobst PM;
  1166.  'DESSIN' (epm) 'TITRE' 'Pression moyenne sur obstacle';
  1167.  
  1168. *  'OPTION' donn 5;
  1169.  
  1170. ** 
  1171. ** EVOLUTIONS LE LONG DE L'AXE
  1172. **-------------------------------
  1173. ** ABSC reelles des points de l'axe
  1174.  Xaxe_12 = 'EVOL' chpo (COOR 1 axe_12) scal axe_12; 
  1175.  Xaxe_34 = 'EVOL' chpo (COOR 1 axe_34) scal axe_34;
  1176.  
  1177.  Xaxe_12 = 'EXTRAIRE' Xaxe_12 ordo;
  1178. * 'LISTE'  Xaxe_12;
  1179.  Xaxe_34 = 'EXTRAIRE' Xaxe_34 ordo;
  1180.  
  1181. ** Pour tracer Uy
  1182.   XUY_34 = ( 'PROG' 5.49) 'ET' Xaxe_34;
  1183.  
  1184. * 'LISTE'  Xaxe_34;
  1185.  
  1186. * 'OPTION' donn 5;
  1187.  
  1188. ** Evolution instantanne
  1189. ** le long de l'axe (fct de l'abscisse curviligne)
  1190.  
  1191.  eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
  1192.  eUYa_12 = EVOL CHPO (rv.inco.'UN') UY axe_12;
  1193.  eKa_12  = EVOL CHPO (rv.inco.'KN') scal axe_12;
  1194.  
  1195.  eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
  1196.  eUYa_34 = EVOL CHPO (rv.inco.'UN') uy axe_34;
  1197.  eKa_34  = EVOL CHPO (rv.inco.'KN') scal axe_34;
  1198.  
  1199. **valeurs le long de l'axe  
  1200.  UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
  1201.  UYa_12 = 'EXTRAIRE' eUYa_12 ORDO;
  1202.  Ka_12  = 'EXTRAIRE' eKa_12 ordo;
  1203.  
  1204.  UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
  1205.  
  1206.  UYa_34 = 'EXTRAIRE' eUYa_34 ordo;
  1207.  UYa_34 = ('PROG' 0.) 'ET' UYa_34;
  1208.  
  1209.  Ka_34  = 'EXTRAIRE' eKa_34  ordo;
  1210.  
  1211. ** Evolution le long de l'axe (fct de l'abscisse relle)
  1212.  
  1213.  eUXa_12 = EVOL VERT 'MANUEL'  Xaxe_12 UXa_12;
  1214.  eUXa_34 = EVOL VERT 'MANUEL'  Xaxe_34 UXa_34;
  1215.  
  1216.  eUXaxe = eUXa_12 'ET' eUXa_34;
  1217.  
  1218.  eUya_12 = EVOL VERT 'MANUEL'  Xaxe_12 Uya_12;
  1219.  eUya_34 = EVOL VERT 'MANUEL'  XUY_34 Uya_34;
  1220.  
  1221.  eUyaxe = eUya_12 'ET' eUya_34;
  1222.  
  1223.  eKa_12 = EVOL VERT 'MANUEL'  Xaxe_12 Ka_12;
  1224.  eKa_34 = EVOL VERT 'MANUEL'  Xaxe_34 Ka_34;
  1225.  
  1226.  eKaxe = eKa_12 'ET' eKa_34;
  1227.  
  1228. ** Ux moyen le long de l'axe
  1229.  
  1230.  UM_12 =  (rv.INCO.'UXa_12') '/' (rv.INCO.'cont');
  1231.  UM_34 =  (rv.INCO.'UXa_34') '/' (rv.INCO.'cont');
  1232.  
  1233.  eUM_12 = EVOL VERT 'MANUEL'  Xaxe_12 UM_12;
  1234.  eUM_34 = EVOL VERT 'MANUEL'  Xaxe_34 UM_34;
  1235.  
  1236.  eUMaxe = eUM_12 'ET' eUM_34;
  1237.  
  1238. ** K moyen le long de l'axe
  1239.  
  1240.  KM_12 =  (rv.INCO.'Ka_12') '/' (rv.INCO.'cont');
  1241.  KM_34 =  (rv.INCO.'Ka_34') '/' (rv.INCO.'cont');
  1242.  
  1243.  eKM_12 = EVOL VERT 'MANUEL'  Xaxe_12 KM_12;
  1244.  eKM_34 = EVOL VERT 'MANUEL'  Xaxe_34 KM_34;
  1245.  
  1246.  eKMaxe = eKM_12 'ET' eKM_34;
  1247.  
  1248. * opti donn 5;
  1249.  
  1250.  TAB2=TABLE;                                                             
  1251. * TAB2.'TITRE'= TABLE ;                                                   
  1252. * TAB2.1='MARQ ETOI NOLI';           
  1253. * TAB2.'TITRE' . 1  = MOT 'Calcul';                                  
  1254. * TAB2.2='MARQ ETOI NOLI';
  1255. * TAB2.'TITRE' . 2  = MOT 'en entree';
  1256.  TAB2.3='TIRL';                                                     
  1257. * TAB2.'TITRE' . 3  = MOT 'en entree';
  1258.  TAB2.4='TIRL';
  1259.  
  1260. ** affichage du cube
  1261.  LX0_12 = 'PROG' 0. 4.5;
  1262.  LX0_34 = 'PROG' 5.5 20.;
  1263.  LY0 = 'PROG' 0. 0.;
  1264.  L0_12 = 'EVOL'  'MANUEL' LX0_12 Ly0;
  1265.  L0_34 = 'EVOL'  'MANUEL' LX0_34 Ly0;
  1266.  L0 = L0_12 'ET' L0_34; 
  1267.  
  1268.  CUBEX  = 'PROG'  4.5 4.5 5.5 5.5 5.5  4.5 4.5;
  1269.  CUBEX2 = 'PROG'  4.49 4.49 5.5 5.5 5.5  4.49 4.49;
  1270.  CUBE1  = 'PROG'   0. 0.05 0.05 0. -0.05  -0.05  0. ;
  1271.  CUBE2  = 'PROG'   0. 0.02 0.02 0. -0.02  -0.02  0. ;
  1272.  CUBE3  = 'PROG'   0. 1e-3 1e-3 0. -1e-3  -1e-3  0. ;
  1273.  
  1274.  CUBEux = 'EVOL' ROUGE 'MANUEL' CUBEx CUBE1;
  1275.  CUBEuy = 'EVOL' ROUGE 'MANUEL' CUBEx2 CUBE2;
  1276.  CUBEK  = 'EVOL' ROUGE 'MANUEL' CUBEx CUBE3;
  1277.  
  1278.  
  1279. ** SORTIES GRAPHIQUES
  1280. **-------------------
  1281. *** VITESSE
  1282. **-------------------------------
  1283.  
  1284.  'DESSIN' (eUMaxe 'ET' L0 'ET' cubeux ) 'TITRE' 'Ux moyen sur axe';
  1285.  
  1286.  'DESSIN'  (eUXaxe  'ET' L0 'ET' cubeux) titr 'PROFIL Ux sur axe'
  1287.  TAB2 tity 'Ux/U0' titx 'X/D';
  1288.  
  1289. * dess (eUyaxe  'ET' L0 'ET' cubeuy) titr 'PROFIL Uy sur axe'
  1290. * TAB2 tity 'Uy/U0' titx 'X/D';
  1291.  
  1292.  dess (eUya_34  'ET' L0_34 'ET' cubeuy) titr 'PROFIL Uy sur axe'
  1293.  tity 'Uy/U0' titx 'X/D';
  1294.  
  1295.  
  1296.  'DESSIN' (eKMaxe 'ET' L0 'ET' cubeK ) 'TITRE' 'K  moyen sur axe';
  1297.  
  1298.  dess (eKaxe 'ET' L0 'ET' cubeK ) titr 'PROFIL K sur axe'
  1299.  TAB2 tity 'K/U0**2' titx 'X/D';
  1300.  
  1301. * opti donn 5;
  1302.  
  1303. fin;
  1304.  
  1305.  
  1306.  
  1307.  
  1308.  
  1309.  
  1310.  
  1311.  
  1312.  
  1313.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales