Télécharger la plaquette Cast3M


fr en

Télécharger recircpp.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  recircpp.dgibi
  2. * tests ajustes par PV le 12/10/08 sur les valeurs obtenues en linux32 
  3. opti dime 2 elem QUA8 ;
  4.  
  5. GRAPH=FAUX ;
  6. *GRAPH=VRAI ;
  7. nbit = 5 ;
  8.  
  9. *
  10. * 1/  Exemple test : canal courbe
  11. *                    Eclt laminaire incompressible
  12. *                    On impose une vitesse a un bout
  13. *                    V normale nulle sur les parois
  14. * On indique comment imposer une Vitesse en entree d'un domaine
  15. * (Vitesse normale non nulle et Vitesse tangentielle nulle !!)
  16. * et comment imposer un glissement sur une paroi (Vitesse normale nulle)
  17. *     Post  : debit sur le contour , Vitesses ,pression
  18. *     Teste : vitesse normale imposee en 2D plan
  19. *             algorithme de projection transitoire
  20. *
  21.  
  22. DISCR='QUAF' ;
  23. KPRES='MSOMMET' ;
  24. R=4. ; n1=5 ; n2=20 ;
  25. *R=4. ; n1=2 ; n2=5  ;
  26. p1=0 0 ; p2= 1 0 ; pc= R 0 ;
  27. entree = p1 d n1 p2 ;
  28. sortie = entree tour pc (-60.) ;
  29. sortie = inve sortie;
  30. q2 = poin 1 sortie ;
  31. q1 = sortie poin final ;
  32. pp1=p1 c n2 pc q1 ;
  33. pp2=p2 c n2 pc q2 ;
  34. paroi= pp1 et pp2 ;
  35. paroi = inve paroi ;
  36.  
  37. cnt= entree et pp2 et sortie et (inve pp1) ;
  38. mt= surf cnt;
  39. Mmt= chan quaf mt ;
  40. Mcnt= chan quaf cnt ;
  41. Mentree= chan quaf entree;
  42. Msortie= chan quaf sortie;
  43. Mparoi = chan quaf paroi ;
  44.  
  45. elim (Mmt et Mcnt et Mentree et Msortie et Mparoi) 1.e-5 ;
  46. $mt = mode Mmt 'NAVIER_STOKES' DISCR ;
  47. mt= doma $mt maillage ;
  48. $entree= mode Mentree 'NAVIER_STOKES' DISCR ;
  49. $sortie= mode Msortie 'NAVIER_STOKES' DISCR ;
  50. $paroi = mode Mparoi  'NAVIER_STOKES' DISCR ;
  51. entrep=doma $entree 'MSOMMET' ;
  52. sortie=doma $sortie maillage ;
  53.  
  54. NU= 1.e-2 ;
  55. U0=1. ;
  56. DT=2. ;
  57.  
  58. * Comment imposer un debit d'entree
  59. * En general il faut imposer la composante normale a
  60. * la vitesse debitante et les composantes tangentielles a ZERO !
  61. * On n'utilise donc pas VNIMP :
  62.  
  63. * On extrait le champ des normales de la surface d'entree
  64. * en prenant garde a l'orientation des normales
  65. nsor=doma $sortie 'NORMALEV';
  66. * On extrait les composantes nx ny etc
  67. nx= exco 'UX' nsor ;
  68. ny= exco 'UY' nsor ;
  69. * on cree les componsantes cartesiennes du champ de vitesse
  70. * qui seront imposees en condition limite type Dirichlet
  71. usx=u0*nx ;
  72. usy=u0*ny ;
  73.  
  74. * Sur les parois la condition de glissement s'obtient
  75. * effectivement en imposant une vitesse normale nulle a l'aide de VNIMP
  76.  
  77.  
  78. rv = eqex 'ITMA' 1  'NITER' 1  'OMEGA' 1.
  79.   OPTI EF 'CENTREE' 'IMPL' KPRES
  80.   'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  81.   'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  82.   'CLIM'
  83.   'UN' UIMP sortie  usx  'UN' VIMP sortie usy
  84.   ;
  85.  
  86. rv.'INCO'= TABLE 'INCO' ;
  87. rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
  88. rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
  89.  
  90.     RVP= EQEX 'OPTI' 'EF' KPRES
  91.    'ZONE' $mt 'OPER' 'KBBT'  -1.  'INCO' 'UN' 'PRES'
  92.    'ZONE' $paroi 'OPER' 'VNIMP' $mt 0. 'INCO' 'UN' 'PRES'
  93.    'CLIM'
  94.    'PRES' 'TIMP' entrep 0. ;
  95.     ;
  96.  
  97. *  rvp. 'METHINV' . 'FCPRECT'=100 ;
  98.  
  99.     RV.'PROJ'= RVP ;
  100.  
  101.  
  102.  exec rv ;
  103.  
  104. * affichage vitesse
  105. un=rv.inco.'UN' ;
  106. ung=vect un 0.1 ux uy jaune ;
  107.  
  108. Si GRAPH ;trace ung mt;Finsi ;
  109. cnt1=chan cnt seg2 ;                                                    
  110. Mcnt1= chan cnt1 quaf;                                                  
  111. $cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
  112. q=dbit (rv.inco.'UN') $cnt ;                                            
  113. Mess ' Debit sur le contour ' q ;
  114.  dq= (abs q) - 3.36856E-02 ;
  115.  list dq ;
  116.  si (dq  > 1.e-5 ) ; erreur 5 ; finsi ;
  117.  
  118. * affichage pression
  119. * On extrait la bonne composante de la pression
  120. * il y a aussi les multiplicateurs de Lagrange vitesse normale
  121. pn=exco rv.inco.'PRESSION' 'PRES' ;
  122.  
  123. * On reaffecte les bon SPG au CHPOINT
  124. pn= kcht $mt scal KPRES pn ;
  125. mtp= doma $mt 'MMAIL' ;
  126.  
  127. Si GRAPH ;trace pn mtp;Finsi ;
  128. mpn=(maxi pn);
  129. dpn= abs (mpn - 1.7705);
  130. mess ' MAXI PN = ' mpn 'dpn=' dpn ;
  131.  si ( dpn > 1.e-4 ) ; erreur 5 ; finsi ;
  132.  
  133.  
  134.  
  135. *
  136. *  2/ Exemple test : recirculation dans une cavite semi circulaire
  137. *                    Eclt laminaire incompressible
  138. *                    On impose une vitesse tg sur la surface libre
  139. *                    V normale nulle sur le fond
  140. *     Voir RECIRC.DGIBI pour l'algorithme Taylor Hood
  141. *     Post  : debit sur le contour , fonction de courant , pression
  142. *     Teste : vitesse normale imposee en 2D plan
  143. *             algorithme de projection transitoire
  144. *
  145.  
  146. p0=0. 0. ;
  147. p1= -1. 0.;
  148. p2= 1. 0. ;
  149. p03= 0. -1. ;
  150. p3= 1. -1. ;
  151. p4= -1. -1. ;
  152.  
  153. n1=  5; n2= 5;
  154. ab= p1 d n1 p2 ;
  155. slib=  (chan ab POI1) ;
  156. nab= nbel slib ;
  157. slib= elem slib (lect 2 pas 1 (nab - 1));
  158.  
  159. bc= p2 c p0 n2 p03 c  p0 n2 p1;
  160. *bc= p2 d n2 p3 d n2 p4 d n2 p1 ;
  161.  
  162. cnt=ab et bc ;
  163. opti elem TRI6 ;
  164. mt= cnt surf ;
  165. Mmt=chan mt quaf ;
  166. Mbc=chan bc quaf ;
  167. elim (Mmt et Mbc) 1.e-5 ;
  168.  
  169. $mt = mode Mmt 'NAVIER_STOKES' DISCR;
  170. $bc = mode Mbc 'NAVIER_STOKES' DISCR;
  171.  
  172. * La cavité est fermée il faut imposer la pression en 1 point !
  173. prep1=doma $mt KPRES ;
  174. bcp=elem  prep1 POI1 (lect 1) ;
  175. *
  176.  
  177. *NU= 1. / 61. ;
  178. NU= 1.e-2 ;
  179. U0=1. ;
  180. DT=1. ;
  181.  
  182. rv = eqex 'ITMA' 10 'NITER' 1  'OMEGA' 1.
  183.   OPTI EF 'CENTREE' 'IMPL' KPRES
  184.   'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  185.   'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  186.   'CLIM'
  187.          'UN' UIMP slib 1. 'UN' VIMP ab 0.
  188.   ;
  189.  
  190. rv.'INCO'= TABLE 'INCO' ;
  191. rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
  192. rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
  193. rv.'INCO'.'VN'=kcht $bc scal sommet 0. ;
  194.  
  195.     RVP= EQEX 'OPTI' 'EF' KPRES
  196.    'ZONE' $mt 'OPER' 'KBBT'  -1.  'INCO' 'UN' 'PRES'
  197.    'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
  198.     CLIM PRES TIMP bcp 0. ;
  199.  
  200. *  rvp. 'METHINV' . 'FCPRECT'=100 ;
  201.  
  202.     RV.'PROJ'= RVP ;
  203.  
  204.  
  205. exec rv ;
  206.  
  207. un = rv.inco.'UN' ;
  208. ung= vect un 0.11  ux uy jaune ;
  209.  
  210. Si GRAPH ;trace ung Mmt ; Finsi ;
  211.  
  212. cnt1=chan cnt seg2 ;                                                    
  213. Mcnt1= chan cnt1 quaf;                                                  
  214. $cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
  215. q=dbit (rv.inco.'UN') $cnt ;                                            
  216. Mess ' Debit sur le contour ' q ;
  217.  dq= (q - .22362) abs; list dq;
  218.  si ( dq > 1.e-5 ) ; erreur 5 ; finsi ;
  219.  
  220. *
  221. * Calcul de la fonction de courant
  222. *
  223.  
  224. sw = 'KOPS' un 'ROT' $mt ;
  225. rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
  226.      'ZONE' $mt 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI'
  227.      'ZONE' $mt 'OPER' 'FIMP' sw      'INCO' 'PSI'
  228.      'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
  229. rk . 'INCO' = 'TABLE' 'INCO' ;
  230. rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;
  231.  
  232. opti isov SULI ;
  233. exec rk ;
  234. psi=rk.inco.'PSI' ;
  235. Si GRAPH ;trace psi Mmt cnt ; Finsi ;
  236. mpsi=(mini psi);
  237. dpsi= abs (mpsi + .17603) ;
  238. mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
  239.  si ( dpsi > 5.e-4 ) ; erreur 5 ; finsi ;
  240.  
  241. *
  242. * Calcul de la pression
  243. *
  244.  
  245.   pn=kcht $mt scal KPRES (exco rv.inco.'PRESSION' 'PRES') ;
  246.   mtp=doma $mt 'MMAIL' ;
  247. Si GRAPH ;  trace pn mtp ; Finsi ;
  248. mpn=(maxi pn);
  249. dpn= abs (mpn - .38065);
  250. mess ' MAXI PN = ' mpn 'dpn=' dpn ;
  251.  si ( dpn > 1.e-3 ) ; erreur 5 ; finsi ;
  252.  
  253.  
  254. *
  255. *  Troisieme cas : Meme probleme geometrie rectangulaire
  256. *                  MACRO CENTRE
  257. *
  258. *
  259. *
  260. *
  261.  
  262. DISCR = 'QUAF';
  263. KPRES= 'MSOMMET' ;
  264.  
  265. p0=0. 0. ;
  266. p1= -1. 0.;
  267. p2= 1. 0. ;
  268. p03= 0. -1. ;
  269. p3= 1. -1. ;
  270. p4= -1. -1. ;
  271.  
  272. n1=  5 ; n2= 5 ;
  273. ab= p1 d n1 p2 ;
  274. slib=  (chan ab POI1) ;
  275. nab= nbel slib ;
  276. slib= elem slib (lect 2 pas 1 (nab - 1));
  277.  
  278. bc= p2 d n2 p3 d n2 p4 d n2 p1 ;
  279.  
  280. cnt=ab et bc ;
  281.  
  282. mt= cnt surf ;
  283. Mmt=chan mt quaf ;
  284. Mbc=chan bc quaf ;
  285. elim (Mmt et Mbc) 1.e-5 ;
  286.  
  287. $mt = mode Mmt 'NAVIER_STOKES' DISCR;
  288. $bc = mode Mbc 'NAVIER_STOKES' DISCR;
  289.  
  290. * La cavité est fermée il faut imposer la pression en 1 point !
  291. prep1=doma $mt KPRES ;
  292. bcp=elem  prep1 POI1 (lect 1) ;
  293. *
  294.  
  295. NU= 1.e-2 ;
  296. U0=1. ;
  297. DT=1. ;
  298.  
  299. rv = eqex 'ITMA' 10 'NITER' 1  'OMEGA' 1.
  300.   OPTI EF 'CENTREE' 'IMPL' KPRES
  301.   'ZONE' $MT 'OPER' 'NS' 1. 'UN'  NU  'INCO' 'UN'
  302.   'ZONE' $MT 'OPER' 'DFDT' 1. 'UN' DT 'INCO' 'UN'
  303.   'CLIM'
  304.          'UN' UIMP slib 1. 'UN' VIMP ab 0.
  305.   ;
  306.  
  307. rv.'INCO'= TABLE 'INCO' ;
  308. rv.'INCO'.'UN'=kcht $mt vect sommet (0. 0.) ;
  309. rv.'INCO'.'PRES'=kcht $mt scal KPRES 0. ;
  310. rv.'INCO'.'VN'=0. ;
  311.  
  312.     RVP= EQEX 'OPTI' 'EF' KPRES
  313.    'ZONE' $mt 'OPER' 'KBBT'  -1.  'INCO' 'UN' 'PRES'
  314.    'ZONE' $bc 'OPER' 'VNIMP' $mt 'VN' 'INCO' 'UN' 'PRES'
  315.     CLIM PRES TIMP bcp 0. ;
  316.     RV.'PROJ'= RVP ;
  317.  
  318. exec rv ;
  319.  
  320.  
  321. un = rv.inco.'UN' ;
  322. ung= vect un 0.11  ux uy jaune ;
  323.  
  324. Si GRAPH ;trace ung Mmt ; Finsi ;
  325.  
  326. cnt1=chan cnt seg2 ;                                                    
  327. Mcnt1= chan cnt1 quaf;                                                  
  328. $cnt=model Mcnt1 'NAVIER_STOKES' line ;                                 
  329. q=dbit (rv.inco.'UN') $cnt ;                                            
  330. Mess ' Debit sur le contour ' q ;
  331.  dq=abs ( q - 2.78540E-03);
  332.  si ( dq  > 1.e-3 ) ; erreur 5 ; finsi ;
  333.  
  334. sw = 'KOPS' un 'ROT' $mt ;
  335. rk = 'EQEX' 'OPTI' 'EF' 'IMPL'
  336.      'ZONE' $mt 'OPER' 'LAPN' 1.0D0   'INCO' 'PSI'
  337.      'ZONE' $mt 'OPER' 'FIMP' sw      'INCO' 'PSI'
  338.      'CLIM' 'PSI' 'TIMP' cnt 0.0D0 ;
  339. rk . 'INCO' = 'TABLE' 'INCO' ;
  340. rk . 'INCO' . 'PSI' = 'KCHT' $mt 'SCAL' 'SOMMET' 0.0D0 ;
  341.  
  342. opti isov SULI ;
  343. exec rk ;
  344. psi=rk.inco.'PSI' ;
  345. mess ' MINI PSI = ' (mini psi) ;
  346. Si GRAPH ;trace psi Mmt cnt ; Finsi ;
  347. mpsi=(mini psi);
  348. dpsi= abs (mpsi + .18561) ;
  349. mess ' MINI PSI = ' mpsi 'dpsi=' dpsi ;
  350.  si ( dpsi > 1.e-3 ) ; erreur 5 ; finsi ;
  351.  
  352.  
  353. *
  354. * Calcul de la pression
  355. *
  356.   pn=kcht $mt scal KPRES (exco rv.inco.'PRESSION' 'PRES') ;
  357.   mtp=doma $mt 'MMAIL';
  358. Si GRAPH ;  trace pn mtp ; Finsi ;
  359. mpn=(maxi pn);
  360. dpn= abs (mpn - .43612);
  361. mess ' MAXI PN = ' mpn 'dpn=' dpn ;
  362.  si ( dpn > 0.1   ) ; erreur 5 ; finsi ;
  363.  
  364. FIN ;
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  

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