Télécharger la plaquette Cast3M


fr en

Télécharger consistence1_Godunov.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  consistence1_Godunov.dgibi
  2. ***********************************************************
  3. **** APPROCHE VF "Cell-Centred Formulation" pour la    ****
  4. **** solution des                                      **** 
  5. **** Equations d'Euler pour un gaz parfait             ****
  6. ****                                                   ****
  7. **** Gaz monoespece  "calorically perfect"             ****
  8. **** Consistence, methode GODUNOV                      ****
  9. ****                                                   ****
  10. **** A. BECCANTINI DRN/DMT/SEMT/TTMF    NOVEMBRE 1998  ****
  11. ***********************************************************
  12.  
  13. 'OPTION'  'DIME' 2 ;
  14. 'OPTION'  'ELEM' 'QUA4' ;
  15. 'OPTION'  'ISOV' 'SULI' ;
  16. 'OPTION'  'ECHO' 0 ;
  17. 'OPTION'  'TRAC' 'X';
  18.  
  19.  
  20. GRAPH = FAUX ;
  21. * GRAPH = VRAI ;
  22.  
  23.  
  24. *
  25. *** Methodes possibles :
  26. *
  27. *   'GODUNOV'
  28. *   'VANLEER'
  29. *   'VLH'
  30. *   'HUSVL'
  31. *   'HUSVLH'
  32. *
  33.  
  34. METO = 'GODUNOV' ;
  35.  
  36. ************ 
  37. * MAILLAGE *
  38. ************
  39.  
  40. NX = 5 ;
  41.  
  42. NY =  1  ;
  43.  
  44. L   = 1.0D0  ;
  45. DX  = L/NX/2.0D0 ;
  46. H   = DX  ;
  47.  
  48. xD =  0.5D0*L   ;
  49. xG = -1.0D0*xD  ;
  50. yH =  0.5D0*H   ;
  51. yB =  -1.0D0*yH ;
  52.  
  53. A1 = xG yB     ;
  54. A2 = 0.0D0 yB  ;
  55. A3 = xD yB     ;
  56. A4 = xD yH     ;
  57. A5 = 0.0D0 yH  ;
  58. A6 = xG yH     ;
  59. VECTG = XG 0.0D0 ;
  60. VECTD = XD 0.0D0 ;
  61.  
  62. xBG = xG-DX ;
  63. XBD = xD+DX ;
  64.  
  65. B1 = xBG yB ;
  66. B2 = xBG yH ;
  67. B3 = xBD yB ;
  68. B4 = xBD yH ;
  69.  
  70. BAS1  = A1 'DROI' NX A2 ;
  71. BAS2  = A2 'DROI' NX A3 ;
  72. HAU2  = A4 'DROI' NX A5 ;
  73. HAU1  = A5 'DROI' NX A6 ;
  74. LATG  = B1 'DROI' NY B2 ;
  75. LAT1  = A1 'DROI' NY A6 ;
  76. LAT12 = A2 'DROI' NY A5 ;
  77. LAT2  = A3 'DROI' NY A4 ;
  78. LATD  = B3 'DROI' NY B4 ;
  79.  
  80.  
  81. DOM1  = LAT12 'TRANSLATION' NX VECTG ;
  82.  
  83. DOM2  = LAT12 'TRANSLATION' NX VECTD ;
  84.  
  85. VECTFG = (-1.0D0*DX) 0.0D0;
  86. VECTFD =  DX 0.0D0;
  87. VECTBH = 0.0D0 (yH '-' yB);
  88. VECTHB = 0.0D0 (yB '-' yH);
  89.  
  90. FRONTG1 = LAT1 'TRANSLATION' 1 VECTFG ;
  91. DCEL1 = FRONTG1 'ET' DOM1;
  92. FRONTG2 = DCEL1 'PLUS' VECTBH;
  93. FRONTG3 = DCEL1 'PLUS' VECTHB;
  94. FRONTG  = (FRONTG1 'ET' FRONTG2 'ET' FRONTG3) 'COULEUR' 'ROUG';
  95. 'ELIMINATION' 0.0001  FRONTG ;
  96.  
  97.  
  98.  
  99. FRONTD1 =  LAT2  'TRANSLATION' 1 VECTFD;
  100. DCEL2 = FRONTD1 'ET' DOM2;
  101. FRONTD2 = DCEL2 'PLUS' VECTBH;
  102. FRONTD3 = DCEL2 'PLUS' VECTHB;
  103. FRONTD  = (FRONTD1 'ET' FRONTD2 'ET' FRONTD3) 'COULEUR' 'VERT';
  104. 'ELIMINATION' 0.0001  FRONTD ;
  105.  
  106.  
  107. *
  108. *** Rotation
  109. *
  110.  
  111. ANGLE = 30.0D0;
  112. ORIG = 0.0D0 0.0D0;
  113.  
  114. 'MESSAGE' ;
  115. 'MESSAGE'  ('CHAINE'  'ANGLE = ' ANGLE);
  116. 'MESSAGE' ;
  117.  
  118. DOM1   = DOM1   'TOURNER' ANGLE ORIG;
  119. DOM2   = DOM2   'TOURNER' ANGLE ORIG;
  120. FRONTG = FRONTG 'TOURNER' ANGLE ORIG;
  121. FRONTD = FRONTD 'TOURNER' ANGLE ORIG;
  122.  
  123. DOMINT = DOM1 'ET' DOM2 ;
  124. 'ELIMINATION' DOMINT 1D-6 ;
  125.  
  126. FRONT = FRONTG ET FRONTD ;
  127. 'ELIMINATION' FRONT 1D-6 ;
  128.  
  129. DOMTOT = DOMINT ET FRONT ;
  130. 'ELIMINATION' DOMTOT 1D-6 ;
  131.  
  132. 'SI' GRAPH;
  133.    'TRACER' DOMTOT 'TITRE' 'Maillage' ;
  134. 'FINSI' ;
  135.  
  136. $DOMTOT = 'MODE' DOMTOT 'EULER'  ;
  137. $DOMINT = 'MODE' DOMINT 'EULER'  ;
  138. $DOM1   = 'MODE' DOM1   'EULER'  ;
  139. $DOM2   = 'MODE' DOM2   'EULER'  ;
  140. $FRONTG = 'MODE' FRONTG 'EULER'  ;
  141. $FRONTD = 'MODE' FRONTD 'EULER'  ;
  142. $FRONT  = 'MODE' FRONT  'EULER'  ;
  143.  
  144. TDOMTOT = 'DOMA' $DOMTOT 'VF'  ;
  145. TDOMINT = 'DOMA' $DOMINT 'VF'  ;
  146. TDOM1   = 'DOMA' $DOM1   'VF'  ;
  147. TDOM2   = 'DOMA' $DOM2   'VF'  ;
  148. TFRONTG = 'DOMA' $FRONTG 'VF'  ;
  149. TFRONTD = 'DOMA' $FRONTD 'VF'  ;
  150. TFRONT  = 'DOMA' $FRONT  'VF'  ;
  151.  
  152. MDOMTOT = TDOMTOT . 'QUAF' ;
  153. MDOM1   = TDOM1   . 'QUAF' ;
  154. MDOM2   = TDOM2   . 'QUAF' ;
  155. MDOMINT = TDOMINT . 'QUAF' ;
  156. MFRONTG = TFRONTG . 'QUAF' ;
  157. MFRONTD = TFRONTD . 'QUAF' ;
  158. MFRONT  = TFRONT  . 'QUAF' ;
  159. 'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2 'ET' MDOMINT 'ET' MFRONTG 
  160.         'ET' MFRONTD 'ET' MFRONT) 1.E-5 ;
  161.  
  162. ************************************************************
  163. * CONDITIONS INITIALES ET LIMITES.                         *
  164. ************************************************************
  165.  
  166. gamgd = 1.4D0;
  167.  
  168. *
  169. *** Etat gauche 
  170. *
  171.  
  172. rog   = 1.11;
  173. ung   = 21.3;
  174. utg   = 11.0;
  175. pg    = 1234.7;
  176.  
  177.  
  178. rouxg  = ((ung '*' ('COS' ANGLE)) '-'
  179.           (utg '*' ('SIN' ANGLE))) '*' rog ;
  180. rouyg  = ((ung '*' ('SIN' ANGLE)) '+'
  181.           (utg '*' ('COS' ANGLE))) '*' rog;
  182.  
  183. recing = 0.5D0 '*' rog '*' ((ung '*' ung) '+' (utg '*' utg));
  184. retg  = (pg '/' (gamgd '-' 1.0)) '+' recing;
  185.  
  186. *
  187. *** Etat droite
  188. *
  189.  
  190. rod   = 1.11;
  191. und   = 21.3;
  192. utd   = 11.0;
  193. pd    = 1234.7;
  194.  
  195.  
  196. rouxd = ((und '*' ('COS' ANGLE)) '-'
  197.          (utd '*' ('SIN' ANGLE))) '*' rod;
  198. rouyd = ((und '*' ('SIN' ANGLE)) '+'
  199.          (utd '*' ('COS' ANGLE))) '*' rod;
  200.  
  201. recind = 0.5D0 '*' rod '*' ((und '*' und) '+' (utd '*' utd));
  202. retd  = (pd '/' (gamgd '-' 1.0)) '+' recind;
  203.  
  204. *
  205. *** gamma
  206. *
  207.  
  208. GAMMA =  'KCHT'  $DOMTOT 'SCAL' 'CENTRE' gamgd;
  209.  
  210. *
  211. *** ro
  212. *
  213.  
  214. RO_f1  = 'KCHT'  $FRONTG  'SCAL' 'CENTRE' rog ;
  215. RO_f2  = 'KCHT'  $FRONTD  'SCAL' 'CENTRE' rod ;
  216. RO_f   = 'KCHT'  $FRONT   'SCAL' 'CENTRE' (RO_f1 'ET' RO_f2) ; 
  217.  
  218. RO_i1  = 'KCHT'  $DOM1   'SCAL' 'CENTRE' rog;
  219. RO_i2  = 'KCHT'  $DOM2   'SCAL' 'CENTRE' rod;
  220. RO_i   = 'KCHT'  $DOMINT 'SCAL' 'CENTRE' (RO_i1 'ET' RO_i2);
  221.  
  222. RN  = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' (RO_i 'ET'  RO_f) ;
  223.  
  224.  
  225. *
  226. *** ro u, ro v
  227. *
  228.  
  229. GN_f1  = 'KCHT'  $FRONTG   'VECT' 'CENTRE' (rouxg rouyg);
  230. GN_f2  = 'KCHT'  $FRONTD   'VECT' 'CENTRE' (rouxd rouyd);
  231. GN_f   = 'KCHT'  $FRONT    'VECT' 'CENTRE' (GN_f1 'ET' GN_f2);
  232.  
  233. GN_i1  = 'KCHT'  $DOM1    'VECT' 'CENTRE' (rouxg rouyg);
  234. GN_i2  = 'KCHT'  $DOM2    'VECT' 'CENTRE' (rouxd rouyd);
  235. GN_i   = 'KCHT'  $DOMINT  'VECT' 'CENTRE' (GN_i1 'ET' GN_i2);
  236.  
  237. GN  = 'KCHT'  $DOMTOT 'VECT' 'CENTRE' (GN_i 'ET'  GN_f) ;
  238.  
  239.  
  240. *
  241. *** ro e 
  242. *
  243.  
  244. RE_f1  = 'KCHT'  $FRONTG  'SCAL' 'CENTRE'  retg ;
  245. RE_f2  = 'KCHT'  $FRONTD  'SCAL' 'CENTRE'  retd ;
  246. RE_f   = 'KCHT'  $FRONT   'SCAL' 'CENTRE' (RE_f1 'ET' RE_f2) ; 
  247.  
  248. RE_i1  = 'KCHT'  $DOM1   'SCAL' 'CENTRE'  retg ;
  249. RE_i2  = 'KCHT'  $DOM2   'SCAL' 'CENTRE'  retd ;
  250. RE_i   = 'KCHT'  $DOMINT 'SCAL' 'CENTRE' (RE_i1 'ET' RE_i2);
  251.  
  252. RETN = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' (RE_i ET RE_f) ;
  253.  
  254. ********************************************************
  255. *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM   ***
  256. ********************************************************
  257.  
  258. MOD1     =  'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE';
  259.  
  260. *
  261. *** GRAPHIQUE DES C.I.
  262. *
  263.  
  264. 'SI' GRAPH;
  265. *
  266. *** CREATION DE CHAMELEM
  267. *
  268.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN;
  269.   CHM_GN   =  'KCHA' $DOMTOT 'CHAM' GN ;
  270.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN;
  271.   TRAC CHM_RN  MOD1 'TITR'  ('CHAINE'  'RHO at t=' 0.0);
  272.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'RHOET at t=' 0.0);
  273.   TRAC CHM_GN  MOD1 'TITR'  ('CHAINE'  'RHOU at t=' 0.0);
  274. 'FINSI' ;
  275.  
  276. ZERO =1.0D-12 ;
  277. TFINAL = 100. ;
  278. NITER = 50 ;
  279. SAFFAC = 0.5 ;
  280. *
  281. RN0 = 'COPIER' RN ;
  282. GN0 = 'COPIER' GN ;
  283. RETN0 = 'COPIER' RETN ;
  284. *
  285. TPS = 0. ;
  286. 'REPETER' BL1 NITER ;
  287. *
  288. **** Primitive variables
  289. *
  290.     VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMMA ;
  291.     ROF VITF PF GAMF = 'PRET' 'PERFMONO'  1 1 
  292.                               $DOMTOT RN VN PN GAMMA;
  293.     RESIDU DELTAT =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
  294.          $DOMTOT ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')
  295.          ROF VITF PF GAMF ;
  296.     DT_CON = SAFFAC '*' DELTAT ;
  297. *
  298. **** The time step linked to tps
  299. *
  300.     DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
  301. *
  302. **** Total time step
  303. *
  304.     DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
  305. *
  306. **** Increment of the variables (convection)
  307. *
  308.     RESIDU = DTMIN '*' RESIDU ;
  309. *    
  310.     DRN = 'REDU' ('EXCO' 'RN' RESIDU 'SCAL') ('DOMA' $DOMINT 'CENTRE');
  311.     DGN = 'REDU' ('EXCO' ('MOTS' 'RUX' 'RUY') RESIDU ('MOTS' 'UX' 'UY'))
  312.           ('DOMA' $DOMINT 'CENTRE') ;
  313.     DRETN = 'REDU' ('EXCO' 'RETN' RESIDU 'SCAL') 
  314.           ('DOMA' $DOMINT 'CENTRE') ;
  315. *    
  316.     TPS  = TPS '+' DTMIN ;
  317.     RN   = RN '+' DRN ;
  318.     GN   = GN '+' DGN ;
  319.     RETN = RETN '+' DRETN ;
  320. *    
  321.     'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
  322.         'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
  323.     'FINSI' ;
  324.  
  325.     'SI' (TPS > TFINAL) ;
  326.        'QUITTER' BL1 ;
  327.     'FINSI' ;
  328.  'FIN' BL1 ;
  329.  
  330. *
  331. *** GRAPHIQUE DES SOLUTIONS
  332. *
  333.  
  334. 'SI' (GRAPH);
  335. *
  336. *** CREATION DE CHAMELEM
  337. *
  338.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  339.   CHM_GN   =  'KCHA' $DOMTOT 'CHAM' GN ;
  340.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  341.   TRAC CHM_RN  MOD1 'TITR'  ('CHAINE'  'RHO at t=' TPS);
  342.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'RHOET at t=' TPS);
  343.   TRAC CHM_GN  MOD1 'TITR'  ('CHAINE'  'RHOU at t=' TPS);
  344. 'FINSI' ;
  345. *
  346. *** Critere de control
  347. *
  348.  ERRO = 'MAXIMUM' (RN '-' RN0) 'ABS' ;
  349.  'SI' (ERRO '>' 1.0D-6) ;
  350.     'ERREUR' 5 ;
  351.  'FINSI' ;
  352.  
  353. 'FIN' ;
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  

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