Télécharger la plaquette Cast3M


fr en

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. PPI = 'POIN' ('DOMA' $mt KPRESS) 'PROC' ('BARY' mt) ;
  165. MPI = 'MANU' 'POI1' PPI ;
  166. RVP = EQEX 'OPTI' 'EF' KPRESS
  167.  'ZONE' $MT  OPER  KBBT  -1. beta  INCO 'UN' 'PRES'
  168.  'CLIM' 'PRES' 'TIMP' MPI 0.
  169. ;
  170.  
  171.     rvp.'METHINV'.TYPINV=1 ;
  172.     rvp.'METHINV'.IMPINV=0 ;
  173.     rvp.'METHINV'.NITMAX=300;
  174.     rvp.'METHINV'.PRECOND=3 ;
  175.     rvp.'METHINV'.RESID  =1.e-8 ;
  176.     rvp.'METHINV' . 'FCPRECT'=100 ;
  177.     rvp.'METHINV' . 'FCPRECI'=100 ;
  178.  
  179.   RV.'PROJ'     =RVP ;
  180.  
  181. **************
  182. * TABLE INCO *
  183. **************
  184. RV.INCO=TABLE INCO ;
  185. RV.INCO.'PRES'=kcht $MT SCAL KPRESS 0. ;
  186. RV.INCO.'UN'=kcht $MT VECT SOMMET (0. 0.) ;
  187. rv.inco.'UNM1' = kcht $MT 'VECT' 'SOMMET' (0. 0.) ;
  188. rv.inco.'IT' = PROG 1 ;
  189. rv.inco.'ER' = PROG 0. ;
  190.  
  191. EXEC RV ;
  192.  
  193. SI ( 'EGA' graph 'O') ;
  194. EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  195.            (RV.INCO.'ER') ;
  196. DESS EVOL4 'XBOR' 0. 2000. 'YBOR' -20.0 0.0 ;
  197. TRACE PEX MT CNT 15 'TITR' 'PRESSION THEORIQUE' ;
  198. pn=elno $MT (exco (rv.'INCO'.'PRESSION') 'PRES') kpress ;
  199. TRACE PN MT CNT 15 'TITR' 'PRESSION' ;
  200. unx= kcht $MT scal sommet (exco 'UX' RV.INCO.'UN') ;
  201. TRACE UEX MT CNT 15 'TITR' 'UX THEORIQUE' ;
  202. TRACE UNX MT CNT 15 'TITR' 'UX' ;
  203. UNCH = VECT RV.INCO.'UN' 0.01 UX UY JAUNE ;
  204. TRACE UNCH MT CNT ;
  205. FINSI ;
  206.  
  207. un = RV.INCO.'UN' ;
  208. unx= kcht $MT scal sommet (exco 'UX' un) ;
  209. uny= kcht $MT scal sommet (exco 'UY' un) ;
  210.  
  211. ERRX = KOPS UNX '-' UEX1 ;
  212. ERRX = KOPS ERRX '*' ERRX ;
  213. ERR2 = 0. ;
  214. NPTD= NBNO (DOMA $MT SOMMET) ;
  215. REPETER BLOC1 NPTD ;
  216. P1 = (DOMA $MT SOMMET) POIN &BLOC1 ;
  217. ERR1 = 'EXTR' ERRX 'SCAL' P1 ;
  218. ERR2 = ERR2 + ERR1 ;
  219. FIN BLOC1 ;
  220. ERR2 = ERR2/NPTD ;
  221. ERR2 = ERR2 ** 0.5 ;
  222. MESSAGE 'ERREUR VITESSE UX EN NORME L2 = ' ERR2 ;
  223.  
  224.  SI ( ERR2 > 6.0E-2 ) ; ERREUR 5 ; FINSI ;
  225.  
  226. pe= (exco (rv.'INCO'.'PRESSION') 'PRES')  ;
  227. PEX = NOEL $MT PEX ;
  228. PPP = ((MAXI PE) + (MINI PE))/2. ;
  229. PN = KOPS PE '-' PPP ;
  230. ERRP = KOPS PEX '-' PN ;
  231. ERRP = KOPS ERRP '*' ERRP ;
  232. ERR2 = 0. ;
  233. NELD= NBEL (DOMA $MT CENTRE) ;
  234. REPETER BLOC1 NELD ;
  235. P1 = (DOMA $MT CENTRE) POIN &BLOC1 ;
  236. ERR1 = 'EXTR' ERRP 'SCAL' P1 ;
  237. ERR2 = ERR2 + ERR1 ;
  238. FIN BLOC1 ;
  239. ERR2 = ERR2/NELD ;
  240. ERR2 = ERR2 ** 0.5 ;
  241. MESSAGE 'ERREUR PRESSION EN NORME L2 = '  ERR2 ;
  242.  
  243.  SI ( ERR2 > 0.80 ) ; ERREUR 5 ; FINSI ;
  244.  
  245.  SI ( (MINI RV.INCO.'ER') > -4.0 ) ; ERREUR 5 ; FINSI ;
  246.  
  247. FIN ;
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  

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