Télécharger la plaquette Cast3M


fr en

Télécharger test1fpu.dgibi

Retour à la liste

Numérotation des lignes : 

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

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