Télécharger la plaquette Cast3M


fr en

Télécharger canal-Chien.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : canal-chien.dgibi
  2. *
  3. * Ce cas teste le modèle K-epsilon Bas Reynolds de Chien
  4. * sur l'écoulement turbulent dans un canal plan.
  5.  opti trace 'PSC';
  6.  COMPLET = FAUX ;
  7.  GRAPH   = FAUX ;
  8. *COMPLET = VRAI ;
  9. *GRAPH   = VRAI ;
  10.  
  11. *-----------------------------------------------------------------------
  12.  opti dime 2 elem TRI6;
  13.  opti dime 2 elem QUA8;
  14. *-----------------------------------------------------------------------
  15. *      |Triangles | Quadrangles
  16. * ---------------------------
  17. * LINE |          | OK
  18. * LINB |          |
  19. * MACRO|   ?      | OK
  20. * QUAF |          | OK
  21. *-----------------------------------------------------------------------
  22.  
  23.  DISCR=MACRO;
  24.  KPRES=CENTREP1;
  25.  KSUPG=SUPG   ;
  26.  
  27. ******=================================================**********
  28. ******=================================================**********
  29.  H=1. ;
  30.  
  31. Si COMPLET ;
  32.  NBIT =900;
  33.  KGEO=3;
  34.  CFLmax=5.;
  35.  DT  = 0.01;
  36.  L=30.*H;
  37. Sinon;
  38.  NBIT = 40;
  39.  KGEO=1;
  40.  CFLmax=20.;
  41.  DT  = 0.15;
  42.  L=3.*H;
  43. Finsi;
  44.  
  45. *--------------------------- procédure ---------------------------------
  46. * Cette procédure calcule à chaque pas de temps le Y+ à partir
  47. * du gradient de vitesse en paroi
  48.  'DEBPROC' CALYPLUS ;
  49.   argu rx*table;
  50. *mess 'DEBUT CALYPLUS';
  51.  yp=rv.inco.'dparoi';
  52.  un=rv.inco.'UN';
  53.  uux=exco un 'UX';
  54.  gru=((Mu*(kops 'GRADS' uux $mt)*(1./Ro)) ABS)**0.5;
  55.  uet1=evol chpo (redu (exco gru 'UY') par1) par1;
  56.  uet2=evol chpo (redu (exco gru 'UY') par2i) par2i;
  57.  xx1=extr uet1 'ABSC';fx1=extr uet1 'ORDO';
  58.  xx2=extr uet2 'ABSC';fx2=extr uet2 'ORDO';
  59.  yplus1=ipol xx1 fx1 (coor 1 mt) ;
  60.  yplus2=ipol xx2 fx2 (coor 1 mt) ;
  61.  
  62.  iy=masq (coor 2 mt) 'SUPERIEUR' (H/2.);
  63.  yplus=(iy*yplus2)+((1.-iy)*yplus1);
  64.  yplus=yplus*yp*(1./nu) ABS ;
  65.  yplus=1.*yplus + 1.e-3;
  66.  rv.inco.'yplus'=yplus;
  67. *trace yplus mt (cont mt) TITR ' Y+' 'NCLK';
  68. * nutsnu=(rv.inco.'MUF')*(1./NU);
  69. * trace nutsnu mt (cont mt) TITR ' Nutsnu' 'NCLK';
  70.  
  71.   nun=psca un (mots 'UX' 'UY')  un (mots 'UX' 'UY');
  72.   nun=nun**0.5;
  73.   Diamin =elno (doma $mt 'DIAMIN') $mt ;
  74.   DTc=Mini (CFLmax*diamin*(inve nun));
  75.  
  76. mess ' dtc' ;  list dtc;
  77.  mess ' dt ' ; list dt;
  78.  
  79.  
  80. * CFLmax=DT * nun * (inve diamin);
  81. * trace CFLmax mt (cont mt) TITR ' CFLmax' 'NCLK';
  82.   rv.inco.'DT'=mini (prog DTc rv.inco.'DTINI');
  83.   mess ' Cflmax=' CFLmax ' DTc=' DTc 'Pdt retenu ' rv.inco.'DT';
  84.  
  85.  st1 mat1 = KOPS 'MATRIK' ;
  86. *mess 'FIN   CALYPLUS';
  87.  FINPROC st1 mat1;
  88. *----------------------- fin procédure ---------------------------------
  89.  
  90. *----------------------- Maillage --------------------------------------
  91. epsi = 1.e-5;
  92.  
  93.  Si(EGA KGEO 1);
  94. *dcl=H/200;
  95. *dm=H/20;
  96. *nm=15;
  97. *dcx=H/12.;
  98. *dmx=H/10.;
  99.  
  100.  dcl=H/130;
  101.  dm=H/13;
  102.  nm=13;
  103.  dcx=H/10.;
  104.  dmx=H/8.;
  105.  Finsi;
  106.  
  107.  Si(EGA KGEO 2);
  108.  dcl=H/300;
  109.  dm=H/30;
  110.  nm=20;
  111.  dcx=H/24.;
  112.  dmx=H/20.;
  113.  Finsi;
  114.  
  115.  Si(EGA KGEO 3);
  116.  dcl=H/400;
  117.  dm=H/40;
  118.  nm=25;
  119.  dcx=H/36.;
  120.  dmx=H/30.;
  121.  Finsi;
  122.  
  123. P0 = 0. 0. ;
  124. P1 = L  0. ;
  125. P2 = L H   ;
  126. P3 = 0. H  ;
  127. PS3 = (P0 PLUS P3)*0.5;
  128. PS31= (P0 PLUS P3)*0.2;
  129. PS32= (P3 MOINS PS31);
  130. entree= P0 d dini dcl dfin dm ps31 d nm ps32
  131.  d dini dm dfin dcl P3;
  132. sortie=(entree plus (L 0.)) inve;
  133. par1  = p0 d dini dcx dfin dmx p1 ;
  134. par2  = (par1 plus (0. H)) inve;
  135.  
  136. pc= 1. 1.;
  137. ang= 43.;
  138. ang= 60.;
  139.  ang=  0. ;
  140. entree=entree TOUR ang  pC;
  141. sortie=sortie TOUR  ang pC;
  142. par1  =par1   TOUR  ang pC;
  143. par2  =par2   TOUR  ang pC;
  144.  
  145. elim (entree et par1 et par2 et sortie) epsi;
  146. mt = daller entree par2 sortie par1;
  147. v1=1. 0. ;
  148. v2=0. -1.;
  149. v3=-1. 0.;
  150. v4=0. 1. ;
  151. vt1= tour v1  ang p0;
  152. vt2= tour v2  ang p0;
  153. vt3= tour v3  ang p0;
  154. vt4= tour v4  ang p0;
  155. list vt1;
  156. list vt2;
  157. list vt3;
  158. list vt4;
  159.  
  160. Mmt    = chan mt     QUAF ;
  161. Mentree= chan entree QUAF ;
  162. Msortie= chan sortie QUAF ;
  163. Mpar1  = chan par1   QUAF ;
  164. Mpar2  = chan par2   QUAF ;
  165.  
  166. elim (Mmt et Mentree et Msortie et Mpar1 et Mpar2) epsi ;
  167.  
  168. $mt     = model Mmt     'NAVIER_STOKES' DISCR;
  169. $entree = model Mentree 'NAVIER_STOKES' DISCR;
  170. $sortie = model Msortie 'NAVIER_STOKES' DISCR;
  171. $par1   = model Mpar1   'NAVIER_STOKES' DISCR;
  172. $par2   = model Mpar2   'NAVIER_STOKES' DISCR;
  173.  
  174. doma $mt 'IMPR';
  175. doma $entree 'IMPR';
  176.  
  177. mt     = doma $mt     maillage;
  178. entree = doma $entree maillage;
  179. sortie = doma $sortie maillage;
  180. par1   = doma $par1   maillage;
  181. par2   = doma $par2   maillage;
  182. par2i  = inve par2;
  183. axe    = par1 plus (0. (H/2.));
  184. elim (axe et mt) epsi;
  185.   entre1= chan entree poi1;
  186.   entre1= elem entre1 (lect 2 pas 1 ((nbel entre1) - 1 ));
  187.  
  188. dparoi=(coor 2 mt);
  189. dparoi=dparoi*(H - dparoi);
  190. dparoi=dparoi+1.e-10;
  191.  
  192. Si Graph;
  193. trace mt TITR ' Maillage ';
  194. dess (evol chpo dparoi sortie) TITR ' Distance à la paroi';
  195.  
  196. * On vérifie l'allongement
  197. diamin=doma $mt 'DIAMIN';
  198. diamax=doma $mt 'DIAMAX';
  199. ksi=diamax*(inve diamin);
  200. ksi=elno $mt ksi centre;
  201. trace ksi mt (cont mt) TITR ' Facteur d allongemet du maillage';
  202. Finsi;
  203.  
  204. *------------------- Fin Maillage --------------------------------------
  205.  
  206. * on impose un profil de vitesse en entrée soit plat (nn=1.e3) soit
  207. * turbulent (nn=7).
  208. nn=1.e3;
  209. nn=7.;
  210. * La vitesse débitante est de 1.
  211. U0  = 1. ;
  212.  
  213. tasser entree;
  214. ye=coor 2 entree;
  215. Dgj=doma $entree 'XXDIAGSI';
  216. a=(SOMT (Dgj*((ye*(H - ye))**(1./nn ))));
  217. Ufr=U0/a*((ye*(H - ye))**(1./nn ));
  218. Si Graph;
  219. dess (evol chpo ufr entree) TITR 'Profil de vitesse en entrée';
  220. Finsi;
  221.  
  222. uux=(cos ang) * U0;
  223. uuy=(sin ang) * U0;
  224. KE=((0.02*Ufr)**2.);
  225. EE=KE**1.5/H;
  226.  
  227. MU  = 3.14E-5;
  228. Ro  = 1.;
  229. YP  = 1.e-2 ;
  230. Nu  = Mu /Ro ;
  231.  
  232. Rey= U0 *H * 2. / NU;
  233. mess ' Reynolds = ' Rey;
  234.  
  235. RV= EQEX  'OMEGA' 1.    'NITER' 1  ITMA NBIT    'FIDT' 1
  236. 'OPTI' 'EF' 'IMPL'   KSUPG  KPRES
  237.   ZONE  $mt  'OPER' 'CALYPLUS'
  238.   ZONE  $mt  'OPER' 'KEPSILON' 1. 'UN' 'NU' 'DT'   'INCO' 'KN' 'EN'
  239.   ZONE  $mt  'OPER' 'NS'       1. 'UN' 'MUF'     'INCO' 'UN'
  240. 'OPTI' 'EFM1' 'IMPL' 'CENTREE'
  241.   ZONE  $mt  'OPER'  DFDT 1. 'UNM' 'DT'  INCO 'UN';
  242.  
  243.  rv.'ALGO_KEPSILON'=mots 'Chien';
  244.  Mdl = chai 'Chien' ;
  245.  
  246.   RV= EQEX RV  CLIM
  247.   UN UIMP entre1 Ufr  UN VIMP entree (uuy + 1.e-4)
  248.   UN UIMP (par1 et par2) 0. UN VIMP (par1 et par2) 0.
  249.   KN TIMP entre1 KE  EN TIMP entre1 EE
  250.   KN TIMP (par1 et par2) 0. EN TIMP (par1 et par2) 0.
  251.   ;
  252.  
  253.  rv.'METHINV'.TYPINV=3 ;
  254.  rv.'METHINV'.IMPINV=0 ;
  255.  rv.'METHINV'.NITMAX=700;
  256.  rv.'METHINV'.PRECOND=3 ;
  257.  rv.'METHINV'.ILUTLFIL=4 ;
  258.  rv.'METHINV'.ILUTDTOL=0.;
  259.  rv.'METHINV'.RESID  =1.e-8;
  260.  rv.'METHINV' . 'FCPRECT'=1   ;
  261.  rv.'METHINV' . 'FCPRECI'=1   ;
  262.  
  263.     RVP= EQEX
  264.     'OPTI' 'EF' KPRES
  265.     ZONE  $mt    'OPER'  KBBT  (-1.) 'INCO' 'UN' 'PRES'
  266.     ;
  267.  
  268.     rvp.'METHINV'.TYPINV=1        ;
  269.     rvp.'METHINV'.IMPINV=0 ;
  270.     rvp.'METHINV'.NITMAX=5000;
  271.     rvp.'METHINV'.PRECOND=3 ;
  272.     rvp.'METHINV'.ILUTLFIL=10;
  273.     rvp.'METHINV'.ILUTDTOL=0.;
  274.     rvp.'METHINV'.RESID  =1.e-10 ;
  275.     rvp.'METHINV' . 'FCPRECT'=50000 ;
  276.     rvp.'METHINV' . 'FCPRECI'=50000 ;
  277.  
  278.  
  279.     RV.'PROJ'= RVP ;
  280.  
  281. rv.'INCO'= table inco ;
  282. rv.'INCO'.'UN'  =  kcht $mt  vect sommet (1.e-5 1.e-5) ;
  283. rv.'INCO'.'KN'  =  kcht $mt  scal sommet 1.e-5 ;
  284. rv.'INCO'.'EN'  =  kcht $mt  scal sommet 1.e-5 ;
  285. rv.'INCO'.'dparoi'=kcht $mt  scal sommet dparoi ;
  286. rv.'INCO'.'UNM' =  kcht $mt  vect sommet (1.e-5 1.e-5) ;
  287. rv.'INCO'.'PRES'=  kcht $mt  scal kpres  0.;
  288. rv.'INCO'.'MUF'=kcht $mt scal sommet nu;
  289. rv.'INCO'.'NU'=NU;
  290. rv.'INCO'.'MU'=MU;
  291. rv.'INCO'.'DT'=DT;
  292. rv.'INCO'.'DTINI'=DT;
  293.  
  294.  rv.'LTPS' = PROG ;
  295.  rv.'Tps' = 1. ;
  296.  rv.'DT'  = DT;
  297.  rv.'NUPAT' =1;
  298.  
  299.  exec rv ;
  300.  
  301.  Si Complet
  302.  opti sauv 'canalKEChien_1';
  303.  sauv ;
  304.  opti sauv 'bidon' ;
  305.  
  306.  exec rv ;
  307.  opti sauv 'canalKEChien_2';
  308.  sauv ;
  309.  opti sauv 'bidon' ;
  310.  
  311.  exec rv ;
  312.  opti sauv 'canalKEChien_3';
  313.  sauv ;
  314.  opti sauv 'bidon' ;
  315.  
  316.  exec rv ;
  317.  opti sauv 'canalKEChien_4';
  318.  sauv ;
  319.  opti sauv 'bidon' ;
  320.  
  321.  exec rv ;
  322.  opti sauv 'canalKEChien_5';
  323.  sauv ;
  324.  opti sauv 'bidon' ;
  325.  
  326. * On purge la mémoire de tous les objets MATRIK pour avoir un petit
  327. * fichier de résultats
  328. rv.'TABRES'=0.;
  329. rv.'PROJ'=0.;
  330. rv.'METHINV'=0;
  331. rv.'METHINV2'=0;
  332. rv.'METHINV3'=0;
  333. rv.'METH_KEP'=0;
  334. rvp=0.;
  335.  
  336.  opti sauv 'canalKEChien';
  337.  sauv ;
  338.  Finsi;
  339.  
  340.  un = rv.inco.'UN' ;
  341.  qs=dbit un $sortie IMPR;
  342.  um=(1./H)*qs;
  343.  unn=(1./um)*un;
  344. *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  345. *                   POST-TRAITEMENT
  346. *-----------------------------------------------------------------------
  347. *--- Edition d'un profil de vitesse asymptotique de Référence ----------
  348. *--- Modèle de Buleev
  349. Debproc ANABU H*'FLOTTANT' Rey*'FLOTTANT' UM*'FLOTTANT' ;
  350. H2=H*H;
  351. p0=0. 0. ; p1=0. H. ;  Pm=0. (H/2.); dpar=H/Rey*0.05;  dctr=H*0.05;
  352. de= p0 d dini dpar dfin dctr pm d dini dctr dfin dpar p1 ;
  353. z= coor 2 de;
  354. Dh=2.*H;
  355. Rho=1.;
  356. mu=um/Rey*Rho*Dh;
  357. m=7.1;
  358. nu=mu/rho;
  359.  
  360. uet=um;
  361. REPETER CCP 15 ;
  362. fp=(uet/Um)**2.*8. ;
  363. F=fp/2/Dh*rho*(um*um);
  364. xi   = z*(1./H);
  365. xiet = (1. - ((1. - (nu*m/h/Uet*4))**0.5))/2.;
  366. u1=H2*F/2./mu*xi*(1. - xi);
  367. ax=xi*(1. - xi)/xiet/(1. - xiet) + 1.e-10;
  368. ax=m*m*F*H/Rho* (log ax);
  369. ax =(abs ax) + 1.e-10;
  370. ax0=(h*h*F/mu/2*xiet*(1. - xiet))**2.;
  371. u2=(ax0+ax)**0.5;
  372. iua=masq xi 'INFERIEUR' xiet;
  373. iub=masq xi 'SUPERIEUR' (1. - xiet);
  374. iu1=iua+iub;
  375.  
  376. uref = (iu1 * u1) + ((1. - iu1)*u2) ;
  377. evuref = evol chpo uref de ;
  378. umn=INTG evuref prog ;
  379. umn=(extr 1 umn)/H; Reyn=rho*umn*Dh/mu;
  380. w=0.9;
  381. uet=(w*(uet*um/umn)) + ((1. - w)*uet);
  382. dumn=abs ((umn - um)/um);
  383. mess 'Iteration ' &ccp '  Umn=' umn ' Reyn=' Reyn ' Uet=' Uet
  384.  ' dumn=' dumn;
  385. *dess evuref ;
  386. Si (dumn < 1.e-4) ; quitter CCP; Finsi ;
  387. Fin CCP;
  388. * Calcul des u+g et u+d
  389. yplus= (extr evuref 'ABSC')*uet/nu;
  390. evuplus=evuref*(1./uet);
  391. yy=extr evuplus ORDO ;
  392. lyz=lect (dime yy) pas -1 1;
  393. yyz=extr yy lyz;
  394. evuplud = (evol manu yplus yyz) coul verte;
  395. evuplug = (evol manu yplus yy) coul rouge ;
  396.  
  397. Resproc evuplug evuplud evuref uet;
  398. Finproc ;
  399.  
  400. DEBPROC KADER Pr*'FLOTTANT' Ret*'FLOTTANT' KOUL*'MOT';
  401. *--- Profil de Kader ---------------------------------------------------
  402.  
  403. yplusr= prog 0.1 pas 0.1 30 pas 1. 100. pas 50. 1.e3
  404.  pas 500. 1.e4 ;
  405. unp = prog (dime yplusr)*1.;
  406. ysr=yplusr/Ret;
  407. Gamma=((Pr*yplusr)**4.)/(unp+(5.*(Pr**3.)*yplusr));
  408. Gamma=Gamma*(-0.01);
  409. tplusr=Pr*yplusr*(exp Gamma);
  410. beta=((3.85*(Pr**(1./3.)) - 1.3)**2.);
  411. beta=beta+(2.12*(log Pr));
  412. aa1=1.5*(2.*unp - ysr);
  413. aa2=2.5*(2.*unp - ysr);
  414. bb1=(unp - ysr)*(unp - ysr)*2. + unp;
  415. bb2=(unp - ysr)*(unp - ysr)*4. + unp;
  416. aa1=aa1/bb1 ;
  417. aa2=aa2/bb2 ;
  418. cc1= ((2.12*(Log ((unp+yplusr)*aa1))) + (beta*unp))
  419.  *(exp(unp/Gamma));
  420. cc2= ((2.12*(Log ((unp+yplusr)*aa2))) + (beta*unp))
  421.  *(exp(unp/Gamma));
  422. tplusr1=tplusr + cc1;
  423. tplusr2=tplusr + cc2;
  424.  
  425. evt1=evol manu   yplusr tplusr1;
  426. evt2=evol manu   yplusr tplusr2;
  427. evt1=evt1 coul KOUL;
  428. evt2=evt2 coul KOUL;
  429. *--- Profil de Kader ---------------------------------------------------
  430. FINPROC evt1 evt2;
  431.  
  432. DEBPROC POST ;
  433.  tabk = table ;
  434. evupd1 evupg1 evuref1 uet1 = anabu 1. 2000. 1. ;
  435. evupd2 evupg2 evuref2 uet2 = anabu 1. 1.e4  1. ;
  436. evupd3 evupg3 evuref3 uet3 = anabu 1. 1.e5  1. ;
  437. evupd4 evupg4 evuref4 uet4 = anabu 1. 1.e6  1. ;
  438. evupd5 evupg5 evuref5 uet5 = anabu 1. 1.e7  1. ;
  439.  
  440. *--- Profil de Reichardt -----------------------------------------------
  441. Ka = 0.41;
  442. yplusr= prog 0.1 pas 0.1 30 pas 1. 100. pas 50. 1.e3
  443.  pas 500. 1.e4 pas 1.e4 2.e5 ;
  444. unp = prog (dime yplusr)*1.;
  445.  
  446. E=9.5;
  447. uplusr=(1./Ka)*(Log (unp + (Ka*yplusr)));
  448. uplusr=uplusr + (7.8* (unp - (exp (yplusr*(-1./11))))
  449.  - ( (yplusr*(-1./11)) * (exp (yplusr*(-1./3.))) ) );
  450.  
  451. evr1=evol manu   yplusr uplusr;
  452. evr1=evr1 coul verte;
  453. *--- Profil de Reichardt -----------------------------------------------
  454.  
  455. *--- Profil de Référence pour le KLbr ----------------------------------
  456.  
  457.  xrr=prog
  458.   2.99401E-06  1.24325E-03  4.78060E-03  1.05194E-02  1.83637E-02
  459.   2.82181E-02  3.99867E-02  5.35741E-02  6.88838E-02  8.58208E-02
  460.   0.10429      0.12419      0.14544      0.16793      0.19156
  461.   0.21625      0.24190      0.26841      0.29568      0.32363
  462.   0.35214      0.38114      0.41052      0.44019      0.47005
  463.   0.50000      0.52995      0.55981      0.58948      0.61886
  464.   0.64786      0.67637      0.70432      0.73159      0.75810 ;
  465.  xrr = xrr et (prog
  466.   0.78375      0.80844      0.83207      0.85456      0.87581
  467.   0.89571      0.91418      0.93112      0.94643      0.96001
  468.   0.97178      0.98163      0.98948      0.99522      0.99875
  469.   1.00000);
  470.  urr = prog
  471.    0.0000       9.36980E-02  0.35215      0.56346      0.68057
  472.    0.75493      0.80827      0.84978      0.88391      0.91302
  473.    0.93845      0.96103      0.98127      0.99955       1.0161
  474.    1.0311       1.0447       1.0569       1.0678       1.0774
  475.    1.0858       1.0927       1.0983       1.1024       1.1049
  476.    1.1058       1.1051       1.1027       1.0988       1.0934
  477.    1.0865       1.0783       1.0688       1.0580       1.0459 ;
  478.  urr = urr et (prog
  479.    1.0325       1.0177       1.0013       0.98319      0.96315
  480.    0.94083      0.91571      0.88702      0.85347      0.81285
  481.    0.76100      0.68940      0.57826      0.37914      0.11368
  482.    5.00378E-04 );
  483.  krr = prog
  484.   1.14242E-08  5.70695E-04  5.85036E-03  9.18812E-03  9.24352E-03
  485.   8.78015E-03  8.29986E-03  7.86419E-03  7.46515E-03  7.08857E-03
  486.   6.72369E-03  6.36405E-03  6.00640E-03  5.64984E-03  5.29542E-03
  487.   4.94455E-03  4.60029E-03  4.26604E-03  3.94641E-03  3.64662E-03
  488.   3.37300E-03  3.13278E-03  2.93380E-03  2.78421E-03  2.69096E-03
  489.   2.65941E-03  2.69096E-03  2.78421E-03  2.93380E-03  3.13278E-03
  490.   3.37300E-03  3.64662E-03  3.94641E-03  4.26604E-03  4.60029E-03 ;
  491.  krr = krr et (prog
  492.   4.94455E-03  5.29542E-03  5.64984E-03  6.00640E-03  6.36405E-03
  493.   6.72369E-03  7.08857E-03  7.46515E-03  7.86419E-03  8.29986E-03
  494.   8.78015E-03  9.24352E-03  9.18812E-03  5.85036E-03  5.70695E-04
  495.   1.14242E-08);
  496.  mrr = prog
  497.    2.0000       2.0031       2.4276       4.7639       9.0737
  498.    15.067       22.284       30.207       38.403       46.558
  499.    54.454       61.929       68.847       75.089       80.550
  500.    85.148       88.826       91.561       93.372       94.328
  501.    94.549       94.214       93.548       92.808       92.241
  502.    92.030       92.241       92.808       93.548       94.214
  503.    94.549       94.328       93.372       91.561       88.826;
  504.  mrr = mrr et (prog
  505.    85.148       80.550       75.089       68.847       61.929
  506.    54.454       46.558       38.403       30.207       22.284
  507.    15.067       9.0737       4.7639       2.4276       2.0031
  508.    2.0000);
  509.  
  510. evkrr=evol manu xrr krr;
  511. evurr=evol manu xrr urr;
  512. evmrr=evol manu xrr mrr;
  513.  
  514. dess evurr TITR ' Profil de vitesse de référence (KL) en sortie';
  515. dess evkrr TITR ' Profil de K de référence (KL) en sortie';
  516. dess evmrr TITR ' Profil de Mueff/mu de référence (KL) en sortie';
  517.  
  518. *--- Profil de Référence pour le KLbr ----------------------------------
  519.  
  520.  tabk.'Rey'=Rey;
  521.  un = rv.inco.'UN' ;
  522.  
  523.  qs=dbit un $sortie IMPR;
  524.  um=(1./H)*qs;
  525.  unn=(1./um)*un;
  526.  
  527. TAB1 = TABLE ;
  528. tab1.1 = 'MOT' '    ';
  529. tab1.2 = 'MOT' 'MARQ CROI REGU';
  530. tab1.3 = 'MOT' 'MARQ CARR REGU';
  531. tab1.4 = 'MOT' 'MARQ LOSA REGU';
  532. tab1.5 = 'MOT' 'MARQ PLUS REGU';
  533. tab1.6 = 'MOT' 'MARQ TRIA REGU';
  534. tab1.'TITRE'    = TABLE         ;
  535.  tab1.'TITRE'. 1 = 'MOT'  'Re 2.e3' ;
  536.  tab1.'TITRE'. 2 = 'MOT'  'Re 1.e4' ;
  537.  tab1.'TITRE'. 3 = 'MOT'  'Re 1.e5' ;
  538.  tab1.'TITRE'. 4 = 'MOT'  'Re 1.e6' ;
  539.  tab1.'TITRE'. 5 = 'MOT'  'Re 1.e7' ;
  540.  tab1.'TITRE'. 6 = 'MOT' 'Reichardt' ;
  541.  
  542. tit1=chai 'Comparaison Profils de vitesse Buleev à différents Reynolds';
  543. tit2=chai tit1 ' et Reichardt';
  544. dess (evupd1 et evupd2 et evupd3 et evupd4 et evupd5 et evr1)
  545.  LOGX XBOR 0.1 5.e5 'GRIL' 'CARR' lege TAB1 titr tit2;
  546.  
  547. dess (evuref1 et evuref2 et evuref3 et evuref4 et evuref5)
  548.  'GRIL' 'CARR' lege TAB1 titr tit1;
  549. *-----------------------------------------------------------------------
  550.  
  551. unx =exco un 'UX';
  552. dudy=exco (kops 'GRADS' unx $mt) 'UY';
  553. uet=(mu*(abs (dudy + 1.e-5)))**0.5;
  554. uet1=redu uet par1;
  555. uet2=redu uet par2;
  556. evuet1= evol chpo uet1 par1;
  557. evuet2= evol chpo uet2 (inve par2);
  558.  
  559. dess (evuet1 et evuet2) TITR ' U* le long des deux parois'
  560. tity 'm/s' titx 'm';
  561.  
  562. mess ' Liste de U*';
  563. yevuet1=extr evuet1 'ORDO';
  564. yevuet2=extr evuet2 'ORDO';
  565. dl=(dime yevuet1) ;
  566. ll=lect (dl - 10) pas 1 dl;
  567. list (extr yevuet1 ll) ;
  568. list (extr yevuet2 ll) ;
  569.  
  570.  
  571.  f1=8.*(((1./um)*uet1)**2.);
  572.  f2=8.*(((1./um)*uet2)**2.);
  573.  fiso = 0.202*(Rey**(-0.196));
  574.  evf1=evol chpo f1 par1;
  575.  evf1=evf1 coul rouge;
  576.  evf2=evol chpo f2 (inve par2);
  577.  evf2=evf2 coul verte;
  578.  xx=extr evf1 'ABSC';
  579.  yfiso= prog (dime xx)*fiso;
  580.  evfiso=evol manu xx yfiso ;
  581.  evfiso=evfiso coul bleue ;
  582.  
  583. TAB1 = TABLE ;
  584. tab1.1 = 'MOT' 'TIRR';
  585. tab1.2 = 'MOT' 'TIRC';
  586. tab1.3 = 'MOT' 'TIRM';
  587. tab1.'TITRE'    = TABLE         ;
  588.  tab1.'TITRE'. 1 = 'MOT' 'Corr' ;
  589.  tab1.'TITRE'. 2 = 'MOT' 'f paroi1' ;
  590.  tab1.'TITRE'. 3 = 'MOT' 'f paroi2' ;
  591.  
  592.  fmax = maxi yfiso;
  593.  
  594.  dess (evfiso et evf2 et evf1) YBOR 0. (5.*fmax)
  595.  TITR ' Coefficient de frottement' 'GRIL' lege TAB1;
  596.  
  597. tabk.'evfiso'=evfiso;
  598. tabk.'evf2'  =evf2  ;
  599. tabk.'evf1'  =evf1  ;
  600. tabk.'fmax'  =fmax  ;
  601.  
  602. yp1=norm ((poin 2 entre1) moins (poin 1 entre1));
  603. yp2=norm ((poin entre1 final) moins (poin ((nbno entre1) - 1) entre1));
  604.  
  605. mess ' YP1=' YP1 '  YP2=' YP2;
  606.  
  607.  yplus1= yp1*uet1*(rv.inco.'MU'**(-1.));
  608.  yplus2= yp2*uet2*(rv.inco.'MU'**(-1.));
  609.  
  610.  evyplus1=evol chpo yplus1 par1;
  611.  evyplus2=evol chpo yplus2 (inve par2);
  612.  
  613.  dess (evyplus1 et evyplus2) TITR ' Y+ le long des deux parois'
  614. tity ' ' titx 'm';
  615.  
  616.  sort2=elem sortie (lect 1 pas 1 ((nbel sortie)/2));
  617.  sort1= syme sort2 'DROIT' (point 1 axe) (axe poin final);
  618.  elim (sort1 et sort2 et sortie) epsi;
  619.  
  620.  uet1=extr (nbno par1) (extr evuet1 'ORDO');
  621.  uet2=extr (nbno par1) (extr evuet2 'ORDO');
  622.  mess ' uet1 et uet2=' uet1 uet2;
  623.  ypls1= evol chpo ((1./mu) * uet1 * (coor 2 sort1)) (inve sort1);
  624.  ypls2= evol chpo ((1./mu) * uet2 * (H - (coor 2 sort2))) (     sort2);
  625.  
  626.  upl1=(evol chpo unx (inve sort1))*(1./uet1);
  627.  upl2=(evol chpo unx sort2)*(1./uet2);
  628.  xypls1= ORDO (extr ypls1 'ORDO');
  629.  yypls1= ORDO (extr upl1  'ORDO');
  630.  
  631.  xypls1=xypls1 enlev 1;
  632.  yypls1=yypls1 enlev 1;
  633.  evupl1=evol manu xypls1 yypls1;
  634.  
  635.  
  636.  xypls2=extr ypls2 'ORDO';
  637.  yypls2=extr upl2 'ORDO';
  638.  xypls2=xypls2 enlev 1;
  639.  yypls2=yypls2 enlev 1;
  640.  evupl2=evol manu xypls2 yypls2 ;
  641.  evupl1=evupl1 coul rouge;
  642.  evupl2=evupl2 coul verte;
  643.  
  644. TAB1 = TABLE ;
  645. tab1.1 = 'MOT' '    ';
  646. tab1.2 = 'MOT' 'TIRC';
  647. tab1.3 = 'MOT' 'TIRM';
  648. tab1.'TITRE'    = TABLE         ;
  649.  tab1.'TITRE'. 1 = 'MOT' 'Reichardt' ;
  650.  tab1.'TITRE'. 2 = 'MOT' 'U+ paroi1' ;
  651.  tab1.'TITRE'. 3 = 'MOT' 'U+ paroi2' ;
  652.  
  653.  dess (evr1 et evupl2 et evupl1) 'LOGX' 'XBOR' 1. 500.
  654.  TITR ' U+ en sortie' 'GRIL' lege TAB1
  655. tity ' ' titx 'm';
  656.  
  657. tabk.'evr1'  =evr1  ;
  658. tabk.'evupl1'=evupl1;
  659. tabk.'evupl2'=evupl2;
  660.  
  661.  nutsnu=(rv.inco.'MUF')*(1./MU);
  662.  trace nutsnu mt (cont mt) 'TITR' ' Nu effectif / nu ';
  663.  
  664.  npl1=(evol chpo nutsnu (inve sort1));
  665.  npl2=(evol chpo nutsnu sort2);
  666.  npl1=ORDO (extr npl1 'ORDO');
  667.  npl2=ORDO (extr npl2 'ORDO');
  668.  npl1=npl1 enlev 1 ;
  669.  npl2=npl2 enlev 1 ;
  670.  
  671.  evnpl1=evol manu xypls1 npl1;
  672.  evnpl2=evol manu xypls2 npl2;
  673.  evnpl1=evnpl1 coul rouge;
  674.  evnpl2=evnpl2 coul verte;
  675.  
  676. TAB1 = TABLE ;
  677. tab1.1 = 'MOT' 'TIRC';
  678. tab1.2 = 'MOT' 'TIRM';
  679. tab1.'TITRE'    = TABLE         ;
  680.  tab1.'TITRE'. 1 = 'MOT' 'nt/n paroi1' ;
  681.  tab1.'TITRE'. 2 = 'MOT' 'nt/n paroi2' ;
  682.  
  683.  dess (evnpl2 et evnpl1) 'LOGX' 'XBOR' 1. 500.
  684.  TITR ' Nut/nu en sortie' 'GRIL' lege TAB1
  685. tity ' ' titx 'm';
  686.  
  687. tabk.'evnpl1'=evnpl1;
  688. tabk.'evnpl2'=evnpl2;
  689.  
  690.  ung = vect un 0.5 'UX' 'UY' jaune;
  691.  trace ung mt 'TITR' ' Vitesses ';
  692.  
  693.  pn = elno $mt rv.inco.PRESSION kpres;
  694.  trace pn mt (cont mt) TITR 'Pression';
  695.  evpaxe= evol chpo pn axe;
  696.  dess evpaxe TITR ' Pression sur l axe';
  697.  
  698. *-----------------------------------------------------------------------
  699.  evax= evol chpo (exco un 'UX') axe   ;
  700.  dess evax TITR 'U sur l axe'
  701. tity 'm/s' titx 'm';
  702. *-----------------------------------------------------------------------
  703.  
  704. *-----------------------------------------------------------------------
  705.  evmuax= evol chpo nutsnu axe   ;
  706.  dess evmuax TITR 'Nut/Mu sur l axe'
  707. tity ' ' titx 'm';
  708. tabk.'evmuax'=evmuax;
  709. *-----------------------------------------------------------------------
  710.  
  711. *$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  712.  
  713. *profils de référence
  714.  evupdr evupgr evusr uetr = anabu 1. Rey 1. ;
  715. mess '***************************************************************';
  716. mess '***************************************************************';
  717. mess ' UET1 = ' UET1  '  UET2 = ' UET2 '  UETR = ' UETR ;
  718. mess '***************************************************************';
  719. mess '***************************************************************';
  720.  
  721.  evus= evol chpo (exco unn 'UX') sortie;
  722.  usix  = extr evusr 'ABSC';
  723.  usi   = ipol usix (extr evus 'ABSC') (extr evus 'ORDO');
  724.  evusi = evol manu usix usi ;
  725.  evdu  = (evusi - evusr)*(evusi - evusr);
  726.  
  727. tab1 = table;
  728. tab1.1 = 'MOT' 'TIRC';
  729. tab1.1 = 'MOT' '    ';
  730. tab1.2 = 'MOT' '    ';
  731. tab1.3 = 'MOT' '    ';
  732. tab1.'TITRE'    = TABLE         ;
  733.  tab1.'TITRE'. 1 = 'MOT' 'KL Référence';
  734.  tab1.'TITRE'. 2 = 'MOT' Mdl ;
  735.  tab1.'TITRE'. 3 = 'MOT' 'Buleev analytique' ;
  736.  evusr=evusr coul rouge;
  737.  evurr=evurr coul vert ;
  738.  
  739.  dess (evurr et evus et evusr) TITR 'U en sortie' 'GRIL' LEGE TAB1
  740. tity 'm/s' titx 'm';
  741.  
  742.  evkrr=evkrr coul verte;
  743.  evks = evol chpo rv.inco.'KN' sortie;
  744.  
  745.  dess (evkrr et evks) TITR 'K en sortie' 'GRIL' LEGE TAB1
  746. tity 'm2/s2' titx 'm';
  747.  
  748.  evmrr=evmrr coul verte;
  749.  evms = evol chpo ((1./MU)*(rv.inco.'MUF')) sortie;
  750.  
  751.  dess (evmrr et evms) TITR 'Mueff en sortie' 'GRIL' LEGE TAB1
  752. tity 'm2/s' titx 'm';
  753.  
  754.  tabk.'evus'=evus;
  755.  tabk.'evusr'=evusr;
  756.  
  757.  dess evdu TITR '|U - uref|**2';
  758.  
  759.  du    = (INTG evdu);
  760.  errn = (abs (du - 4.06931E-03));
  761.  mess ' delta u = ' du ' errn = ' errn;
  762.  Si (NON COMPLET);
  763.  Si(errn > 0.0001); erreur 5 ; Finsi ;
  764.  Finsi;
  765.  
  766. FINPROC tabk;
  767.  
  768. ******* Profils de T+ de Kader ***********************************
  769. Si FAUX;
  770.  
  771.  evt2p1 evt1p1=kader 0.7 1.e4 'VERTE';
  772.  evt2p2 evt1p2=kader 0.7 1.e5 'ROUGE';
  773.  evt2p3 evt1p3=kader 0.7 1.e6 'BLEUE';
  774.  
  775. tab1 = table;
  776. tab1.1 = 'MOT' '    ';
  777. tab1.2 = 'MOT' 'TIRC';
  778. tab1.3 = 'MOT' '    ';
  779. tab1.4 = 'MOT' 'TIRC';
  780. tab1.5 = 'MOT' '    ';
  781. tab1.6 = 'MOT' 'TIRC';
  782. tab1.'TITRE'    = TABLE         ;
  783.  tab1.'TITRE'. 1 = 'MOT' 'T+1 Ret=1.e4' ;
  784.  tab1.'TITRE'. 2 = 'MOT' 'T+2 Ret=1.e4' ;
  785.  tab1.'TITRE'. 3 = 'MOT' 'T+1 Ret=1.e5' ;
  786.  tab1.'TITRE'. 4 = 'MOT' 'T+2 Ret=1.e5' ;
  787.  tab1.'TITRE'. 5 = 'MOT' 'T+1 Ret=1.e6' ;
  788.  tab1.'TITRE'. 6 = 'MOT' 'T+2 Ret=1.e6' ;
  789.  
  790.  dess (evt1p1 et evt2p1 et evt1p2 et evt2p2 et evt1p3 et evt2p3)
  791.   LOGX XBOR 0.1 1.E5 tity 'T+' titx 'Y+' 'GRIL' LEGE TAB1
  792.  TITR ' Teta+ Cas 1 et 2 Fonction du Reynolds Turbulent';
  793.  
  794.  evt2p1 evt1p1=kader 0.1 1.e4 'VERTE';
  795.  evt2p2 evt1p2=kader 1.  1.e5 'ROUGE';
  796.  evt2p3 evt1p3=kader 10. 1.e6 'BLEUE';
  797.  
  798.  tab1.'TITRE'. 1 = 'MOT' 'T+1 Pr=0.1' ;
  799.  tab1.'TITRE'. 2 = 'MOT' 'T+2 Pr=0.1' ;
  800.  tab1.'TITRE'. 3 = 'MOT' 'T+1 Pr=1.' ;
  801.  tab1.'TITRE'. 4 = 'MOT' 'T+2 Pr=1.' ;
  802.  tab1.'TITRE'. 5 = 'MOT' 'T+1 Pr=10.' ;
  803.  tab1.'TITRE'. 6 = 'MOT' 'T+2 Pr=10.' ;
  804.  
  805.  dess (evt1p1 et evt2p1 et evt1p2 et evt2p2 et evt1p3 et evt2p3)
  806.   LOGX XBOR 0.1 1.E5 tity 'T+' titx 'Y+' 'GRIL' LEGE TAB1
  807.  TITR ' Teta+ Cas 1 et 2 Fonction du Prandtl';
  808. Finsi;
  809. ******* Profils de T+ de Kader  FIN ******************************
  810.  
  811.  
  812. *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  813. *profils de référence
  814.  evupdr evupgr evusr uetr = anabu 1. Rey 1. ;
  815. mess '***************************************************************';
  816. mess '***************************************************************';
  817. mess ' UET1 = ' UET1  '  UET2 = ' UET2 '  UETR = ' UETR ;
  818. mess '***************************************************************';
  819. mess '***************************************************************';
  820.  
  821.  evus= evol chpo (exco unn 'UX') sortie;
  822.  usix  = extr evusr 'ABSC';
  823.  usi   = ipol usix (extr evus 'ABSC') (extr evus 'ORDO');
  824.  evusi = evol manu usix usi ;
  825.  evdu  = (evusi - evusr)*(evusi - evusr);
  826.  
  827.  Si GRAPH ;
  828.   t1=POST;
  829.  
  830. tab1 = table;
  831. tab1.1 = 'MOT' 'TIRC';
  832. tab1.2 = 'MOT' '    ';
  833. tab1.'TITRE'    = TABLE         ;
  834.  tab1.'TITRE'. 1 = 'MOT' 'modèle KL' ;
  835.  tab1.'TITRE'. 2 = 'MOT' 'Buleev analytique' ;
  836.  evusr=evusr coul rouge;
  837.  
  838.  dess (evus et evusr) TITR 'U en sortie' 'GRIL' LEGE TAB1
  839. tity 'm/s' titx 'm';
  840.  
  841.  dess evdu TITR '|U - uref|**2';
  842.  Finsi;
  843.  
  844.  du    = (INTG evdu);
  845.  errn = (abs (du - 2.597E-4));
  846.  mess ' delta u = ' du ' errn = ' errn;
  847.  Si (NON COMPLET);
  848.  Si(errn > 0.0001); erreur 5 ; Finsi ;
  849.  Finsi;
  850.  
  851.  FIN;
  852.  
  853.  
  854.  
  855.  
  856.  
  857.  
  858.  
  859.  
  860.  
  861.  
  862.  
  863.  
  864.  
  865.  
  866.  

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