Télécharger smithhutton.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : smithhutton.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. *-------------------------------------------------------------------------------
  6. * Convection/Diffusion : CAS SMITH ET HUTTON
  7. * REFERENCE : NUMERICAL HEAT TRANSFER, VOL.5, p.439, 1982
  8. *-------------------------------------------------------------------------------
  9. * Les faces latérales et supérieure d'une boite rectangulaire sont
  10. * imperméables. Le fluide rentre et sort de la boite par la face
  11. * inférieure. Il rentre par la moitié gauche de la face inférieure
  12. * et sort par la moitié droite de cette même face.
  13. *
  14. * Le champ de vitesse est connu et donné par :
  15. * v(x,y) = 2y(1-x2) i - 2x(1-y2) j
  16. * la taille de la boite étant [-1,1]x[0,1]
  17. *
  18. * L'écoulement transporte un scalaire passif qui est injecté en entrée
  19. * suivant le profil :
  20. * c(x,0) = 1 + tanh(10(2x+1))
  21. * avec x variant de -1 à 0 (entrée du domaine).
  22. * Sur les frontières imperméables du domaine, la concentration est
  23. * constante et égale à 1-tanh(10).
  24. *
  25. * On cherche la solution stationnaire du transport par diffusion et
  26. * convection du champ scalaire passif. On compare la solution en sortie
  27. * avec la solution de référence. Plusieurs solutions sont calculées
  28. * suivant le Peclet (rapport entre la convection et la diffusion).
  29. *-------------------------------------------------------------------------------
  30. *
  31. 'OPTI' 'DIME' 2 'ISOV' 'SULI' ;
  32. *
  33. * Pe : Peclet (convection sur diffusion)
  34. * Les cas traités dans la référence sont 10, 100, 500, 1000 et 1000000
  35. Pe = 10 ;
  36. *
  37. * KSUPG : Option de décentrement (CENTREE/SUPG/SUPGDC)
  38. * GRAPH : Booleen pour l'affichage des tracés à l'issue du calcul
  39. * POST1 : Booleen pour affichage Résidu à chaque pas
  40. * POST2 : Booleen pour affichage C(x,t) à chaque pas
  41. * COMPLET : Booleen modifiant la finesse du maillage et la précision des calculs
  42. * NX : Nombre de maille suivant x
  43. * NY : Nombre de maille suivant y
  44. * nbiter : Nombre de pas de temps
  45. * tolera : Tolérance pour la comparaison avec la solution de référence
  46. KSUPG = 'SUPG' ;
  47. GRAPH = FAUX ;
  48. POST1 = FAUX ;
  49. POST2 = FAUX ;
  50. COMPLET = faux ;
  51. 'OPTI' 'ELEM' 'QUA4' ;
  52. 'SI' ( COMPLET ) ;
  53. NY = 20 ;
  54. NX = 2 * NY ;
  55. nbiter = 2000 ;
  56. tolera = 0.01 ;
  57. 'SINON' ;
  58. NY = 10 ;
  59. NX = 2 * NY ;
  60. nbiter = 100 ;
  61. tolera = 0.2 ;
  62. 'FINSI' ;
  63. *
  64. *==================================================================
  65. * Calcul du résidu en température et arrêt suivant un critère
  66. *==================================================================
  67. * E/ : RVX : TABLE : TABLE des données créées par EQEX
  68. * ARG1 : Fréquence d'impression
  69. * ARG2 : Critère d'arrêt
  70. * /S : MAT1 : MATRIK : Objet vide
  71. * /S : CHP1 : CHPO : Objet vide
  72. *------------------------------------------------------------------
  73. *
  74. 'DEBPROC' residu rvx*table ;
  75. RV = rvx . 'EQEX' ;
  76. FREQ = RVX . 'ARG1' ;
  77. EPS0 = RVX . 'ARG2' ;
  78. NITER = RV . 'NITER' ;
  79. DD = RV . 'PASDETPS' . 'NUPASDT' ;
  80. NN = DD '/' FREQ ;
  81. L0 = 'EGA' (DD '-' (FREQ*NN)) 0 ;
  82. 'SI' L0 ;
  83. RANG0 = RV . 'PASDETPS' . 'NUPASDT' ;
  84. TIME0 = RV . 'PASDETPS' . 'TPS' ;
  85. CN0 = RV . 'INCO' . 'CN' ;
  86. CNM0 = RV . 'INCO' . 'CN2' ;
  87. ERR0 = ('MAXIMUM' ('ABS' (CN0 '-' CNM0))) '+' 1.D-20 ;
  88. ERR10 = ('LOG' ERR0 ) '/' ('LOG' 10.) ;
  89. 'MESSAGE' 'Résidu en Température au pas'
  90. RANG0 '(' TIME0 ') :' ERR0 ':' ERR10 ;
  91. RV . 'INCO' . 'IT' = RV . 'INCO' . 'IT' 'ET' ('PROG' RANG0) ;
  92. RV . 'INCO' . 'TI' = RV . 'INCO' . 'TI' 'ET' ('PROG' TIME0) ;
  93. RV . 'INCO' . 'ER' = RV . 'INCO' . 'ER' 'ET' ('PROG' ERR10) ;
  94. EV1 = 'EVOL' 'MANUEL' (RV . 'INCO' . 'IT') (RV . 'INCO' . 'ER') ;
  95. Y1 = ('LOG' EPS0) '/' ('LOG' 10) ;
  96. 'SI' POST1 ;
  97. X1 = 0. ; X2 = RV . 'ITMA' ;
  98. 'DESSIN' EV1 'YBOR' Y1 0. 'NCLK'
  99. 'TITR' 'Evolution du résidu' ;
  100. 'FINSI' ;
  101. 'SI' POST2 ;
  102. L1 = (PROG 0. PAS 100. 2000.)* 1.D-3 ;
  103. trace L1 cn0 DOMTOT CNT1 'TITR' 'Concentration' 'NCLK' ;
  104. 'FINSI' ;
  105. 'SI' ((ERR10 < Y1) 'ET' (DD > 10)) ;
  106. RV . 'TFINAL' = RV . 'PASDETPS' . 'TPS' ;
  107. 'FINSI' ;
  108. 'FINSI' ;
  109. RV . 'INCO' . 'CN2' = 'COPIER' RV . 'INCO' . 'CN' ;
  110. mat1 chp1 = 'KOPS' 'MATRIK' ;
  111. 'FINP' mat1 chp1 ;
  112. *------------------------------------------------------------------
  113. *
  114. *
  115. *- MAILLAGE
  116. *
  117. *
  118. * Points
  119. A1 = -1.0 0.0 ;
  120. A2 = 1.0 0.0 ;
  121. A3 = 1.0 1.0 ;
  122. A4 = -1.0 1.0 ;
  123. A0 = 0.0 0.0 ;
  124. *
  125. * Lignes
  126. LIN = A1 'DROI' NY A0 ;
  127. LOUT = A0 'DROI' NY A2 ;
  128. FBAS = LIN 'ET' LOUT ;
  129. FDRO = A2 'DROI' NY A3 ;
  130. FHAU = A3 'DROI' NX A4 ;
  131. FGAU = A4 'DROI' NY A1 ;
  132. LIMP = FDRO 'ET' FHAU 'ET' FGAU ;
  133. *
  134. * Maillage
  135. DOMTOT = 'DALL' FBAS FDRO FHAU FGAU 'PLAN' ;
  136. CNT1 = 'CONT' DOMTOT ;
  137. *
  138. * Modèles et sous-modèles
  139. DOM2 = 'CHAN' 'QUAF' DOMTOT ;
  140. LIN2 = 'CHAN' 'QUAF' LIN ;
  141. LIMP2 = 'CHAN' 'QUAF' LIMP ;
  142. $DOMTOT = 'MODE' DOM2 'NAVIER_STOKES' 'LINE' ;
  143. $LIN = 'MODE' LIN2 'NAVIER_STOKES' 'LINE' ;
  144. $LIMP = 'MODE' LIMP2 'NAVIER_STOKES' 'LINE' ;
  145. *
  146. * Récupération des maillages et fusion des supports
  147. DOMTOT = 'DOMA' $DOMTOT 'MAILLAGE' ;
  148. LIN = 'DOMA' $LIN 'MAILLAGE' ;
  149. LIMP = 'DOMA' $LIMP 'MAILLAGE' ;
  150. FDOMTOT = 'DOMA' $DOMTOT 'FACE' ;
  151. CLIN = 'DOMA' $LIN 'CENTRE' ;
  152. CLIMP = 'DOMA' $LIMP 'CENTRE' ;
  153. 'ELIM' FDOMTOT (CLIN 'ET' CLIMP) 1.D-3 ;
  154. *
  155. * Champ de vitesse
  156. XX YY = 'COOR' DOMTOT ;
  157. VXSH = (2.*YY)*(1.0-(XX*XX)) ;
  158. VYSH = (-2.*XX)*(1.0-(YY*YY)) ;
  159. VX0 = 'NOMC' 'UX' VXSH ;
  160. VY0 = 'NOMC' 'UY' VYSH ;
  161. VX = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UX' VX0 ;
  162. VY = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 'COMP' 'UY' VY0 ;
  163. CHVIT = 'KCHT' $DOMTOT 'VECT' 'SOMMET' 'COMP' 'UX' 'UY' (VX 'ET' VY) ;
  164. *
  165. * Diffusion (1/Pe)
  166. DIF = 1.0 / ('FLOT' Pe) ;
  167. *
  168. * Profil de concentration à l'entrée
  169. XBAS = 'COOR' 1 LIN ;
  170. TOTO = 2.0*XBAS + 1.0 * 10. ;
  171. SOLUTION = 'TANH' TOTO ;
  172. SOLUTION = 1.0 + SOLUTION ;
  173. CHP1 = 'KCHT' $LIN 'SCAL' 'SOMMET' 0. SOLUTION ;
  174. CHP1 = 'NOMC' 'CN' CHP1 ;
  175. *
  176. * Conditions aux limites en concentration sur les frontières imperméables
  177. C1 = (1.0 - (TANH 10.0)) ;
  178. *
  179. * Description du problème de transport
  180. RV1 = 'EQEX' $DOMTOT 'ITMA' nbiter 'ALFA' 0.7
  181. 'ZONE' $DOMTOT 'OPER' residu 1 1.D-7
  182. 'OPTI' KSUPG 'EFM1' 'EXPL'
  183. 'ZONE' $DOMTOT 'OPER' 'TSCAL' DIF 'VITESSE' 0. 'INCO' 'CN'
  184. 'OPTI' 'CENTREE'
  185. 'ZONE' $DOMTOT 'OPER' 'DFDT' 1. 'CN' 'DELTAT' 'INCO' 'CN' ;
  186. *
  187. * Description des conditions aux limites
  188. RV1 = 'EQEX' RV1
  189. 'CLIM' 'CN' 'TIMP' LIN CHP1
  190. 'CLIM' 'CN' 'TIMP' LIMP C1 ;
  191. *
  192. * Description des conditions initiales
  193. RV1 . 'INCO' = TABLE 'INCO' ;
  194. RV1 . 'INCO' . 'CN' = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 0. ;
  195. *
  196. * Autres data
  197. RV1 . 'INCO' . 'VITESSE' = CHVIT ;
  198. RV1 . 'INCO' . 'CN2' = 'KCHT' $DOMTOT 'SCAL' 'SOMMET' 0. ;
  199. RV1 . 'INCO' . 'IT' = 'PROG' ;
  200. RV1 . 'INCO' . 'TI' = 'PROG' ;
  201. RV1 . 'INCO' . 'ER' = 'PROG' ;
  202. *
  203. *
  204. *- CALCUL
  205. *
  206. *
  207. EXEC RV1 ;
  208. *
  209. *
  210. *- ANALYSE DES RESULTATS
  211. *
  212. *
  213. EVOL1 = 'EVOL' 'CHPO' (RV1 . 'INCO' . 'CN') 'SCAL' FBAS ;
  214. EVOL2 = 'EVOL' 'CHPO' XX 'SCAL' FBAS ;
  215. LIX = 'EXTR' EVOL2 'ORDO' ;
  216. LIU = 'EXTR' EVOL1 'ORDO' ;
  217. EVOL3 = 'EVOL' 'MANU' 'X' LIX 'C(X,0)' LIU ;
  218. *
  219. * Solutions de référence
  220. LIXT = PROG 0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ;
  221. SI ( Pe 'EGA' 10 ) ;
  222. LIUT = PROG 1.989 1.402 1.146 0.946 0.775 0.621 0.480 0.349
  223. 0.227 0.111 0.000 ;
  224. FINSI ;
  225. SI ( Pe 'EGA' 100 ) ;
  226. LIUT = PROG 2.000 1.940 1.836 1.627 1.288 0.869 0.480 0.209
  227. 0.070 0.017 0.000 ;
  228. FINSI ;
  229. SI ( Pe 'EGA' 500 ) ;
  230. LIUT = PROG 2.000 2.000 1.998 1.965 1.702 0.947 0.242 0.023
  231. 0.001 0.000 0.000 ;
  232. FINSI ;
  233. SI ( Pe 'EGA' 1000 ) ;
  234. LIUT = PROG 2.000 2.000 2.000 1.985 1.841 0.951 0.154 0.001
  235. 0.000 0.000 0.000 ;
  236. FINSI ;
  237. =div st}le="font: normal normal 1em/1.2em monospace; margin:0; padding:0; background:none; vertical-align:top;">SI ( Pe 'EGA' 1000000 ) ;
  • LIUT = PROG 2.000 2.000 2.000 1.999 1.964 1.000 0.036 0.001
  • 0.000 0.000 0.000 ;
  • FINSI ;
  •  
  • EVOL4 = EVOL 'MANU' 'X' LIXT 'Uref(X,0)' LIUT ;
  •  
  • LIUC = IPOL LIXT LIX LIU ;
  • NP = DIME LIXT ;
  • ERR0 = 0. ;
  • REPETER BLOC1 NP ;
  • UCAL = EXTRAIRE LIUC &BLOC1 ;
  • UREF = EXTRAIRE LIUT &BLOC1 ;
  • ERR0 = ERR0 + ((UCAL-UREF)*(UCAL-UREF)) ;
  • FIN BLOC1 ;
  • ERR0 = ERR0/NP ;
  • ERR0 = ERR0 '**' 0.5 ;
  • MESSAGE 'ERREUR ' ERR0 ;
  • *
  • *
  • *- POST-TRAITEMENT
  • *
  • *
  • 'SI' GRAPH ;
  • *
  • * Maillage
  • 'TRAC' DOMTOT 'TITR' 'Maillage' ;
  • *
  • * Vitesse
  • UNCH = 'VECT' CHVIT 0.1 'UX' 'UY' 'ROUGE' ;
  • 'TRAC' UNCH DOMTOT CNT1 'TITR' 'Vitesse transportante' ;
  • *
  • * Concentration
  • 'TRAC' DOMTOT (RV1 . 'INCO' . 'CN') CNT1 'TITR' 'Concentration' ;
  • *
  • * Convergence vers la solution stationnaire
  • EV1 = 'EVOL' 'MANU' 'ITERATIONS' (RV1 . 'INCO' . 'IT')
  • 'LOG|E|inf' (RV1 . 'INCO' . 'ER') ;
  • 'DESS' EV1 'XBOR' 0. ('FLOT' NBITER) 'YBOR' -20.0 0.0
  • 'TITR' 'Evolution du résidu' ;
  • *
  • * Comparaison avec la solution analytique
  • TAB1 = TABLE ;
  • TAB1 . TITRE = TABLE ;
  • TAB1. TITRE . 1 = 'MOT' 'CAST3M' ;
  • TAB1. TITRE . 2 = 'MOT' 'REFERENCE' ;
  • TAB1. 2 = 'MARQ LOSA NOLI' ;
  • 'TITR' 'Comparaison CAST3M/Référence' ;
  • 'DESS' (EVOL3 ET EVOL4) 'LEGE' TAB1
  • 'MIMA' 'XBOR' -1.0 1.0 'YBOR' -1.0 3.0 ;
  • *
  • FINSI ;
  • *
  • * Arret si problème
  • *
  • SI ( ERR0 > tolera ) ;
  • 'MESS' 'Erreur de ' err0 ' > ' tolera ;
  • ERREUR 5 ;
  • FINSI ;
  • 'MESS' 'Erreur de ' err0 ' < ' tolera ;
  • FIN ;
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales