Télécharger Comte-Bellot.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : Comte-Bellot.dgibi
  2. ************************************************************************************
  3. * *
  4. * Ecoulement entre deux plaques planes infinies (Symétrie) *
  5. * *
  6. * Penser a augmenter YP qd Rey diminue *
  7. * Rey = 4.E6 -> YP = 1.e-3 Rey = 4.e4 -> YP=1.E2 *
  8. * Nu = 1.e-7 Rey=4.e6 *
  9. ************************************************************************************
  10. OPTI TRAC 'PSC' ;
  11. GRAPH = 'N' ;
  12. OPTION DIME 2 ELEM QUA8;
  13. COMPLET=VRAI;
  14. COMPLET=FAUX;
  15.  
  16. sqrt2 = 2.**0.5 ;
  17. tole = 1e-7;
  18.  
  19. ******************************************************************
  20. * Paramètres géométriques du sous-canal
  21. ******************************************************************
  22.  
  23. * Demi-espacement entre les deux plaques
  24. h = 0.25;
  25.  
  26.  
  27. * Longueur de l'espace simulé (20 fois le diamètre hydraulique)
  28. Si COMPLET;
  29. L = 20;
  30. Sinon;
  31. L = 10;
  32. Finsi;
  33. ******************************************************************
  34. * Paramètres du maillage
  35. ******************************************************************
  36.  
  37. * nombre de points sur la verticale
  38. nv0 = 15 ;
  39.  
  40. * nombre de points sur la longueur
  41. nh0 = 35;
  42. nh0 = 70;
  43.  
  44.  
  45.  
  46. ******************************************************************
  47. * Paramètres lié à l'écoulement
  48. ******************************************************************
  49.  
  50. * Constante du modèle de turbulence
  51. * CMU = 0.09 ;
  52.  
  53. * Viscosité cinématique du fluide
  54. NU = 1.e-7;
  55.  
  56. * Diamètre hydraulique du sous-canal (DH = 1 car adim)
  57. DH = 4.*h ;
  58.  
  59. * Vitesse initiale
  60. UI = 0.4626;
  61.  
  62. * Constante cnu
  63. cnu = 0.09 ;
  64.  
  65. ******************************************************************
  66. * Paramètres liés à la résolution numérique
  67. ******************************************************************
  68.  
  69. KPRESS = 'CENTREP1' ;
  70. DISCR = 'MACRO';
  71.  
  72. * Nombre d'itérations
  73. * Pas de temps
  74. Si Complet;
  75. NITMA =800;
  76. DT = 0.1 ;
  77. Sinon;
  78. NITMA =10;
  79. DT = 2. ;
  80. Finsi ;
  81.  
  82.  
  83. ******************************************************************
  84. * Construction du maillage
  85. ******************************************************************
  86.  
  87. * Points pour construire le maillage élémentaire
  88. p0 = 0. 0. ;
  89. p1 = L 0. ;
  90. p2 = L h ;
  91. p3 = 0. h ;
  92.  
  93.  
  94.  
  95. * Construction des cotés élémentaires
  96. bsym = DROI nh0 p0 p1 ;
  97. sortie = DROI (-1*nv0) p1 p2 'DINI' 0.15 'DFIN' 0.1 ;
  98. bwall = DROI nh0 p2 p3 ;
  99. cot34 = DROI (-1*nv0) p0 p3 'DINI' 0.15 'DFIN' 0.1 ;
  100. entree = INVE cot34 ;
  101.  
  102.  
  103. dom = daller bsym sortie bwall entree ;
  104. *tracer dom;
  105.  
  106. * Description du domaine de calcul et des frontières
  107.  
  108.  
  109. NN = 'NBEL' dom;
  110. 'MESSAGE' 'Nombre d éléments du maillage' NN;
  111.  
  112.  
  113.  
  114. ********************************************************************
  115. * Changement des éléments du maillage en DISCR:
  116. ********************************************************************
  117.  
  118. £dom = 'CHANGER' dom QUAF ;
  119. £entree = 'CHANGER' entree QUAF ;
  120. £sortie = 'CHANGER' sortie QUAF ;
  121. £bsym = 'CHANGER' bsym QUAF ;
  122. £bwall = 'CHANGER' bwall QUAF ;
  123.  
  124.  
  125. elim tole (£dom et £entree et £sortie et £bsym et £bwall) ;
  126.  
  127.  
  128. **************************************************************
  129. * Formulation du domaine Navier Stokes:
  130. **************************************************************
  131.  
  132. $dom = 'MODE' £dom 'NAVIER_STOKES' DISCR ;
  133. $entree = 'MODE' £entree 'NAVIER_STOKES' DISCR ;
  134. $sortie = 'MODE' £sortie 'NAVIER_STOKES' DISCR ;
  135. $bsym = 'MODE' £bsym 'NAVIER_STOKES' DISCR ;
  136. $bwall = 'MODE' £bwall 'NAVIER_STOKES' DISCR ;
  137.  
  138. DOMA $dom 'IMPR' ;
  139. entree = 'DOMA' $entree 'MAILLAGE' ;
  140. sortie = 'DOMA' $sortie 'MAILLAGE' ;
  141. bwall = 'DOMA' $bwall 'MAILLAGE' ;
  142. bsym = 'DOMA' $bsym 'MAILLAGE' ;
  143.  
  144.  
  145. *****************************************************************
  146. * Création de la table des inconnues et initialisation:
  147. *****************************************************************
  148. * Rey = 4.E6 -> YP = 1.e-3 Rey = 4.e4 -> YP=1.E2
  149. YP = 1.e-2 ;
  150. YP = 1.e-3 ;
  151.  
  152. KE = (0.02 * UI)*(0.02 * UI)/2.;
  153. EE = (KE**1.5)/h ;
  154.  
  155.  
  156. * Vitesse de frottement :
  157.  
  158. * Viscosité turbulente :
  159. NUT = 0.09 * KE * KE / EE ; Rey=UI*DH/NU;
  160. 'MESSAGE' 'NUT =' NUT ' Rey=' Rey;
  161.  
  162.  
  163. RV = EQEX 'ITMA' NITMA 'OPTI' 'EF' 'IMPL' 'SUPG'
  164. 'ZONE' $bwall 'OPER' 'FPU' 1. 'UN' NU 'UET' YP 'INCO' 'UN' 'KN' 'EN'
  165. 'ZONE' $dom 'OPER' 'KEPSILON' 1. 'UN' NU 'DT' 'INCO' 'KN' 'EN'
  166. 'ZONE' $dom 'OPER' 'NS' 1. 'UN' 'MUF' 'INCO' 'UN'
  167. ;
  168.  
  169.  
  170.  
  171. * Condition limite sur le mur
  172. *RV = EQEX RV 'CLIM' 'UN' 'UIMP' bwall 0. 'UN' 'VIMP' bwall 0.;
  173.  
  174. * Condition limite de symétrie
  175. RV = EQEX RV 'CLIM' 'UN' 'VIMP' bsym 0.;
  176.  
  177.  
  178. * Condition limite sur les frontières libres
  179. * Remarque : pas besoin d'imposer une condition à la sortie car en utilisant la fonction 'CENTREP1'
  180. * pour la pression, si aucune valeur n'est précisé, son intégrale sur cette frontière est maintenue nulle.
  181. RV = EQEX RV 'CLIM' 'UN' 'UIMP' entree UI 'UN' 'VIMP' entree 0.;
  182. RV = EQEX RV 'CLIM' 'KN' 'TIMP' entree KE;
  183. RV = EQEX RV 'CLIM' 'EN' 'TIMP' entree EE;
  184.  
  185.  
  186. RV = EQEX RV 'OPTI' 'EFM1' 'CENTREE'
  187. 'ZONE' $dom 'OPER' DFDT 1. 'UN' DT 'INCO' 'UN';
  188.  
  189.  
  190. RVP = EQEX 'OPTI' 'EF' KPRESS
  191. 'ZONE' $dom 'OPER' 'KBBT' -1. 'INCO' 'UN' 'PRES'
  192. 'ZONE' $bwall 'OPER' 'VNIMP' $dom 0. 'INCO' 'UN' 'PRES';
  193.  
  194.  
  195. RV.'PROJ'= RVP ;
  196.  
  197.  
  198. RV.INCO= 'TABLE' INCO;
  199. RV.INCO.'UN' = 'KCHT' $dom 'VECT' 'SOMMET' (UI 0.);
  200. RV.INCO.'UNM1' = 'KCHT' $dom 'VECT' 'SOMMET' (0.D0 0.D0);
  201. RV.INCO.'KN' = KCHT $dom SCAL SOMMET KE;
  202. RV.INCO.'EN' = KCHT $dom SCAL SOMMET EE;
  203. RV.INCO.'PRES' = 'KCHT' $dom 'SCAL' KPRESS 0.;
  204. RV.INCO.'DT' = DT ;
  205. *RV.INCO.'UET' = KCHT $bwall 'SCAL' 'SOMMET' 2.E-2;
  206.  
  207. rv.inco.'MUF' = 200.*NU;
  208.  
  209. temp place;
  210. exec rv ;
  211. temp place;
  212.  
  213.  
  214. * RESULTATS
  215.  
  216. un= rv.inco.'UN' ;
  217. kn = rv.inco.'KN' ;
  218. en = rv.inco.'EN' ;
  219. prodt = rv.inco.'PRODT' ;
  220. nut =(rv.inco.'MUF')*(1./NU) ;
  221. pres = ELNO $dom (exco 'PRES' rv.inco.'PRESSION') KPRESS;
  222. nun = (pscal un (mots 'UX' 'UY') un (mots 'UX' 'UY'))**0.5 ;
  223.  
  224. USORTIE = 'EVOL' 'CHPO' un 'UX' sortie;
  225. KSORTIE = EVOL 'CHPO' RV.INCO.'KN' sortie;
  226. ESORTIE = EVOL 'CHPO' RV.INCO.'EN' sortie;
  227. NSORTIE = EVOL 'CHPO' nut sortie;
  228. PTSORTIE = EVOL 'CHPO' prodt sortie;
  229.  
  230. YPLUSI = YP * (rv.inco.'UET') / NU;
  231. evyplus=evol chpo yplusi (inve bwall);
  232. YSORTIE = COOR 2 sortie;
  233.  
  234. * Vitesse de frottement à la sortie du domaine (au point p2 plus précisément) :
  235. UFROT = 'EXTR' (rv.inco.'UET') 'SCAL' p2;
  236.  
  237. * Coordonnées yplus de la sortie
  238. YPLUS = ((0.25+YP)-YSORTIE) * UFROT / NU;
  239.  
  240. * Champ de vitesse normé
  241. UPLUS = USORTIE / UFROT;
  242.  
  243. Ksi= rv.inco.'MUF' * prodt *(inve rv.inco.'EN');
  244.  
  245. *--------- Référence calcul fin ----------------------------------------
  246. *-----------------------------------------------------------------------
  247. si faux;
  248. nx = extr nsortie 'ABSC' ;
  249. ny = extr nsortie 'ORDO' ;
  250. ky = extr ksortie 'ORDO' ;
  251. uy = extr usortie 'ORDO' ;
  252. xypl=extr evyplus 'ABSC' ;
  253. yypl=extr evyplus 'ORDO' ;
  254.  
  255. list nx;
  256. list ny;
  257. list ky;
  258. list uy;
  259. list xypl;
  260. list yypl;
  261. finsi;
  262.  
  263. xnuts=prog
  264. 0.0000 2.38164E-02 4.76327E-02 6.71384E-02 8.66440E-02
  265. 0.10262 0.11859 0.13168 0.14476 0.15548
  266. 0.16619 0.17497 0.18375 0.19093 0.19812
  267. 0.20401 0.20989 0.21472 0.21954 0.22348
  268. 0.22743 0.23067 0.23390 0.23655 0.23920
  269. 0.24137 0.24354 0.24531 0.24709 0.24854
  270. 0.25000 ;
  271. ynuts=prog
  272. 260.31 366.25 901.62 1383.4 1818.7
  273. 2121.0 2354.5 2485.6 2555.9 2565.4
  274. 2532.7 2464.9 2358.3 2239.4 2091.6
  275. 1947.3 1781.9 1629.6 1462.3 1313.2
  276. 1153.8 1015.0 869.28 744.33 615.34
  277. 506.32 395.64 303.88 212.84 139.97
  278. 63.065 ;
  279. yky=prog
  280. 1.08988E-05 2.03752E-05 6.73979E-05 1.24382E-04 1.88152E-04
  281. 2.44708E-04 3.01914E-04 3.49488E-04 3.96502E-04 4.35483E-04
  282. 4.75142E-04 5.07189E-04 5.38433E-04 5.63879E-04 5.88796E-04
  283. 6.09357E-04 6.29581E-04 6.46532E-04 6.63290E-04 6.77661E-04
  284. 6.92061E-04 7.04813E-04 7.17758E-04 7.29813E-04 7.42301E-04
  285. 7.55025E-04 7.68504E-04 7.84668E-04 8.00892E-04 8.24540E-04
  286. 7.63839E-04 ;
  287. yuy=prog
  288. 0.50623 0.50569 0.50187 0.49693 0.49080
  289. 0.48504 0.47862 0.47285 0.46655 0.46091
  290. 0.45477 0.44928 0.44327 0.43787 0.43190
  291. 0.42647 0.42043 0.41488 0.40862 0.40279
  292. 0.39613 0.38983 0.38249 0.37538 0.36686
  293. 0.35831 0.34760 0.33619 0.32065 0.30182
  294. 0.26896 ;
  295. xypl=prog
  296. 0.0000 0.14286 0.28571 0.42857 0.57143
  297. 0.71429 0.85714 1.0000 1.1429 1.2857
  298. 1.4286 1.5714 1.7143 1.8571 2.0000
  299. 2.1429 2.2857 2.4286 2.5714 2.7143
  300. 2.8571 3.0000 3.1429 3.2857 3.4286
  301. 3.5714 3.7143 3.8571 4.0000 4.1429
  302. 4.2857 4.4286 4.5714 4.7143 4.8571;
  303.  
  304. xypl=xypl et (prog
  305. 5.0000 5.1429 5.2857 5.4286 5.5714
  306. 5.7143 5.8571 6.0000 6.1429 6.2857
  307. 6.4286 6.5714 6.7143 6.8571 7.0000
  308. 7.1429 7.2857 7.4286 7.5714 7.7143
  309. 7.8571 8.0000 8.1429 8.2857 8.4286
  310. 8.5714 8.7143 8.8571 9.0000 9.1429
  311. 9.2857 9.4286 9.5714 9.7143 9.8571);
  312.  
  313. xypl=xypl et (prog
  314. 10.0000 10.143 10.286 10.429 10.571
  315. 10.714 10.857 11.000 11.143 11.286
  316. 11.429 11.571 11.714 11.857 12.000
  317. 12.143 12.286 12.429 12.571 12.714
  318. 12.857 13.000 13.143 13.286 13.429
  319. 13.571 13.714 13.857 14.000 14.143
  320. 14.286 14.429 14.571 14.714 14.857);
  321.  
  322.  
  323. xypl=xypl et (prog
  324. 15.000 15.143 15.286 15.429 15.571
  325. 15.714 15.857 16.000 16.143 16.286
  326. 16.429 16.571 16.714 16.857 17.000
  327. 17.143 17.286 17.429 17.571 17.714
  328. 17.857 18.000 18.143 18.286 18.429
  329. 18.571 18.714 18.857 19.000 19.143
  330. 19.286 19.429 19.571 19.714 19.857 20.000) ;
  331.  
  332. yypl=prog
  333. 227.14 200.86 186.99 180.11 175.86
  334. 172.94 170.42 169.16 167.07 166.45
  335. 164.82 164.58 163.12 162.93 161.80
  336. 161.67 160.69 160.50 159.77 159.58
  337. 158.94 158.71 158.22 158.01 157.57
  338. 157.33 156.99 156.79 156.46 156.25
  339. 155.99 155.81 155.56 155.37 155.17;
  340.  
  341. yypl=yypl et (prog
  342. 155.01 154.81 154.65 154.49 154.35
  343. 154.19 154.05 153.92 153.80 153.67
  344. 153.55 153.44 153.34 153.22 153.13
  345. 153.03 152.95 152.85 152.77 152.68
  346. 152.62 152.53 152.47 152.39 152.34
  347. 152.26 152.21 152.14 152.10 152.03
  348. 151.99 151.93 151.90 151.84 151.80);
  349.  
  350. yypl=yypl et (prog
  351. 151.75 151.73 151.67 151.65 151.60
  352. 151.58 151.54 151.52 151.48 151.46
  353. 151.43 151.41 151.38 151.37 151.33
  354. 151.33 151.30 151.29 151.26 151.26
  355. 151.23 151.23 151.20 151.21 151.18
  356. 151.19 151.16 151.17 151.14 151.16
  357. 151.13 151.15 151.12 151.14 151.12);
  358.  
  359. yypl=yypl et (prog
  360. 151.13 151.11 151.13 151.12 151.13
  361. 151.12 151.13 151.13 151.14 151.14
  362. 151.14 151.15 151.15 151.17 151.15
  363. 151.19 151.16 151.21 151.17 151.23
  364. 151.18 151.26 151.19 151.29 151.21
  365. 151.33 151.23 151.36 151.26 151.39
  366. 151.30 151.41 151.35 151.42 151.42 151.38) ;
  367.  
  368. evyplusr=evol manu xypl yypl ;
  369. evnutsr=evol manu xnuts ynuts ;
  370. evuyr=evol manu xnuts yuy ;
  371. evkyr=evol manu xnuts yky ;
  372. *-----------------------------------------------------------------------
  373. ERR=0;
  374. TAB1=TABLE;
  375. TAB1.'TITRE'= table ;
  376. TAB1.1 =' MARQ CROI REGU';
  377. TAB1.2 ='TIRR MARQ TRIA REGU';
  378. TAB1.'TITRE' . 1 ='Référence';
  379. TAB1.'TITRE' . 2 ='Calcul grossier';
  380.  
  381. m=(INTG evyplusr 'ABSO') SOMM ;
  382. dyplus=(SOMM (INTG (evyplus - evyplusr) 'ABSO'))/m;
  383. mess ' YPLUS : m=' m ' dyplus=' dyplus;
  384. Si (dyplus > 0.07);ERR=ERR + 1;Finsi;
  385.  
  386. m=(INTG evuyr 'ABSO') SOMM ;
  387. duyr=(SOMM (INTG (evuyr - usortie) 'ABSO'))/m;
  388. mess ' Usortie : m=' m ' duyr=' duyr;
  389. Si (duyr > 0.07);ERR=ERR + 1;Finsi;
  390.  
  391. m=(INTG evnutsr 'ABSO') SOMM ;
  392. dnut=(SOMM (INTG (evnutsr - nsortie) 'ABSO'))/m;
  393. mess ' Nu sortie : m=' m ' dnut=' dnut;
  394. Si (dnut > 0.3);ERR=ERR + 1;Finsi;
  395.  
  396. m=(INTG evkyr 'ABSO') SOMM ;
  397. dks=(SOMM (INTG (evkyr - ksortie) 'ABSO'))/m;
  398. mess ' K sortie : m=' m ' dks=' dks;
  399. Si (dks > 0.3);ERR=ERR + 1;Finsi;
  400.  
  401.  
  402. Si (EGA GRAPH 'O') ;
  403. DESS USORTIE TITR ' USORTIE ' ;
  404. DESS evyplus TITR 'Y+';
  405. DESS KSORTIE TITR ' KSORTIE ' ;
  406. DESS ESORTIE TITR ' ESORTIE ' ;
  407. DESS NSORTIE TITR ' NUTSORTIE / NU ' ;
  408. DESS PTSORTIE TITR ' ProdtSORTIE ' ;
  409.  
  410. trace kn dom (cont dom) TITR 'KN' ;
  411. trace en dom (cont dom) TITR 'EN' ;
  412. ung = vect un 0.1 ux uy jaune ;
  413. trace ung dom;
  414. trace rv.inco.'TKTE' dom (cont dom) TITR 'TKTE';
  415. trace nut dom (cont dom) TITR 'nut/nu ' ;
  416. trace Ksi dom (cont dom) TITR ' Ksi';
  417. trace prodt dom (cont dom) TITR ' P';
  418.  
  419. dess (evyplusr et evyplus) TITR ' Y+' LEGE TAB1 'GRIL';
  420. dess (evuyr et usortie) TITR 'Usortie' LEGE TAB1 'GRIL';
  421. dess (evnutsr et nsortie) TITR 'Nut/nu sortie' LEGE TAB1 'GRIL';
  422. dess (evkyr et ksortie) TITR 'Ksortie' LEGE TAB1 'GRIL';
  423. FINSI;
  424.  
  425. Si (ERR > 0);
  426. mess 'NOMBRE d ERREURS ' ERR;
  427. ERREUR 5;
  428. FINSI;
  429. FIN ;
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  

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