Télécharger la plaquette Cast3M


fr en

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. * CHANGEMENT DE SUPPORT pour les operations
  374. absci   = extr evyplus  'ABSC' ;
  375. ordon   ='IPOL' absci evyplusr ;
  376. evyplusr= evol manu absci ordon;
  377.  
  378. absci   = extr usortie  'ABSC' ;
  379. ordon   ='IPOL' absci evuyr ;
  380. evuyr   = evol manu absci ordon;
  381.  
  382. absci   = extr nsortie  'ABSC' ;
  383. ordon   ='IPOL' absci evnutsr ;
  384. evnutsr = evol manu absci ordon;
  385.  
  386. absci   = extr ksortie  'ABSC' ;
  387. ordon   ='IPOL' absci evkyr ;
  388. evkyr   = evol manu absci ordon;
  389.  
  390. *-----------------------------------------------------------------------
  391. TAB1='TABL';
  392. TAB1.'TITRE'= table ;
  393. TAB1.1 ='     MARQ CROI REGU';
  394. TAB1.2 ='TIRR MARQ TRIA REGU';
  395. TAB1.'TITRE' . 1 ='Référence';
  396. TAB1.'TITRE' . 2 ='Calcul grossier';
  397.  
  398. m=(INTG evyplusr 'ABSO') ;  
  399. dyplus=((INTG (evyplusr - evyplus) 'ABSO'))/m;
  400. mess ' YPLUS : m=' m ' dyplus=' dyplus;
  401. SI (dyplus > 0.016);
  402.  'MESS' 'ERREUR calcul dyplus ' dyplus;
  403.   ERREUR 5;
  404. FINSI;
  405. dess (evyplus et (evyplusr coul bleu));
  406.  
  407. m=(INTG evuyr 'ABSO') ;
  408. duyr=((INTG (evuyr - usortie) 'ABSO'))/m;
  409. mess ' Usortie : m=' m ' duyr=' duyr;
  410. SI (duyr > 0.032);
  411.  'MESS' 'ERREUR calcul duyr ' duyr;
  412.   ERREUR 5;
  413. FINSI;
  414. dess (usortie et (evuyr coul bleu));
  415.  
  416. m=(INTG evnutsr 'ABSO') ;
  417. dnut=((INTG (evnutsr - nsortie) 'ABSO'))/m;
  418. mess ' Nu sortie : m=' m ' dnut=' dnut;
  419. SI (dnut > 0.75);
  420.  'MESS' 'ERREUR calcul dnut ' dnut;
  421.   ERREUR 5;
  422. FINSI;
  423. dess (nsortie et (evnutsr coul bleu));
  424.  
  425. m=(INTG evkyr 'ABSO') ;
  426. dks=((INTG (evkyr - ksortie) 'ABSO'))/m;
  427. mess ' K sortie : m=' m ' dks=' dks;
  428. SI (dks > 0.65);
  429.  'MESS' 'ERREUR calcul dks ' dks;
  430.   ERREUR 5;
  431. FINSI;
  432. dess (ksortie et (evkyr coul bleu));
  433.  
  434. Si (EGA GRAPH 'O') ;
  435.   DESS USORTIE TITR ' USORTIE ' ;
  436.   DESS evyplus TITR 'Y+';
  437.   DESS KSORTIE TITR ' KSORTIE ' ;
  438.   DESS ESORTIE TITR ' ESORTIE ' ;
  439.   DESS NSORTIE TITR ' NUTSORTIE / NU ' ;
  440.   DESS PTSORTIE TITR ' ProdtSORTIE ' ;
  441.  
  442.   trace kn dom (cont dom) TITR 'KN' ;
  443.   trace en dom (cont dom) TITR 'EN' ;
  444.   ung = vect un 0.1 ux uy jaune ;
  445.   trace ung dom;
  446.   trace rv.inco.'TKTE' dom (cont dom) TITR 'TKTE';
  447.   trace nut dom (cont dom) TITR 'nut/nu ' ;
  448.   trace Ksi dom (cont dom) TITR ' Ksi';
  449.   trace prodt dom (cont dom) TITR ' P';
  450.  
  451.   dess (evyplusr et evyplus) TITR ' Y+' LEGE TAB1 'GRIL';
  452.   dess (evuyr et usortie) TITR 'Usortie' LEGE TAB1 'GRIL';
  453.   dess (evnutsr et nsortie) TITR 'Nut/nu sortie' LEGE TAB1 'GRIL';
  454.   dess (evkyr et ksortie) TITR 'Ksortie' LEGE TAB1 'GRIL';
  455. FINSI;
  456.  
  457. FIN ;
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  

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