Télécharger la plaquette Cast3M


fr en

Télécharger canalKL.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : canalKL.dgibi
  2. * Ce cas teste le modèle K-L bas Reynolds (Wolfshtein Yap)
  3. * sur l'écoulement turbulent dans un canal plan.
  4. *  opti trace 'PSC';
  5. mess ' Configuration Courte    ' ; COMPLET = FAUX ;
  6. mess ' Configuration SANS TRACE' ; GRAPH   = FAUX ;
  7. *mess ' Configuration Complete  ' ; COMPLET = VRAI ;
  8. *mess ' Configuration POST      ' ; GRAPH   = VRAI ;
  9.  
  10. *-----------------------------------------------------------------------
  11.  opti dime 2 elem TRI6;
  12.  opti dime 2 elem QUA8;
  13. *-----------------------------------------------------------------------
  14. *      |Triangles | Quadrangles
  15. * ---------------------------
  16. * LINE |          | OK
  17. * LINB |          |
  18. * MACRO|   ?      | OK
  19. * QUAF |          | OK
  20. *-----------------------------------------------------------------------
  21.  
  22. *-----------------------------------------------------------------------
  23. *--- Relecture d'un profil de vitesse asymptotique de Référence --------
  24.  Si (NON COMPLET);
  25.  zsr=prog 0.
  26.   2.22209E-03  9.51930E-03  2.15866E-02  3.81009E-02  5.87392E-02
  27.   8.31785E-02  0.11110      0.14217      0.17607      0.21248
  28.   0.25108      0.29155      0.33355      0.37676      0.42088
  29.   0.46556      0.51049      0.55534      0.59980      0.64353
  30.   0.68623      0.72755      0.76718      0.80480      0.84008
  31.   0.87270      0.90234      0.92868      0.95138      0.97014
  32.   0.98462      0.99450      0.99946      1.       ;
  33.  usr=prog 0.
  34.   0.25980      0.60654      0.74057      0.81674      0.87078
  35.   0.91336      0.94860      0.97849       1.0042       1.0263
  36.   1.0454       1.0617       1.0752       1.0858       1.0935
  37.   1.0979       1.0988       1.0962       1.0903       1.0812
  38.   1.0692       1.0544       1.0369       1.0164       0.99266
  39.   0.96512      0.93291      0.89458      0.84739      0.78511
  40.   0.68969      0.48982      6.49844E-02  0.       ;
  41.  Sinon;
  42. option acquerir 'ucanal_data.dgibi';
  43. acque aa*mot;
  44. acque aa*mot;
  45. acque aa*mot;
  46. acque aa*mot;
  47. acque aa*mot;
  48. acque aa*mot;
  49. aa=1.;
  50. acque z*flottant u*flottant v*flottant;
  51. usr=prog u; zsr = prog z;
  52.  
  53.  REPETER BLOC1 501;
  54.  acque z*flottant u*flottant v*flottant;
  55.  usr=usr et (prog u);
  56.  zsr=zsr et (prog z);
  57.  FIN BLOC1 ;
  58.  
  59.  Finsi;
  60.  
  61.  evusr = evol manu zsr usr;
  62.  su= (INTG evusr);
  63.  ep = maxi zsr;
  64.  um=su/ep;
  65. mess ' Vitesse moyenne ' um;
  66. usr = usr * (1./um);
  67. zsr = zsr * (1./ep);
  68. evusr = evol manu zsr usr ;
  69. Si GRAPH ; dess evusr TITR ' Profil de vitesse de référence'; Finsi;
  70. *-----------------------------------------------------------------------
  71.  
  72. *-----------------------------------------------------------------------
  73. *--- Profil de Reichardt -----------------------------------------------
  74. Ka = 0.41;
  75. yplusr= prog 1. pas 0.5 1500.;
  76. unp = prog (dime yplusr)*1.;
  77.  
  78. E=9.5;
  79.  
  80. uplusr=(1./Ka)*(Log (unp + (Ka*yplusr)));
  81. uplusr=uplusr + (7.8* (unp - (exp (yplusr*(-1./11))))
  82.  - ( (yplusr*(-1./11)) * (exp (yplusr*(-1./3.))) ) );
  83.  
  84. evr1=evol manu   yplusr uplusr;
  85. Si GRAPH ;dess evr1 'LOGX' TITR 'Profil de Reichardt'; Finsi;
  86.  
  87. *--- Profil de Reichardt -----------------------------------------------
  88. *-----------------------------------------------------------------------
  89.  
  90.  DISCR=QUAF;
  91. *DISCR=MACRO;
  92. *DISCR=LINE;
  93. KPRES=CENTREP1;
  94. *KPRES=MSOMMET ;
  95. *KPRES=CENTRE  ;
  96. KSUPG=SUPG;
  97. KSUPG=CENTREE;
  98. KSUPG=SUPG   ;
  99. epsi = 1.e-5;
  100.  
  101. ******=================================================**********
  102. ******=================================================**********
  103. ******=================================================**********
  104. ******=================================================**********
  105. ******=================================================**********
  106. Si COMPLET;
  107.  NBIT =400 ;
  108.  DT  = 0.01;
  109.  H=1. ;
  110.  L=40.*H;
  111.  dcl=H/120;
  112.  dm=H/20;
  113.  dcx=H/15.;
  114.  dmx=H/10.;
  115. Sinon;
  116.  NBIT =7;
  117.  DT  = 0.1;
  118.  H=1. ;
  119.  L=10.*H;
  120.  dcl=H/20;
  121.  dm=H/5;
  122.  dcx=H/10.;
  123.  dmx=H/5.;
  124.  
  125. P0 = 0. 0. ;
  126. P1 = L  0. ;
  127. P2 = L H   ;
  128. P3 = 0. H  ;
  129. PS3 = (P0 PLUS P3)*0.5;
  130.  
  131. entree= P0 d dini dcl dfin dm ps3 d dini dm dfin dcl P3;
  132. *entree= P0 d 12 P3;
  133. sortie=(entree plus (L 0.)) inve;
  134. par1  = p0 d dini dcx dfin dmx p1 ;
  135. *par1  = p0 d dini dcl dfin dcl p1 ;
  136. par2  = (par1 plus (0. H)) inve;
  137.  
  138. pc= 1. 1.;
  139. ang= 43.;
  140. ang= 60.;
  141.  ang=  0. ;
  142. entree=entree TOUR ang  pC;
  143. sortie=sortie TOUR  ang pC;
  144. par1  =par1   TOUR  ang pC;
  145. par2  =par2   TOUR  ang pC;
  146.  
  147. elim (entree et par1 et par2 et sortie) epsi;
  148. mt = daller entree par2 sortie par1;
  149. v1=1. 0. ;
  150. v2=0. -1.;
  151. v3=-1. 0.;
  152. v4=0. 1. ;
  153. vt1= tour v1  ang p0;
  154. vt2= tour v2  ang p0;
  155. vt3= tour v3  ang p0;
  156. vt4= tour v4  ang p0;
  157. list vt1;
  158. list vt2;
  159. list vt3;
  160. list vt4;
  161.  
  162. Mmt    = chan mt     QUAF ;
  163. Mentree= chan entree QUAF ;
  164. Msortie= chan sortie QUAF ;
  165. Mpar1  = chan par1   QUAF ;
  166. Mpar2  = chan par2   QUAF ;
  167.  
  168. elim (Mmt et Mentree et Msortie et Mpar1 et Mpar2) epsi ;
  169.  
  170. $mt     = model Mmt     'NAVIER_STOKES' DISCR;
  171. $entree = model Mentree 'NAVIER_STOKES' DISCR;
  172. $sortie = model Msortie 'NAVIER_STOKES' DISCR;
  173. $par1   = model Mpar1   'NAVIER_STOKES' DISCR;
  174. $par2   = model Mpar2   'NAVIER_STOKES' DISCR;
  175.  
  176. doma $mt 'IMPR';
  177. doma $entree 'IMPR';
  178. mt     = doma $mt     maillage;
  179. entree = doma $entree maillage;
  180. sortie = doma $sortie maillage;
  181. par1   = doma $par1   maillage;
  182. par2   = doma $par2   maillage;
  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. lm1=0.41 *(coor 2 mt);
  189. ilm1= lm1 masq 'INFERIEUR' (0.4 *H/2.) ;
  190. lm2=0.41*(H - (coor 2 mt));
  191. ilm2= lm2 masq 'INFERIEUR' (0.4 *H/2.) ;
  192. lma=0.41 *(coor 2 axe);
  193. lm = (lm1*ilm1)+(lm2*ilm2) + lma;
  194. dparoi=(coor 2 mt);
  195. dparoi=dparoi*(H - dparoi);
  196. dparoi=dparoi+1.e-10;
  197. Si GRAPH ;
  198. dess (evol chpo dparoi sortie) TITR ' Distance à la paroi';
  199. trace mt TITR ' Maillage ';
  200. Finsi ;
  201.  
  202. ye= coor 2 entre1;
  203. Rey = 10000.;
  204. U0  = 1. ;
  205. Ufr = U0 ;
  206. *Ufr = ye*(H - ye)*U0*1.5*4. ;
  207. uux=(cos ang) * U0;
  208. uuy=(sin ang) * U0;
  209.  
  210. NU  = 1./Rey;
  211. NU  = 3.14E-5;
  212. YP  = 1.e-2 ;
  213. fnu=1.e10;
  214. fnu=Nu ;
  215.  
  216. RV= EQEX  'OMEGA' 1.    'NITER' 1  ITMA NBIT    'FIDT' 1
  217. 'OPTI' 'EF' 'IMPL'   KSUPG  KPRES
  218.   ZONE  $mt  'OPER' 'KEPSILON' 1. 'UN' 'NU' DT   'INCO' 'KN' 'EN'
  219.   ZONE  $mt  'OPER' 'NS'       1. 'UN' 'MUF'     'INCO' 'UN'
  220. 'OPTI' 'EFM1' 'IMPL' 'CENTREE'
  221.   ZONE  $mt  'OPER'  DFDT 1. 'UNM' DT  INCO 'UN';
  222.   ;
  223.  
  224. rv.'ALGO_KEPSILON'=mots 'KLbr';
  225.  
  226.   RV= EQEX RV  CLIM
  227. * UN UIMP entre1 Ufr  UN VIMP entree 1.e-4
  228.   UN UIMP entre1 uux  UN VIMP entree (uuy + 1.e-4)
  229.   UN UIMP (par1 et par2) 0. UN VIMP (par1 et par2) 0.
  230.   KN TIMP entre1 ((0.02*Ufr)**2.)
  231.   KN TIMP (par1 et par2) 0. EN TIMP (par1 et par2) 0.
  232.   ;
  233.  
  234.  rv.'METHINV'.TYPINV=3 ;
  235.  rv.'METHINV'.IMPINV=0 ;
  236.  rv.'METHINV'.NITMAX=5000;
  237.  rv.'METHINV'.PRECOND=3 ;
  238. *rv.'METHINV'.PRECOND=2 ;
  239.  rv.'METHINV'.ILUTLFIL=4 ;
  240.  rv.'METHINV'.ILUTDTOL=0.;
  241.  rv.'METHINV'.RESID  =1.e-8;
  242.  rv.'METHINV' . 'FCPRECT'=1   ;
  243.  rv.'METHINV' . 'FCPRECI'=1   ;
  244.  
  245.     RVP= EQEX
  246.     'OPTI' 'EF' KPRES
  247.     ZONE  $mt    'OPER'  KBBT  (-1.) 'INCO' 'UN' 'PRES'
  248.     ;
  249.  
  250. *   rvp.'METHINV'.TYPINV=TYPINVPR ;
  251.     rvp.'METHINV'.TYPINV=3        ;
  252.     rvp.'METHINV'.IMPINV=0 ;
  253.     rvp.'METHINV'.NITMAX=2500;
  254.     rvp.'METHINV'.PRECOND=3 ;
  255.     rvp.'METHINV'.ILUTLFIL=10;
  256.     rvp.'METHINV'.ILUTDTOL=0.;
  257.     rvp.'METHINV'.RESID  =1.e-10 ;
  258.     rvp.'METHINV' . 'FCPRECT'=50000 ;
  259.     rvp.'METHINV' . 'FCPRECI'=50000 ;
  260.  
  261.  
  262.     RV.'PROJ'= RVP ;
  263.  
  264. rv.'INCO'= table inco ;
  265. rv.'INCO'.'UN'  =  kcht $mt  vect sommet (1.e-5 1.e-5) ;
  266. rv.'INCO'.'KN'  =  kcht $mt  scal sommet 1.e-5 ;
  267. rv.'INCO'.'Echl'=  kcht $mt  scal sommet dparoi ;
  268. rv.'INCO'.'dparoi'=  kcht $mt  scal sommet dparoi ;
  269. rv.'INCO'.'EN'  =  kcht $mt  scal sommet 1.e-5 ;
  270. rv.'INCO'.'UNM' =  kcht $mt  vect sommet (1.e-5 1.e-5) ;
  271. rv.'INCO'.'UNMM'=  kcht $mt  vect sommet (1.e-5 1.e-5) ;
  272. rv.'INCO'.'PRES'=  kcht $mt  scal kpres  0.;
  273. rv.'INCO'.'NU'=NU;
  274.  
  275.  
  276.  rv.'LTPS' = PROG ;
  277.  rv.'Tps' = 1. ;
  278.  rv.'DT'  = DT;
  279.  rv.'NUPAT' =1;
  280.  
  281. mess ' EXECUTION '; exec rv ;
  282.  
  283.  un = rv.inco.'UN' ;
  284.  
  285. Si GRAPH;
  286.  
  287. unx =exco un 'UX';
  288. dudy=exco (kops 'GRADS' unx $mt) 'UY';
  289. uet=(nu*(abs (dudy + 1.e-5)))**0.5;
  290. uet1=redu uet par1;
  291. uet2=redu uet par2;
  292. evuet1= evol chpo uet1 par1;
  293. evuet2= evol chpo uet2 (inve par2);
  294.  
  295. dess (evuet1 et evuet2) TITR ' U* le long des deux parois';
  296. mess ' Nb de pts ' (nbno entre1);
  297.  
  298. yp1=norm ((poin 2 entre1) moins (poin 1 entre1));
  299. yp2=norm ((poin entre1 final) moins (poin ((nbno entre1) - 1) entre1));
  300.  
  301. mess ' YP1=' YP1 '  YP2=' YP2;
  302.  
  303.  yplus1= yp1*uet1*(1./ rv.inco.'NU');
  304.  yplus2= yp2*uet2*(1./ rv.inco.'NU');
  305.  evyplus1=evol chpo yplus1 par1;
  306.  evyplus2=evol chpo yplus2 (inve par2);
  307.  
  308.  dess (evyplus1 et evyplus2) TITR ' Y+ le long des deux parois';
  309.  
  310.  sort2=elem sortie (lect 1 pas 1 ((nbel sortie)/2));
  311.  sort1= syme sort2 'DROIT' (point 1 axe) (axe poin final);
  312.  elim (sort1 et sort2 et sortie) epsi;
  313.  
  314.  uet1=extr (nbno par1) (extr evuet1 'ORDO');
  315.  uet2=extr 1 (extr evuet2 'ORDO');
  316.  mess ' uet1 et uet2=' uet1 uet2;
  317.  ypls1= evol chpo ((1./nu) * uet1 * (coor 2 sort1)) (inve sort1);
  318.  ypls2= evol chpo ((1./nu) * uet2 * (H - (coor 2 sort2))) (     sort2);
  319.  
  320.  upl1=(evol chpo unx (inve sort1))*(1./uet1);
  321.  upl2=(evol chpo unx sort2)*(1./uet2);
  322.  xypls1= ORDO (extr ypls1 'ORDO');
  323.  yypls1= ORDO (extr upl1  'ORDO');
  324.  
  325.  evupl1=evol manu xypls1 yypls1;
  326.  
  327.  evupl2=evol manu (extr ypls2 'ORDO')(extr upl2 'ORDO');
  328.  evupl1=evupl1 coul rouge;
  329.  evupl2=evupl2 coul verte;
  330.  
  331. TAB1 = TABLE ;
  332. tab1.1 = 'MOT' 'TIRC';
  333. tab1.2 = 'MOT' 'TIRR';
  334. tab1.3 = 'MOT' 'TIRM';
  335. tab1.'TITRE'    = TABLE         ;
  336.  tab1.'TITRE'. 1 = 'MOT' 'Reichardt' ;
  337.  tab1.'TITRE'. 2 = 'MOT' 'U+ paroi1' ;
  338.  tab1.'TITRE'. 3 = 'MOT' 'U+ paroi2' ;
  339.  
  340.  dess (evupl2 et evupl1 et evr1) 'LOGX' 'XBOR' 1. 500.
  341.  TITR ' U+ en sortie' 'GRIL' ;
  342.  
  343.  nutsnu=(rv.inco.'MUF')*(1./NU);
  344.  trace nutsnu mt (cont mt) 'TITR' ' Nu effectif / nu ';
  345.  
  346.  ung = vect un 0.5 'UX' 'UY' jaune;
  347.  trace ung mt 'TITR' ' Vitesses ';
  348.  
  349.  pn = elno $mt rv.inco.PRESSION kpres;
  350.  trace pn mt (cont mt) TITR 'Pression';
  351.  evpaxe= evol chpo pn axe;
  352.  dess evpaxe TITR ' Pression sur l axe';
  353.  
  354.  evax= evol chpo (exco un 'UX') axe   ;
  355.  dess evax TITR 'U sur l axe';
  356.  
  357. Finsi ;
  358.  
  359.  qs=dbit un $sortie IMPR;
  360.  um=(1./H)*qs;
  361.  unn=(1./um)*un;
  362.  
  363.  evus= evol chpo (exco unn 'UX') sortie;
  364.  usix  = extr evusr 'ABSC';
  365.  usi   = ipol usix (extr evus 'ABSC') (extr evus 'ORDO');
  366.  evusi = evol manu usix usi ;
  367.  evdu  = (evusi - evusr)*(evusi - evusr);
  368.  
  369.  Si GRAPH ;
  370.  dess (evus et evusr) TITR 'U en sortie';
  371.  dess evdu TITR '|U - uref|**2';
  372.  Finsi;
  373.  
  374.  du    = (INTG evdu);
  375.  errn = (abs (du - 5.89880E-03));
  376.  mess ' delta u = ' du ' errn = ' errn;
  377.  Si (NON COMPLET);
  378.  Si(errn > 0.0011); erreur 5 ; Finsi ;
  379.  Finsi;
  380.  
  381.  FIN;
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  

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