Télécharger la plaquette Cast3M


fr en

Télécharger recircp.dgibi

Retour à la liste

Numérotation des lignes : 

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

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