Télécharger vortex.dgibi

Retour à la liste

Numérotation des lignes :

  1. *****************************************************
  2. * fichier : vortex.dgibi *
  3. ** modifie le 15/06/2014 passage EQPR -> EQEX *
  4. *****************************************************
  5. COMPLET = FAUX ;
  6. SI ( COMPLET ) ;
  7. NX = 20 ;
  8. NITER = 2000 ;
  9. CFL = 1.0 ;
  10. SINON ;
  11. NX = 5 ;
  12. NITER = 250 ;
  13. CFL = 1.0 ;
  14. FINSI ;
  15. GRAPH = 'N' ;
  16. KPRESS = 'CENTRE';
  17. BETA=1. ;
  18. * pour triangles:
  19. * BETA = 1000. ;
  20.  
  21.  
  22. *********************************************************************
  23. * SIMULATION D'UN TOURBILLON - COMPARAISON AVEC SOLUTION ANALYTIQUE *
  24. * DES EQUATIONS DE NAVIER-STOKES *
  25. * OPERATEUR NS - OPTION SUPG *
  26. * H. Paillere - Avril 1997 *
  27. *********************************************************************
  28.  
  29. ******************************************************************
  30. * PROCEDURE POUR ESTIMER LA CONVERGENCE VERS L'ETAT STATIONNAIRE *
  31. ******************************************************************
  32.  
  33. DEBPROC CALCUL ;
  34. ARGU RVX*'TABLE' ;
  35.  
  36. RV = RVX.'EQEX' ;
  37.  
  38. UN = RV.INCO.'UN' ;
  39. UNM1 = RV.INCO.'UNM1' ;
  40.  
  41. DD = RV.PASDETPS.'NUPASDT' ;
  42. NN = DD/5 ;
  43. LO = (DD-(5*NN)) EGA 0 ;
  44. SI ( LO ) ;
  45.  
  46. unx= kcht $MT scal sommet (exco 'UX' un) ;
  47. unm1x = kcht $MT scal sommet (exco 'UX' unm1) ;
  48. uny= kcht $MT scal sommet (exco 'UY' un) ;
  49. unm1y = kcht $MT scal sommet (exco 'UY' unm1) ;
  50.  
  51. ERRX = KOPS unx '-' unm1x ;
  52. ERRY = KOPS uny '-' unm1y ;
  53. ELIX = MAXI ERRX 'ABS' ;
  54. ELIY = MAXI ERRY 'ABS' ;
  55. ELIX = (LOG (ELIX + 1.0E-20))/(LOG 20.) ;
  56. ELIY = (LOG (ELIY + 1.0E-20))/(LOG 20.) ;
  57.  
  58. pn=elno $MT (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  59. PMAX = MAXI PN ;
  60. PMIN = MINI PN ;
  61. DP = PMAX - PMIN ;
  62. MESSAGE 'ITER ' RV.PASDETPS.'NUPASDT' ' ERREUR LINF ' ELIY
  63. 'PMAX = ' PMAX 'PMIN = ' PMIN 'DP = ' DP ;
  64. IT = PROG RV.PASDETPS.'NUPASDT' ;
  65. ER = PROG ELIY ;
  66. RV.INCO.'IT' = (RV.INCO.'IT') ET IT ;
  67. RV.INCO.'ER' = (RV.INCO.'ER') ET ER ;
  68.  
  69. FINSI;
  70.  
  71. RV.INCO.'UNM1' = KCHT $MT 'VECT' 'SOMMET' (RV.INCO.'UN') ;
  72. as2 ama1 = 'KOPS' 'MATRIK' ;
  73. FINPROC as2 ama1 ;
  74.  
  75. **************************
  76. * DEFINITION DU MAILLAGE *
  77. **************************
  78.  
  79. option dime 2 elem qua8 ;
  80.  
  81. P1 = 0.0 (0.5*(-1.)) ;
  82. p2 = 1.0 (0.5*(-1.)) ;
  83. p3 = 1.0 0.5 ;
  84. p4 = 0.0 0.5 ;
  85.  
  86. NY = NX ;
  87.  
  88. bas = p1 d NX p2 ;
  89. cdro = p2 d NY p3 ;
  90. haut = p3 d NX p4 ;
  91. cgau = p4 d NY p1 ;
  92. cnt = bas et cdro et haut et cgau ;
  93. mt = bas cdro haut cgau daller ;
  94.  
  95. tass MT ;
  96. orienter mt ;
  97. mt= chan mt quaf ;
  98.  
  99. $mt = MODE mt 'NAVIER_STOKES' 'MACRO' ;
  100.  
  101. mt2 = mt ;
  102. mt = DOMA $mt maillage ;
  103.  
  104. ***********************
  105. * SOLUTION ANALYTIQUE *
  106. ***********************
  107.  
  108. PI = 3.141592654 ;
  109. B = 10. ;
  110. NU = 1.0/B ;
  111. ALPHA = -1.0 ;
  112. XX YY = COOR (DOMA $MT SOMMET) ;
  113.  
  114. CX = COS ((KOPS XX '*' 180.)) ;
  115. CY = COS ((KOPS YY '*' 180.)) ;
  116. SX = SIN ((KOPS XX '*' 180.)) ;
  117. SY = SIN ((KOPS YY '*' 180.)) ;
  118.  
  119. UEX = KOPS ( KOPS SX '*' SY ) '*' B ;
  120. VEX = KOPS ( KOPS CX '*' CY ) '*' B ;
  121. PEX = KOPS ( KOPS CX '*' SY ) '*' (PI*(-2.)*ALPHA) ;
  122.  
  123. UEX1 = NOMC 'UX' UEX 'NATU' 'DISCRET' ;
  124. VEX1 = NOMC 'UY' VEX 'NATU' 'DISCRET' ;
  125. UUX = KCHT $MT VECT SOMMET (UEX1 ET VEX1) ;
  126.  
  127. B1 = B * CY ;
  128. B2 = B * SX ;
  129.  
  130. XXC YYC = COOR (DOMA $MT CENTRE) ;
  131. CX = COS ((KOPS XXC '*' 180.)) ;
  132. CY = COS ((KOPS YYC '*' 180.)) ;
  133. SX = SIN ((KOPS XXC '*' 180.)) ;
  134. SY = SIN ((KOPS YYC '*' 180.)) ;
  135.  
  136. F1_1 = KOPS ( KOPS SX '*' CX ) '*' (PI*B*B) ;
  137. F1_2 = KOPS ( KOPS SX '*' SY ) '*' (2.0*PI*PI*(ALPHA+1.0)) ;
  138. F1 = KOPS F1_1 '+' F1_2 ;
  139. F2_1 = KOPS ( KOPS SY '*' CY ) '*' (PI*B*B*(-1.0)) ;
  140. F2_2 = KOPS ( KOPS CX '*' CY ) '*' (2.0*PI*PI*(1.0-ALPHA)) ;
  141. F2 = KOPS F2_1 '+' F2_2 ;
  142.  
  143. F1 = nomc 'UX' F1 'NATU' 'DISCRET' ;
  144. F2 = nomc 'UY' F2 'NATU' 'DISCRET' ;
  145. SS = KCHT $MT VECT CENTRE (F1 ET F2) ;
  146.  
  147. ***************************
  148. * TABLES EQEX ET PRESSION *
  149. ***************************
  150. RV = EQEX $MT 'ITMA' NITER 'ALFA' CFL
  151. 'ZONE' $MT 'OPER' CALCUL
  152. 'OPTI' 'SUPG' 'MMPG'
  153. 'ZONE' $MT 'OPER' 'NS' NU SS 'INCO' 'UN'
  154. 'OPTI' 'CENTREE' EFM1
  155. ZONE $MT 'OPER' 'DFDT' 1. 'UN' 'DELTAT' INCO 'UN'
  156. ;
  157.  
  158. RV = EQEX RV
  159. 'CLIM' 'UN' 'UIMP' cgau 0. 'UN' 'VIMP' cgau B1
  160. 'UN' 'UIMP' cdro 0. 'UN' 'VIMP' cdro (B1*(-1.))
  161. 'UN' 'UIMP' bas (B2*(-1.)) 'UN' 'VIMP' bas 0.
  162. 'UN' 'UIMP' haut B2 'UN' 'VIMP' haut 0. ;
  163.  
  164. RVP = EQEX 'OPTI' 'EF' KPRESS
  165. 'ZONE' $MT OPER KBBT -1. beta INCO 'UN' 'PRES'
  166. ;
  167.  
  168. rvp.'METHINV'.TYPINV=1 ;
  169. rvp.'METHINV'.IMPINV=0 ;
  170. rvp.'METHINV'.NITMAX=300;
  171. rvp.'METHINV'.PRECOND=3 ;
  172. rvp.'METHINV'.RESID =1.e-8 ;
  173. rvp.'METHINV' . 'FCPRECT'=100 ;
  174. rvp.'METHINV' . 'FCPRECI'=100 ;
  175.  
  176. RV.'PROJ' =RVP ;
  177.  
  178. **************
  179. * TABLE INCO *
  180. **************
  181. RV.INCO=TABLE INCO ;
  182. RV.INCO.'PRES'=kcht $MT SCAL KPRESS 0. ;
  183. RV.INCO.'UN'=kcht $MT VECT SOMMET (0. 0.) ;
  184. rv.inco.'UNM1' = kcht $MT 'VECT' 'SOMMET' (0. 0.) ;
  185. rv.inco.'IT' = PROG 1 ;
  186. rv.inco.'ER' = PROG 0. ;
  187.  
  188. EXEC RV ;
  189.  
  190. SI ( 'EGA' graph 'O') ;
  191. EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  192. (RV.INCO.'ER') ;
  193. DESS EVOL4 'XBOR' 0. 2000. 'YBOR' -20.0 0.0 ;
  194. TRACE PEX MT CNT 15 'TITR' 'PRESSION THEORIQUE' ;
  195. pn=elno $MT (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  196. TRACE PN MT CNT 15 'TITR' 'PRESSION' ;
  197. unx= kcht $MT scal sommet (exco 'UX' RV.INCO.'UN') ;
  198. TRACE UEX MT CNT 15 'TITR' 'UX THEORIQUE' ;
  199. TRACE UNX MT CNT 15 'TITR' 'UX' ;
  200. UNCH = VECT RV.INCO.'UN' 0.01 UX UY JAUNE ;
  201. TRACE UNCH MT CNT ;
  202. FINSI ;
  203.  
  204. un = RV.INCO.'UN' ;
  205. unx= kcht $MT scal sommet (exco 'UX' un) ;
  206. uny= kcht $MT scal sommet (exco 'UY' un) ;
  207.  
  208. ERRX = KOPS UNX '-' UEX1 ;
  209. ERRX = KOPS ERRX '*' ERRX ;
  210. ERR2 = 0. ;
  211. NPTD= NBNO (DOMA $MT SOMMET) ;
  212. REPETER BLOC1 NPTD ;
  213. P1 = (DOMA $MT SOMMET) POIN &BLOC1 ;
  214. ERR1 = 'EXTR' ERRX 'SCAL' P1 ;
  215. ERR2 = ERR2 + ERR1 ;
  216. FIN BLOC1 ;
  217. ERR2 = ERR2/NPTD ;
  218. ERR2 = ERR2 ** 0.5 ;
  219. MESSAGE 'ERREUR VITESSE UX EN NORME L2 = ' ERR2 ;
  220.  
  221. SI ( ERR2 > 6.0E-2 ) ; ERREUR 5 ; FINSI ;
  222.  
  223. pe= (exco (rv.'INCO'.'PRESSION') 'PRES') ;
  224. PEX = NOEL $MT PEX ;
  225. PPP = ((MAXI PE) + (MINI PE))/2. ;
  226. PN = KOPS PE '-' PPP ;
  227. ERRP = KOPS PEX '-' PN ;
  228. ERRP = KOPS ERRP '*' ERRP ;
  229. ERR2 = 0. ;
  230. NELD= NBEL (DOMA $MT CENTRE) ;
  231. REPETER BLOC1 NELD ;
  232. P1 = (DOMA $MT CENTRE) POIN &BLOC1 ;
  233. ERR1 = 'EXTR' ERRP 'SCAL' P1 ;
  234. ERR2 = ERR2 + ERR1 ;
  235. FIN BLOC1 ;
  236. ERR2 = ERR2/NELD ;
  237. ERR2 = ERR2 ** 0.5 ;
  238. MESSAGE 'ERREUR PRESSION EN NORME L2 = ' ERR2 ;
  239.  
  240. SI ( ERR2 > 0.80 ) ; ERREUR 5 ; FINSI ;
  241.  
  242. SI ( (MINI RV.INCO.'ER') > -4.0 ) ; ERREUR 5 ; FINSI ;
  243.  
  244. FIN ;
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  

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