Télécharger carre_expl.dgibi

Retour à la liste

Numérotation des lignes :

  1. *****************************************************
  2. ************************************************************************
  3. ************************************************************************
  4. * fichier : carre_expl.dgibi *
  5. ** modifie le 15/06/2014 passage EQPR -> EQEX *
  6. *****************************************************
  7. ** ECOULEMENT AUTOUR D'UNE CYLINDRE DE SECTION CAREE*
  8. ** Gregory Turbelin 29/12/1998
  9. ** modifie le 4/02/2000
  10. **
  11. ** CAS TEST DE LYN ET RODI
  12. ** RESULTATS EXPERIMENTAUX REFERENCES A L_ADRESSE
  13. ** http://vortex.mech.surrey.ac.uk/database/test43/test43.html
  14. *****************************************************
  15.  
  16. GRAPH='N' ;
  17. COMPLET=FAUX ;
  18. 'OPTION' 'DIME' 2 'ELEM' qua8;
  19.  
  20. DISCR = 'MACRO' ;
  21. KPRESS = 'CENTRE';
  22. BETAP = 1. ;
  23.  
  24. **************************
  25. ** PROCEDURES DE CALCUL
  26. **************************
  27. **--------------------------
  28. ***************************************
  29. **Calcul des valeurs moyennes sur l'axe
  30. ***************************************
  31. ** En chaque point de l'axe, et a chaque iteration,
  32. ** on additionne les valeurs de la variable. La moyenne, en un point donne,
  33. ** sera la somme des valeurs, divisee par le nbre de valeurs stockees
  34.  
  35. 'DEBPROC' VMPM;
  36. 'ARGUMENT' rvx*'TABLE';
  37. rv= rvx.'EQEX';
  38.  
  39. DD = rv.PASDETPS.'NUPASDT' ;
  40. NN = DD/20;
  41. T0 = (DD '-' (20*NN)) 'EGA' 0;
  42. 'SI' T0;
  43.  
  44. PP = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  45. UN = RV.INCO.'UN' ;
  46.  
  47. unx= kcht $surfa scal sommet (exco 'UX' un) ;
  48.  
  49. ** compteur pour connaitre le nombre de valeurs stockées
  50. cont = 1.;
  51. rv.INCO.'cont'= (rv.INCO.'cont' '+' cont);
  52.  
  53. ** VITESSE LE LONG DE L'AXE
  54. eUXa_12 = EVOL CHPO UNX axe_12;
  55. UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
  56.  
  57. eUXa_34 = EVOL CHPO UNX axe_34;
  58. UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
  59.  
  60. rv.INCO.'UXa_12'= (rv.INCO.'UXa_12' '+' UXa_12 );
  61. rv.INCO.'UXa_34'= (rv.INCO.'UXa_34' '+' UXa_34 );
  62.  
  63. ** PRESSION SUR L'OBSTACLE
  64.  
  65. ePobst = 'EVOL' chpo PP scal ('INVERSE' obstacle);
  66. Pobst = 'EXTRAIRE' ePobst ordo;
  67. rv.INCO.'Pobst'= (rv.INCO.'Pobst' '+' Pobst);
  68.  
  69. ** ENERGIE LE LONG DE L'AXE
  70. eka_12 = EVOL CHPO (rv.inco.'KN') scal axe_12;
  71. ka_12 = 'EXTRAIRE' eKa_12 ORDO;
  72.  
  73. eKa_34 = EVOL CHPO (rv.inco.'KN') scal axe_34;
  74. ka_34 = 'EXTRAIRE' eka_34 ordo;
  75.  
  76. rv.INCO.'Ka_12'= (rv.INCO.'Ka_12' '+' Ka_12 );
  77. rv.INCO.'Ka_34'= (rv.INCO.'Ka_34' '+' Ka_34 );
  78.  
  79. 'FINSI' ;
  80. as2 ama1 = 'KOPS' 'MATRIK' ;
  81. 'FINPROC' as2 ama1 ;
  82.  
  83. ***************************************
  84. **CALCUL DES COEFF AERODYNAMIQUES
  85. **PAR INTEGRATION DU CHAMP DE PRESSION
  86. **SUR LES PAROIS DE L_OBSTACLE
  87. ***************************************
  88. **------------------------------------------------
  89. 'DEBPROC' CLCD;
  90. 'ARGUMENT' rvx*'TABLE';
  91. rv = rvx.'EQEX';
  92.  
  93. DD = rv.PASDETPS.'NUPASDT' ;
  94. NN = DD/20;
  95. T0 = (DD '-' (20*NN)) 'EGA' 0;
  96. 'SI' T0;
  97.  
  98. PP = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  99.  
  100. ** PRESSION FACE1
  101. PP1 = 'REDU' PP face1 ;
  102. ** PRESSION FACE2
  103. PP2 = 'REDU' PP face2;
  104. **PRESSION FACE3
  105. PP3 = 'REDU' PP face3;
  106. **PRESSION FACE4
  107. PP4 = 'REDU' PP face4;
  108.  
  109. **COORD. face1
  110. XX1 = REDU CHPX face1 ;
  111. **COORD. face2
  112. YY2 = REDU CHPY face2 ;
  113. **COORD. face3
  114. XX3 = REDU CHPX face3 ;
  115. **COORD. face4
  116. YY4 = REDU CHPY face4 ;
  117.  
  118. ** CHPO ---> EVOL (CHPO)
  119. ePP1 = EVOL CHPO PP1 face1 ;
  120. ePP2 = EVOL CHPO PP2 face2 ;
  121. ePP3 = EVOL CHPO PP3 face3 ;
  122. ePP4 = EVOL CHPO PP4 face4 ;
  123.  
  124. eXX1 = EVOL CHPO XX1 face1 ;
  125. eXX3 = EVOL CHPO XX3 face3 ;
  126. eYY2 = EVOL CHPO YY2 face2 ;
  127. eYY4 = EVOL CHPO YY4 face4 ;
  128.  
  129. ** EVOL (CHPO) --> LISTREEL
  130.  
  131. P1 = EXTR ePP1 ORDO;
  132. P2 = EXTR ePP2 ORDO;
  133. P3 = EXTR ePP3 ORDO;
  134. P4 = EXTR ePP4 ORDO;
  135.  
  136. X1 = EXTR eXX1 ORDO;
  137. X3 = EXTR eXX3 ORDO;
  138. Y2 = EXTR eYY2 ORDO;
  139. Y4 = EXTR eYY4 ORDO;
  140.  
  141. ** CALCUL DE LA PORTANCE ET TRAINEE
  142. *
  143. ** LISTREEL --> EVOL (MANU)
  144. evXP1 = EVOL MANU X1 P1;
  145. evXP3 = EVOL MANU X3 P3;
  146.  
  147. evyP2 = EVOL MANU Y2 P2 ;
  148. evYP4 = EVOL MANU Y4 P4 ;
  149.  
  150. ** CALCUL INTEGRALE (Pdx)
  151. intXP1 = INTG evXP1 ;
  152. intXP3 = INTG evXP3 ;
  153.  
  154. ** CALCUL INTEGRALE (Pdy)
  155. intYP2 = INTG evYP2 ;
  156. intYP4 = INTG evYP4 ;
  157.  
  158. CL = ((-1)*intXP3)+ intXP1 ;
  159. CD = (intYP4)+((-1)*intYP2) ;
  160.  
  161. rv.inco.'CL' = (rv.inco.'CL') et (CL);
  162. rv.inco.'CD' = (rv.inco.'CD') et (CD);
  163.  
  164. PT = 'PROG' rv.'PASDETPS'.'TPS' ;
  165. rv.INCO.'TP'= rv.INCO.'TP' ET PT;
  166.  
  167. 'FINSI' ;
  168. as2 ama1 = 'KOPS' 'MATRIK' ;
  169.  
  170. 'FINPROC' as2 ama1 ;
  171.  
  172. ****************************************
  173. ** PROCEDURE POUR ESTIMER LA CONVERGENCE
  174. ****************************************
  175. ** -----------------------------------------------
  176. 'DEBPROC' ERR;
  177. 'ARGUMENT' rvx*'TABLE';
  178. rv= rvx.'EQEX';
  179.  
  180. DD = rv.PASDETPS.'NUPASDT' ;
  181. NN = DD/20;
  182. T0 = (DD '-' (20*NN)) 'EGA' 0;
  183. 'SI' T0;
  184.  
  185. UN = RV.INCO.'UN' ;
  186. UNM1 = RV.INCO.'UNM1' ;
  187.  
  188. unx= kcht $surfa scal sommet (exco 'UX' un) ;
  189. unm1x = kcht $surfa scal sommet (exco 'UX' unm1) ;
  190. uny= kcht $surfa scal sommet (exco 'UY' un) ;
  191. unm1y = kcht $surfa scal sommet (exco 'UY' unm1) ;
  192.  
  193. ERRX = KOPS unx '-' unm1x ;
  194. ERRY = KOPS uny '-' unm1y ;
  195. ELIX = MAXI ERRX 'ABS' ;
  196. ELIY = MAXI ERRY 'ABS' ;
  197. ELIX = (LOG (ELIX + 1.0E-20))/(LOG 20.) ;
  198. ELIY = (LOG (ELIY + 1.0E-20))/(LOG 20.) ;
  199.  
  200. **! pour afficher NUT/NU, attention aux dimensions
  201.  
  202. MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' ' ERREUR LINF ' ELIX ELIY
  203. 'MAX NUT/NU = ' ((MAXI RV.INCO.'NUT')*Re) ;
  204. IT = PROG RV.PASDETPS.'NUPASDT' ;
  205.  
  206. ERY = PROG ELIY ;
  207. ERX = PROG ELIX ;
  208.  
  209. RV.INCO.'IT' = (RV.INCO.'IT') ET IT ;
  210. RV.INCO.'ERY' = (RV.INCO.'ERY') ET ERY ;
  211. RV.INCO.'ERX' = (RV.INCO.'ERX') ET ERX ;
  212.  
  213. 'FINSI' ;
  214. RV.INCO.'UNM1'= KCHT $surfa vect sommet (RV.INCO.'UN') ;
  215.  
  216. as2 ama1 = 'KOPS' 'MATRIK' ;
  217. 'FINPROC' as2 ama1 ;
  218.  
  219. **-------------------------------------------
  220. ** Procedures pour estimer la valeur moyenne
  221. ** et l_ecart type d'une liste de relles
  222. **-------------------------------------------
  223. 'DEBPROC' MOYL LR*'LISTREEL';
  224. LSi =0.;
  225. i=1;
  226. 'REPETER' boucle ('DIME' LR);
  227. ielem = 'EXTRAIRE' LR i;
  228. LSi = LSi '+' ielem;
  229. moyenne = LSi/i;
  230. i =i+1;
  231. 'FIN' boucle;
  232. 'FINPROC' moyenne ;
  233. ** Ecart 'TYPE'
  234. 'DEBPROC' RMS LR*'LISTREEL';
  235. LMi = 0.;
  236. LEi = 0.;
  237. i=1;
  238. 'REPETER' boucle ('DIME' LR);
  239. ielem = 'EXTRAIRE' LR i;
  240. ielem2 = (ielem**2);
  241.  
  242. LMi = LMi '+' ielem;
  243. LEi = LEi '+' ielem2;
  244.  
  245. moyenne = LMi/i;
  246. moyenne2 = LEi/i;
  247. i =i+1;
  248. 'FIN' boucle;
  249. ecart = (moyenne2 '-' ((moyenne)**2))**(0.5);
  250. 'FINPROC' ecart ;
  251. ***********************************
  252. ** CONSTRUCTION DU MAILLAGE (ADIM)
  253. *************************************
  254. **--------------------------------
  255. *** Diametre du cylindre (longueur de reference)
  256. *** D_ref = 0.04 m <==> X2 (= Y2+Y3)
  257.  
  258. **Longueurs sans dim
  259.  
  260. X1 = 4.5 ; X2 = 1. ; X3 = 14.5;
  261. Y1 = 6.5 ; Y2 = 0.5 ; Y3 = 0.5 ; Y4=6.5 ;
  262. *
  263. **Points caracteristiques
  264. *
  265. P11=0. 0.;
  266. P12=X1 0.;
  267. P13=(X1+X2) 0.;
  268. P14=(X1+X2+X3) 0.;
  269. *
  270. P21=0. Y1;
  271. P22=X1 Y1;
  272. P23=(X1+X2) Y1;
  273. P24=(X1+X2+X3) Y1;
  274. *
  275. P31=0. (Y1+Y2);
  276. P32=X1 (Y1+Y2);
  277. P33=(X1+X2) (Y1+Y2);
  278. P34=(X1+X2+X3) (Y1+Y2);
  279. *
  280. P41=0. (Y1+Y2+Y3);
  281. P42=X1 (Y1+Y2+Y3);
  282. P43=(X1+X2) (Y1+Y2+Y3);
  283. P44=(X1+X2+X3) (Y1+Y2+Y3);
  284.  
  285. P51=0. (Y1+Y2+Y3+Y4);
  286. P52=X1 (Y1+Y2+Y3+Y4);
  287. P53=(X1+X2) (Y1+Y2+Y3+Y4);
  288. P54=(X1+X2+X3) (Y1+Y2+Y3+Y4);
  289. *
  290. **Nb de mailles horizontales, sur obstacle
  291. si complet ;
  292. Nx = 6;
  293. sinon ;
  294. Nx = 2;
  295. finsi ;
  296. **Nb de mailles verticales 1/2 obstacle
  297. si complet ;
  298. Ny = 3;
  299. sinon ;
  300. Ny = 1;
  301. finsi ;
  302. **Nb de mailles horizontales, avant obstacle
  303. Nx1= 3*Nx;
  304. **Nb de mailles horizontales, sur obstacle
  305. Nx2 = Nx ;
  306. **Nb de mailles horizontales, après obstacle
  307. Nx3 =7*Nx;
  308. *
  309. **Nb de mailles verticales sous obstacle
  310. Ny1 = 8*Ny;
  311. **Nb de mailles verticales 1/2 obstacle inf
  312. Ny2 = Ny;
  313. **Nb de mailles verticales 1/2 obstacle sup
  314. Ny3 = Ny;
  315. **Nb de mailles verticales au dessus obstacle
  316. Ny4 = 8*Ny;
  317. *
  318. *
  319. ** Definition des densites
  320. **-------------------------
  321. DX1 = (X1/Nx1);
  322. Dx1_1 = Dx1 ; Dx1_2 = 0.3*Dx1;
  323. *
  324. DX2 = (X2/Nx2);
  325. Dx2_2 = Dx2 ; Dx2_3 = Dx2;
  326. *
  327. DX3 = (X3/Nx3);
  328. Dx3_3 = 0.3*Dx3 ; Dx3_4 = 1.5*Dx3;
  329. *
  330. Dy1 = (Y1/NY1);
  331. Dy1_1 =1.*DY1 ; Dy1_2 = 0.3*DY1;
  332. *
  333. DY2 = (Y2/NY2);
  334. DY2_2 = DY2; DY2_3 = DY2;
  335. *
  336. DY3 = (Y3/NY3);
  337. DY3_3 =0.5*DY3 ; DY3_4 = DY3;
  338. *
  339. DY4 = (Y4/NY4);
  340. DY4_4 =0.3*DY4 ; DY4_5 = 1.*DY4;
  341. *
  342.  
  343. ********************************************
  344. **LIGNES POUR UN PARAMETRAGE ZONES PAR ZONES
  345. **--------------------------------
  346. L1x_12 = P11 droite (-1*Nx1) P12 'DINI' Dx1_1 'DFIN' Dx1_2;
  347. L1x_23 = P12 droite (-1*Nx2) P13 'DINI' Dx2_2 'DFIN' Dx2_3;
  348. L1x_34 = P13 droite (-1*Nx3) P14 'DINI' Dx3_3 'DFIN' Dx3_4;
  349.  
  350. L2x_12 = P21 droite (-1*Nx1) P22 'DINI' Dx1_1 'DFIN' Dx1_2;
  351. L2x_23 = P22 droite (-1*Nx2) P23 'DINI' Dx2_2 'DFIN' Dx2_3;
  352. L2x_34 = P23 droite (-1*Nx3) P24 'DINI' Dx3_3 'DFIN' Dx3_4;
  353.  
  354. L3x_12 = P31 droite (-1*Nx1) P32 'DINI' Dx1_1 'DFIN' Dx1_2;
  355. L3x_23 = P32 droite (-1*Nx2) P33 'DINI' Dx2_2 'DFIN' Dx2_3;
  356. L3x_34 = P33 droite (-1*Nx3) P34 'DINI' Dx3_3 'DFIN' Dx3_4;
  357.  
  358. L4x_12 = P41 droite (-1*Nx1) P42 'DINI' Dx1_1 'DFIN' Dx1_2;
  359. L4x_23 = P42 droite (-1*Nx2) P43 'DINI' Dx2_2 'DFIN' Dx2_3;
  360. L4x_34 = P43 droite (-1*Nx3) P44 'DINI' Dx3_3 'DFIN' Dx3_4;
  361.  
  362. L5x_12 = P51 droite (-1*Nx1) P52 'DINI' Dx1_1 'DFIN' Dx1_2;
  363. L5x_23 = P52 droite (-1*Nx2) P53 'DINI' Dx2_2 'DFIN' Dx2_3;
  364. L5x_34 = P53 droite (-1*Nx3) P54 'DINI' Dx3_3 'DFIN' Dx3_4;
  365.  
  366. L1x = L1x_12 'ET' L1x_23 'ET' L1x_34;
  367. L2x = L2x_12 'ET' L2x_23 'ET' L2x_34;
  368. L3x = L3x_12 'ET' L3x_23 'ET' L3x_34;
  369. L4x = L4x_12 'ET' L4x_23 'ET' L4x_34;
  370. L5x = L5x_12 'ET' L5x_23 'ET' L5x_34;
  371.  
  372. L1y_12 = P11 droite (-1*Ny1) P21 'DINI' Dy1_1 'DFIN' Dy1_2;
  373. L1y_23 = P21 droite (-1*Ny2) P31 'DINI' Dy2_2 'DFIN' Dy2_3;
  374. L1y_34 = P31 droite (-1*Ny3) P41 'DINI' Dy3_3 'DFIN' Dy3_4;
  375. L1y_45 = P41 droite (-1*Ny4) P51 'DINI' Dy4_4 'DFIN' Dy4_5;
  376.  
  377. L2y_12 = P12 droite (-1*Ny1) P22 'DINI' Dy1_1 'DFIN' Dy1_2;
  378. L2y_23 = P22 droite (-1*Ny2) P32 'DINI' Dy2_2 'DFIN' Dy2_3;
  379. L2y_34 = P32 droite (-1*Ny3) P42 'DINI' Dy3_3 'DFIN' Dy3_4;
  380. L2y_45 = P42 droite (-1*Ny4) P52 'DINI' Dy4_4 'DFIN' Dy4_5;
  381.  
  382. L3y_12 = P13 droite (-1*Ny1) P23 'DINI' Dy1_1 'DFIN' Dy1_2;
  383. L3y_23 = P23 droite (-1*Ny2) P33 'DINI' Dy2_2 'DFIN' Dy2_3;
  384. L3y_34 = P33 droite (-1*Ny3) P43 'DINI' Dy3_3 'DFIN' Dy3_4;
  385. L3y_45 = P43 droite (-1*Ny4) P53 'DINI' Dy4_4 'DFIN' Dy4_5;
  386.  
  387. L4y_12 = P14 droite (-1*Ny1) P24 'DINI' Dy1_1 'DFIN' Dy1_2;
  388. L4y_23 = P24 droite (-1*Ny2) P34 'DINI' Dy2_2 'DFIN' Dy2_3;
  389. L4y_34 = P34 droite (-1*Ny3) P44 'DINI' Dy3_3 'DFIN' Dy3_4;
  390. L4y_45 = P44 droite (-1*Ny4) P54 'DINI' Dy4_4 'DFIN' Dy4_5;
  391.  
  392. L1y = L1y_12 'ET' L1y_23 'ET' L1y_34 'ET' L1y_45;
  393. L2y = L2y_12 'ET' L2y_23 'ET' L2y_34 'ET' L2y_45;
  394. L3y = L3y_12 'ET' L3y_23 'ET' L3y_34 'ET' L3y_45;
  395. L4y = L4y_12 'ET' L4y_23 'ET' L4y_34 'ET' L4y_45;
  396.  
  397. bas = (L1x) 'COULEUR' vert ; 'ELIMINATION' bas 1.e-5;
  398. sortie = (L4y) 'COULEUR' jaune; 'ELIMINATION' sortie 1.e-5;
  399. haut = (L5x) 'COULEUR' bleu ; 'ELIMINATION' haut 1.e-5;
  400. entree = (L1y) 'COULEUR' rouge; 'ELIMINATION' entree 1.e-5;
  401. axe_12 = L3x_12; axe_34=L3x_34;
  402. axe = (axe_12 'ET' axe_34) 'COULEUR' turq; 'ELIMINATION' axe 1.e-5;
  403.  
  404. face1 = L2x_23 'COULEUR' vert ;
  405. face2 = (L3y_23 'ET' L3y_34) 'COULEUR' jaune;
  406. face3 = L4x_23 'COULEUR' bleu ;
  407. face4 = (L2y_23 'ET' L2y_34) 'COULEUR' rouge ;
  408.  
  409. obstacle = FACE1 'ET' FACE2 'ET' FACE3 'ET' FACE4 ;
  410. 'ELIMINATION' obstacle 1.e-5; ;
  411.  
  412. **Creation de sous-zones
  413. ZX1Y1 = 'DALLER' L1x_12 L2y_12 L2x_12 L1y_12 PLAN ;
  414. ZX1Y2 = 'DALLER' L2x_12 L2y_23 L3x_12 L1y_23 PLAN ;
  415. ZX1Y3 = 'DALLER' L3x_12 L2y_34 L4x_12 L1y_34 PLAN ;
  416. ZX1Y4 = 'DALLER' L4x_12 L2y_45 L5x_12 L1y_45 PLAN ;
  417.  
  418. ZX1 = ZX1Y1 'ET' ZX1Y2 'ET' ZX1Y3 'ET' ZX1Y4;
  419. 'ELIMINATION' 1.e-5 ZX1 ;
  420.  
  421. ZX2Y1 = 'DALLER' L1x_23 L3y_12 L2x_23 L2y_12 PLAN ;
  422. ZX2Y4 = 'DALLER' L4x_23 L3y_45 L5x_23 L2y_45 PLAN ;
  423.  
  424. ZX2 = ZX2Y1 'ET' ZX2Y4;
  425. 'ELIMINATION' ZX2 1.e-5;
  426. orienter Zx2;
  427.  
  428. ZX3Y1 = 'DALLER' L1x_34 L4y_12 L2x_34 L3y_12 PLAN ;
  429. ZX3Y2 = 'DALLER' L2x_34 L4y_23 L3x_34 L3y_23 PLAN ;
  430. ZX3Y3 = 'DALLER' L3x_34 L4y_34 L4x_34 L3y_34 PLAN ;
  431. ZX3Y4 = 'DALLER' L4x_34 L4y_45 L5x_34 L3y_45 PLAN ;
  432.  
  433. ZX3 = ZX3Y1 'ET' ZX3Y2 'ET' ZX3Y3 'ET' ZX3Y4;
  434. 'ELIMINATION' ZX3 1.e-5;
  435. orienter Zx3;
  436.  
  437. SURFA = Zx1 'ET' ZX2 'ET' Zx3 ;
  438. 'ELIMINATION' surfa 1.e-5;
  439. orienter surfa;
  440.  
  441. surfa = 'CHANGER' surfa 'QUAF' ;
  442. bas = 'CHANGER' bas 'QUAF' ;
  443. entree = 'CHANGER' entree 'QUAF' ;
  444. haut = 'CHANGER' haut 'QUAF' ;
  445. sortie = 'CHANGER' sortie 'QUAF' ;
  446. face1 = 'CHANGER' face1 'QUAF' ;
  447. face2 = 'CHANGER' face2 'QUAF' ;
  448. face3 = 'CHANGER' face3 'QUAF' ;
  449. face4 = 'CHANGER' face4 'QUAF' ;
  450. axe_12 = 'CHANGER' axe_12 'QUAF' ;
  451. axe_34 = 'CHANGER' axe_34 'QUAF' ;
  452.  
  453. 'ELIMINATION' 1.D-6 (bas 'ET' entree 'ET' haut 'ET' sortie
  454. et face1 et face2 et face3 et face4 et axe_12 et axe_34
  455. 'ET' surfa) ;
  456.  
  457. *************************
  458. ** DEFINITION DES MODELES
  459. *************************
  460. **------------------------
  461. $surfa = 'MODE' surfa 'NAVIER_STOKES' DISCR ;
  462. $bas = 'MODE' bas 'NAVIER_STOKES' DISCR ;
  463. $entree = 'MODE' entree 'NAVIER_STOKES' DISCR ;
  464. $haut = 'MODE' haut 'NAVIER_STOKES' DISCR ;
  465. $sortie = 'MODE' sortie 'NAVIER_STOKES' DISCR ;
  466. $face1 = 'MODE' face1 'NAVIER_STOKES' DISCR ;
  467. $face2 = 'MODE' face2 'NAVIER_STOKES' DISCR ;
  468. $face3 = 'MODE' face3 'NAVIER_STOKES' DISCR ;
  469. $face4 = 'MODE' face4 'NAVIER_STOKES' DISCR ;
  470. $axe_12 = 'MODE' axe_12 'NAVIER_STOKES' DISCR ;
  471. $axe_34 = 'MODE' axe_34 'NAVIER_STOKES' DISCR ;
  472.  
  473. surfa = 'DOMA' $surfa 'MAILLAGE' ;
  474. 'ORIENTER' surfa;
  475.  
  476. entree = 'DOMA' $entree 'MAILLAGE' ;
  477. haut = 'DOMA' $haut 'MAILLAGE' ;
  478. bas = 'DOMA' $bas 'MAILLAGE' ;
  479. sortie = 'DOMA' $sortie 'MAILLAGE' ;
  480. face1 = 'DOMA' $face1 'MAILLAGE' ;
  481. face2 = 'DOMA' $face2 'MAILLAGE' ;
  482. face3 = 'DOMA' $face3 'MAILLAGE' ;
  483. face4 = 'DOMA' $face4 'MAILLAGE' ;
  484. axe_34 = 'DOMA' $axe_34 'MAILLAGE' ;
  485. axe_12 = 'DOMA' $axe_12 'MAILLAGE' ;
  486.  
  487. si ('EGA' graph 'O' );
  488. 'TRACER' surfa;
  489. finsi ;
  490.  
  491. ** ABSC ET ORDO des points du maillage
  492. CHPX CHPY = COOR SURFA;
  493.  
  494. *****************************
  495. *** PROPRIETES PHYSIQUES
  496. ***************************
  497.  
  498. ** Hauteur du domaine (sans dim)
  499. Hda = (Y1 + Y2 + Y3 + Y4);
  500.  
  501. ** Longueur du domaine (sans dim)
  502. Lda = (X1+X2+X3) ;
  503.  
  504. ** Reynolds de l'écoulement
  505. Re = 2.2e4;
  506. iRe = 1/Re ;
  507.  
  508. ** Distances Normales aux parois de l'obstacle
  509. ** exprimees en fonction de la taille des mailles
  510. ** sur la face_1
  511. Xface1 = coor 1 face1;
  512. eXface1 = evol chpo Xface1 face1;
  513. eXface1 = 'EXTRAIRE' eXface1 ABSC;
  514.  
  515. ** Yp = 1/5 hauteur de la premiere maille (sans dim)
  516. YP = ('EXTR' eXface1 2)/5;
  517.  
  518. ** vitesse de frottement sans dim
  519. ufini = 150./(yp*Re);
  520.  
  521. **constante de Von Karman
  522. C1 = 0.4;
  523.  
  524. ** constante du modele RNG k-e
  525. CNU = 0.0845;
  526.  
  527. ** taux de turbulence
  528. INT = 0.02 ;
  529.  
  530. alpha = 1.5 ;
  531.  
  532. ** Energie cinetique
  533. ** Kadim = K/(U_ref**2)
  534.  
  535. Kk = alpha*((INT)**2) ;
  536.  
  537. ** Taille des tourbillons energetiques
  538. ** On suppose que k**0.5/Lk = U0/L0
  539. ** ou L0 est une dimension caracteristique
  540. ** associe a l_ecoulement moyen
  541.  
  542. L0 = 1. ;
  543.  
  544. Lk = (alpha**0.5) * INT * L0;
  545.  
  546. ** Dissipation
  547.  
  548. Ek = (Kk**1.5)'/'(Lk);
  549.  
  550. **Reynolds associe aux tourbillons energetiques
  551.  
  552. Rek = alpha *((INT)**2) * (Re*L0) ;
  553.  
  554. **NUT/NU
  555.  
  556. NUTSNU = cnu * ((Kk**2)/Ek) * Re ;
  557.  
  558. ** On doit, normalement, avoir NUTSNU = cnu *Rek
  559.  
  560. **Verification
  561. Liste Lda ;
  562. List Hda ;
  563. 'LISTE' Re;
  564. 'LISTE' nu;
  565. 'LISTE' yp;
  566. 'LISTE' ufini;
  567. list lk ;
  568. list kk ;
  569. list Ek ;
  570. list Rek ;
  571. list (CNU*Rek) ;
  572. list NUTSNU;
  573.  
  574. ** Definition des profils d'entrees et initiaux
  575. **----------------------------------------------
  576. ** Definition du domaine d'imposition
  577. M_imp =entree;
  578.  
  579. ** VITESSES
  580. **-----------
  581. U1 = 1. 0.;
  582. U1 = 'KCHT' $surfa vect sommet U1;
  583.  
  584. ** Valeurs IMPOSEES
  585. UX_IMP = 'EXCO' ('REDU' U1 M_imp) 'UX' SCAL;
  586. UY_IMP = 'EXCO' ('REDU' U1 M_imp) 'UY' SCAL;
  587.  
  588. ** Valeurs INITIALES (sans dim)
  589. U_ini = 'MANU' 'CHPO' surfa 1 scal 0.;
  590. U_ini = U_ini 'NOMC' 'UX' ;
  591. U_ini = 'KCHT' $surfa vect sommet U_ini;
  592.  
  593. **ENERGIE CINETIQUE
  594. **-----------------
  595. K1 = 'MANU' 'CHPO' surfa 1 scal Kk;
  596.  
  597. ***INITIALE
  598. K_ini = 'KCHT' $surfa scal sommet K1;
  599. ***IMPOSEE
  600. K_IMP = 'KCHT' $surfa scal sommet K1;
  601. K_IMP= 'EXCO' ('REDU' K_IMP M_imp) SCAL;
  602.  
  603. **DISSIPATION
  604. **--------------
  605.  
  606. E1 = 'MANU' 'CHPO' surfa 1 scal Ek;
  607.  
  608. ** INITIALE
  609. E_ini = 'KCHT' $surfa scal sommet E1;
  610. **IMPOSEE
  611. E_IMP = 'KCHT' $surfa scal sommet E1;
  612. E_IMP='EXCO' ('REDU' E_IMP M_imp) SCAL;
  613.  
  614. **NUT INITIALE
  615. **------------------
  616. *fonction de K et E
  617. N1=CNU*(K_ini**2)*(E_ini**-1);
  618. N1 = 'KCHT' $surfa scal sommet N1;
  619. N_ini = 'NOEL' $surfa N1 ;
  620.  
  621. ** U* INITIALES
  622. **-------------
  623. UET_ini1 = 'KCHT' $face1 scal centre Ufini;
  624. UET_ini2 = 'KCHT' $face2 scal centre Ufini;
  625. UET_ini3 = 'KCHT' $face3 scal centre Ufini;
  626. UET_ini4 = 'KCHT' $face4 scal centre Ufini;
  627.  
  628. ** Module de resolution
  629. **-------------------
  630.  
  631. RV = 'EQEX' $surfa ITMA 42 ALFA 0.8
  632. OPTI 'RNG'
  633. ZONE $surfa OPER FILTREKE 1. Hda iRe 'UN' INCO 'KN' 'EN'
  634. ZONE $surfa OPER NSKE iRe 'NUT' INCO 'UN' 'KN' 'EN'
  635. ZONE $face4 OPER FPU iRe 'UET4' YP INCO 'UN' 'KN' 'EN'
  636. ZONE $face2 OPER FPU iRe 'UET2' YP INCO 'UN' 'KN' 'EN'
  637. ZONE $face3 OPER FPU iRe 'UET3' YP INCO 'UN' 'KN' 'EN'
  638. ZONE $face1 OPER FPU iRe 'UET1' YP INCO 'UN' 'KN' 'EN' ;
  639.  
  640. RV = EQEX RV 'OPTI' 'CENTREE'
  641. 'ZONE' $surfa 'OPER' 'DFDT' 1. 'UN' 'DELTAT' INCO 'UN'
  642. 'ZONE' $surfa 'OPER' 'DFDT' 1. 'KN' 'DELTAT' INCO 'KN'
  643. 'ZONE' $surfa 'OPER' 'DFDT' 1. 'EN' 'DELTAT' INCO 'EN' ;
  644.  
  645. RV = 'EQEX' RV
  646. ZONE $surfa OPER ERR
  647. ZONE $surfa OPER VMPM
  648. ZONE $surfa OPER CLCD ;
  649.  
  650. RV = 'EQEX' RV
  651. CLIM 'UN' UIMP M_imp UX_imp
  652. CLIM 'UN' VIMP M_imp UY_imp
  653. CLIM 'UN' VIMP (haut 'ET' bas) 0.
  654. CLIM 'KN' TIMP M_imp K_imp
  655. CLIM 'EN' TIMP M_imp E_imp ;
  656.  
  657. RVP = EQEX 'OPTI' 'EF' KPRESS
  658. 'ZONE' $surfa OPER KBBT -1. betap 'INCO' 'UN' 'PRES'
  659. 'ZONE' $face4 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  660. 'ZONE' $face2 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  661. 'ZONE' $face3 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  662. 'ZONE' $face1 OPER VNIMP $surfa 0. 'INCO' 'UN' 'PRES'
  663. ;
  664.  
  665. RV.'PROJ' = RVP ;
  666.  
  667.  
  668. RV.INCO = TABLE INCO;
  669. RV.INCO.'PRES'= kcht $surfa SCAL KPRESS 0.;
  670. RV.INCO.'UN' = U_ini;
  671. RV.INCO.'KN' = K_ini;
  672. RV.INCO.'EN' = E_ini;
  673. RV.INCO.'NUT' = N_ini;
  674. RV.INCO.'UET4' = UET_ini4;
  675. RV.INCO.'UET3' = UET_ini3;
  676. RV.INCO.'UET2' = UET_ini2;
  677. RV.INCO.'UET1' = UET_ini1;
  678. RV.INCO.'CL' = 'PROG' ;
  679. RV.INCO.'CD' = 'PROG' ;
  680. RV.INCO.'TP' = 'PROG' ;
  681. RV.INCO.'cont' = 0.;
  682. RV.INCO.'UXa_12'= 'PROG' ('NBNO' axe_12)*0. ;
  683. RV.INCO.'UXa_34'= 'PROG' ('NBNO' axe_34)*0. ;
  684. RV.INCO.'Pobst' = 'PROG' (('NBNO' obstacle) '+' 1 )*0. ;
  685. rv.INCO.'Ka_12'= 'PROG' ('NBNO' axe_12)*0. ;
  686. rv.INCO.'Ka_34'= 'PROG' ('NBNO' axe_34)*0. ;
  687. rv.inco.'IT' = PROG 1 ;
  688. RV.INCO.'UNM1'= KCHT $surfa vect sommet (1.E-3 1.E-3) ;
  689. rv.inco.'ERY' = PROG 0. ;
  690. rv.inco.'ERX' = PROG 0. ;
  691.  
  692. EXEC RV;
  693.  
  694. * 'OPTION' donn 5;
  695. * 'FIN' ;
  696.  
  697. rt=rv.pasdetps;
  698. tps1 = rt.tps;
  699. it1 = rt.nupasdt;
  700. dt = rt.'DELTAT-1';
  701.  
  702. 'MESSAGE' 'No ITERATION ' it1;
  703. 'MESSAGE' 'DERNIER PAS ' DT;
  704. 'MESSAGE' 'TEMPS PHYSIQUE' tps1;
  705.  
  706. ** PROFILS SUR LE DOMAINE TOTAL
  707. **--------------------------------
  708. vv = 'VECTEUR' rv.inco.'UN' .1 ux uy ;
  709. Vx = rv.inco.'UN' exco UX SCAL ;
  710. Vy = rv.inco.'UN' exco UY SCAL ;
  711. kk = rv.inco.'KN' ;
  712. EE = rv.inco.'EN';
  713. nuts = 'ELNO' ('KCHT' $surfa scal centre (RV.INCO.'NUT')) $surfa;
  714. NUTSNU = nuts*Re ;
  715. PP = 'ELNO' $surfa (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  716.  
  717. si ('EGA' graph 'O' );
  718. *'OPTION' isov ligne;
  719. titre 'Vecteur Vitesse à t=' tps1;
  720. 'TRACER' vv surfa ('CONTOUR' surfa);
  721. titre 'vecteur vitesse composante X à t=' tps1;
  722. 'TRACER' Vx surfa ('CONTOUR' surfa) ;
  723. titre 'vecteur vitesse composante Y à t=' tps1;
  724. 'TRACER' VY surfa ('CONTOUR' surfa) ;
  725. titre 'CHAMP DE PRESSION';
  726. 'TRACER' PP surfa ('CONTOUR' surfa);
  727. titre 'Energie Cinetique Turbulente à t=' tps1;
  728. 'TRACER' kk surfa ('CONTOUR' surfa);
  729. titre 'dissipation à t=' tps1;
  730. 'TRACER' ee surfa ('CONTOUR' surfa);
  731. titre 'NUT/NU à t=' tps1;
  732. 'TRACER' nutsnu surfa ('CONTOUR' surfa);
  733. finsi ;
  734.  
  735. ** FONCTION DE COURANT
  736. **---------------------
  737.  
  738. un1='EXCO' (rv.inco.'UN') ux;
  739. un2='EXCO' (rv.inco.'UN') uy;
  740.  
  741. un1='KCHT' $surfa SCAL SOMMET un1;
  742. un2='KCHT' $surfa SCAL SOMMET un2;
  743.  
  744. ** CALCUL DU ROTATIONNEL DU CHAMP DE VITESSE
  745.  
  746. rt2d='KOPS' (rv.inco.'UN') 'ROT' $surfa;
  747. sw= (-1*rt2d);
  748.  
  749. ** DEFINITION DE Utangentielle SUR CHAQUE FACE
  750. ** ! la direction de Ut depend de la définition
  751.  
  752. unnh=KCHT $HAUT SCAL SOMMET (un1);
  753. unnh=NOEL $HAUT unnh;
  754. unnh=KCHT $HAUT SCAL CENTRE unnh;
  755.  
  756. unnb=KCHT $bas SCAL SOMMET (-1*un1);
  757. unnb=NOEL $bas unnb;
  758. unnb=KCHT $bas SCAL CENTRE unnb;
  759.  
  760. unne=KCHT $entree SCAL SOMMET (un2);
  761. unne=NOEL $entree unne;
  762. unne=KCHT $entree SCAL CENTRE unne;
  763.  
  764. unns=KCHT $sortie SCAL SOMMET (-1*un2);
  765. unns=NOEL $sortie unns;
  766. unns=KCHT $sortie SCAL CENTRE unns;
  767.  
  768. rk=EQEX $SURFA 'OPTI' 'EF' 'IMPL'
  769. ZONE $SURFA OPER LAPN 1. inco 'psi'
  770. ZONE $SURFA OPER FIMP sw inco 'psi'
  771. ZONE $HAUT OPER FIMP unnh inco 'psi'
  772. ZONE $bas OPER FIMP unnb inco 'psi'
  773. ZONE $sortie OPER FIMP unns inco 'psi'
  774. ZONE $entree OPER FIMP unne inco 'psi'
  775.  
  776. CLIM 'psi' TIMP obstacle 0. ;
  777.  
  778. rk.inco.'psi'=KCHT $SURFA SCAL SOMMET 0.;
  779. exec rk;
  780.  
  781. si ('EGA' graph 'O' );
  782. psi=rk.inco.'psi';
  783. 'OPTION' ISOV LIGNE;
  784. toto = 'PROG' -0.51 PAs 0.02 Npas 50;
  785. 'TRACER' psi surfa toto ('CONTOUR' (surfa)) titr 'LIGNES DE COURANT';
  786. finsi ;
  787.  
  788. ** CALCUL DES COEFFS AERODYNAMIQUES
  789. **------------------------------
  790. ** Trainee 'ET' portance
  791. cd = 2*(rv.inco.'CD');
  792. cl = 2*(rv.inco.'CL');
  793. tp = RV.INCO.'TP' ;
  794.  
  795. Ntp = ('DIME' tp);
  796.  
  797. ecd = 'EVOL' 'MANUEL' tp cd;
  798. ecl = 'EVOL' 'MANUEL' tp cl;
  799. moycd = MOYL cd;
  800. moycl = MOYL cl;
  801. rmscd = RMS cd;
  802. rmscl = RMS cl;
  803.  
  804. 'LISTE' moycd; 'LISTE' rmscd; 'LISTE' moycl; 'LISTE' rmscl;
  805.  
  806. si ('EGA' graph 'O' );
  807. 'DESSIN' ecl 'TITRE' 'cl';
  808. 'DESSIN' ecd 'TITRE' 'cd';
  809.  
  810. DSPCL = DSPR 8 ecl Fmax 0.8;
  811. 'DESSIN' DSPCL TITR 'DSP';
  812. finsi ;
  813.  
  814. ** Evolution du champ de pression autour de l'obstacle
  815. ** Presion instantanee
  816. Pref = 'EXTRAIRE' PP Scal P32;
  817. Pref = Pref '-' 0.5;
  818.  
  819. Pobst = 'REDU' PP obstacle;
  820. CP = 2*(Pobst '-' Pref);
  821. eCP = 'EVOL' chpo CP scal ('INVERSE' obstacle) ;
  822.  
  823. ** On enleve les valeurs aux extremites
  824. Oobst = 'EXTRAIRE' eCP ordo;
  825. Oobst = 'ENLEVER' Oobst 1 ;
  826. Aobst = 'EXTRAIRE' eCP absc;
  827. Aobst = 'ENLEVER' Aobst 1 ;
  828. NbAobst = 'DIME' Aobst ;
  829. Aobst = 'ENLEVER' Aobst NbAobst;
  830. Oobst = 'ENLEVER' Oobst NbAobst;
  831. eCP = 'EVOL' 'MANUEL' Aobst Oobst;
  832.  
  833. si ('EGA' graph 'O' );
  834. 'DESSIN' (eCP) 'TITRE' 'pression instantanee sur obstacle' ;
  835. finsi ;
  836.  
  837. ** Pression moyenne le long de l'obstacle
  838. PM =(rv.INCO.'Pobst') '/' (rv.INCO.'cont');
  839. PM = 2*(PM '-'('PROG'('DIME' PM)*Pref ));
  840. PM = 'ENLEVER' PM 1;
  841. PM = 'ENLEVER' PM ('DIME' Pm);
  842. ePM = 'EVOL' 'MANUEL' Aobst PM;
  843. si ('EGA' graph 'O' );
  844. 'DESSIN' (epm) 'TITRE' 'Pression moyenne sur obstacle';
  845. finsi ;
  846.  
  847. ** EVOLUTIONS LE LONG DE L'AXE
  848. **-------------------------------
  849. ** ABSC reelles des points de l'axe
  850. Xaxe_12 = 'EVOL' chpo (COOR 1 axe_12) scal axe_12;
  851. Xaxe_34 = 'EVOL' chpo (COOR 1 axe_34) scal axe_34;
  852. Xaxe_12 = 'EXTRAIRE' Xaxe_12 ordo;
  853. Xaxe_34 = 'EXTRAIRE' Xaxe_34 ordo;
  854.  
  855. ** Pour tracer Uy
  856. XUY_34 = ( 'PROG' 5.49) 'ET' Xaxe_34;
  857.  
  858. ** Evolution instantanne le long de l'axe (fct de l'abscisse curviligne)
  859.  
  860. eUXa_12 = EVOL CHPO (rv.inco.'UN') UX axe_12;
  861. eUYa_12 = EVOL CHPO (rv.inco.'UN') UY axe_12;
  862. eKa_12 = EVOL CHPO (rv.inco.'KN') scal axe_12;
  863.  
  864. eUXa_34 = EVOL CHPO (rv.inco.'UN') ux axe_34;
  865. eUYa_34 = EVOL CHPO (rv.inco.'UN') uy axe_34;
  866. eKa_34 = EVOL CHPO (rv.inco.'KN') scal axe_34;
  867.  
  868. ** valeurs le long de l'axe
  869. UXa_12 = 'EXTRAIRE' eUXa_12 ORDO;
  870. UYa_12 = 'EXTRAIRE' eUYa_12 ORDO;
  871. Ka_12 = 'EXTRAIRE' eKa_12 ordo;
  872.  
  873. UXa_34 = 'EXTRAIRE' eUXa_34 ordo;
  874. UYa_34 = 'EXTRAIRE' eUYa_34 ordo;
  875. UYa_34 = ('PROG' 0.) 'ET' UYa_34;
  876. Ka_34 = 'EXTRAIRE' eKa_34 ordo;
  877.  
  878. ** Evolution le long de l'axe (fct de l'abscisse relle)
  879.  
  880. eUXa_12 = EVOL VERT 'MANUEL' Xaxe_12 UXa_12;
  881. eUXa_34 = EVOL VERT 'MANUEL' Xaxe_34 UXa_34;
  882. eUXaxe = eUXa_12 'ET' eUXa_34;
  883.  
  884. eUya_12 = EVOL VERT 'MANUEL' Xaxe_12 Uya_12;
  885. eUya_34 = EVOL VERT 'MANUEL' XUY_34 Uya_34;
  886. eUyaxe = eUya_12 'ET' eUya_34;
  887.  
  888. eKa_12 = EVOL VERT 'MANUEL' Xaxe_12 Ka_12;
  889. eKa_34 = EVOL VERT 'MANUEL' Xaxe_34 Ka_34;
  890. eKaxe = eKa_12 'ET' eKa_34;
  891.  
  892. ** Ux moyen le long de l'axe
  893.  
  894. UM_12 = (rv.INCO.'UXa_12') '/' (rv.INCO.'cont');
  895. UM_34 = (rv.INCO.'UXa_34') '/' (rv.INCO.'cont');
  896. eUM_12 = EVOL VERT 'MANUEL' Xaxe_12 UM_12;
  897. eUM_34 = EVOL VERT 'MANUEL' Xaxe_34 UM_34;
  898. eUMaxe = eUM_12 'ET' eUM_34;
  899.  
  900. ** K moyen le long de l'axe
  901. KM_12 = (rv.INCO.'Ka_12') '/' (rv.INCO.'cont');
  902. KM_34 = (rv.INCO.'Ka_34') '/' (rv.INCO.'cont');
  903. eKM_12 = EVOL VERT 'MANUEL' Xaxe_12 KM_12;
  904. eKM_34 = EVOL VERT 'MANUEL' Xaxe_34 KM_34;
  905. eKMaxe = eKM_12 'ET' eKM_34;
  906.  
  907. ** Axes
  908. LX0_12 = 'PROG' 0. 4.5; LX0_34 = 'PROG' 5.5 20.;
  909. LY0 = 'PROG' 0. 0.;
  910. L0_12 = 'EVOL' 'MANUEL' LX0_12 Ly0;
  911. L0_34 = 'EVOL' 'MANUEL' LX0_34 Ly0;
  912. L0 = L0_12 'ET' L0_34;
  913.  
  914. si ('EGA' graph 'O' );
  915. 'DESSIN' (eUMaxe 'ET' L0 ) 'TITRE' 'Ux moyen sur axe';
  916. 'DESSIN' (eUXaxe 'ET' L0) titr 'PROFIL Ux sur axe'
  917. tity 'Ux/U0' titx 'X/D';
  918. 'DESSIN' (eUya_34 'ET' L0_34) titr 'PROFIL Uy sur axe'
  919. tity 'Uy/U0' titx 'X/D';
  920. 'DESSIN' (eKMaxe 'ET' L0 ) 'TITRE' 'K moyen sur axe';
  921. 'DESSIN' (eKaxe 'ET' L0) titr 'PROFIL K sur axe'
  922. tity 'K/U0**2' titx 'X/D';
  923. finsi ;
  924.  
  925. Fin;
  926.  
  927.  
  928.  
  929.  
  930.  
  931.  
  932.  
  933.  
  934.  

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