Télécharger la plaquette Cast3M


fr en

Télécharger tube2D.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  tube2D.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ******************************************************************
  5. * CALCUL DU TUBE A CHOC 2D                                       *
  6. * FORMULATION VF COMPRESSIBLE EXPLICITE                          *
  7. * DIFFERENTS SOLVEURS                                            *
  8. *                                                                *
  9. * A. BECCANTINI TTMF    MARS    1998                             *
  10. ******************************************************************
  11.  
  12. 'OPTION'  'DIME' 2 ;
  13. 'OPTION'  'ELEM' 'QUA4' ;
  14. 'OPTION'  'ISOV' 'SULI' ;
  15. 'OPTION'  'ECHO' 1 ;
  16. 'OPTION' 'TRAC' 'X' ;
  17.  
  18.  GRAPH  = VRAI ;
  19.  GRAPH  = FAUX ;
  20.  
  21. ************ 
  22. * MAILLAGE *
  23. ************
  24.  
  25.  RAF = 4 ;
  26.  NX=25 '*' RAF ;
  27.  NY= RAF ;
  28.  L1 = 1. ;
  29.  H1 = (L1 '/' NX) * NY ;
  30.  
  31.  A3 = 0.0 0.0 ;
  32.  A4 = 0.0 H1 ;
  33.  A3A4 = A3 'DROIT' NY A4 ;
  34.  A1A2 = A3A4 'PLUS' ((-1 * L1) 0.0) ;
  35.  A5A6 = A3A4 'PLUS' (L1 0.0) ;
  36.  DOM1 = A1A2 'REGLER' NX A3A4 ;
  37.  DOM2 = A3A4 'REGLER' NX A5A6 ;
  38.  DOMTOT = DOM1 'ET' DOM2  ;
  39.  
  40.  $DOMTOT = 'MODE' DOMTOT 'EULER'  ;
  41.  $DOM1   = 'MODE' DOM1   'EULER'  ;
  42.  $DOM2   = 'MODE' DOM2   'EULER'  ;
  43.  
  44.  TDOMTOT = 'DOMA' $DOMTOT 'VF'  ;
  45.  TDOM1   = 'DOMA' $DOM1   'VF'  ;
  46.  TDOM2   = 'DOMA' $DOM2   'VF'  ;
  47.  
  48.  MDOMTOT = TDOMTOT . 'QUAF' ;
  49.  MDOM1   = TDOM1   . 'QUAF' ;
  50.  MDOM2   = TDOM2   . 'QUAF' ;
  51.  'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) ((H1 '/' NY) '/' 100.) ;
  52.  
  53. *
  54. ******* Creation de la ligne Utilisée pour le Post-Traitement
  55. *       reliant les points centres
  56. *
  57.  
  58.  PINI  = (TDOMTOT . 'CENTRE') 'POIN' 'PROC' ((-1 * L1) 0.0) ;
  59.  PFIN  = (TDOMTOT . 'CENTRE') 'POIN' 'PROC' (L1  0.0) ;
  60.  
  61.  IAUX  = (2 * NX) - 1 ;
  62.  COURB = PINI 'DROIT' IAUX PFIN;
  63.  'ELIMINATION'  Courb ('DOMA' $DOMTOT 'CENTRE')
  64.     ((H1 '/' NY) '/' 100.) ;
  65.  
  66. ************************************************************
  67. * CONDITIONS INITIALES ET LIMITES.                         *
  68. ************************************************************
  69.  
  70.  
  71. gamgd = 1.4D0;
  72.  
  73. *
  74. *** Etat gauche 
  75. *
  76.  
  77. rog   = 1.0 ;
  78. pg    = 1.0 ;
  79. reg    = pg '/' (gamgd '-' 1.0) ;
  80.  
  81. *
  82. *** Etat droite
  83. *
  84.  
  85. rod   =  1.0D-1 ;
  86. pd    =  1.0D-1;
  87. red   = (pd '/' (gamgd '-' 1.0))  ;
  88.  
  89. *
  90.  
  91. *
  92. *** gamma
  93. *
  94.  
  95. GAMN =  'KCHT'  $DOMTOT 'SCAL' 'CENTRE' gamgd;
  96.  
  97. *
  98. *** ro
  99. *
  100.  
  101. RO_1  = 'KCHT'  $DOM1   'SCAL' 'CENTRE' rog;
  102. RO_2  = 'KCHT'  $DOM2   'SCAL' 'CENTRE' rod;
  103. RN  = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' (RO_1 'ET'  RO_2) ;
  104.  
  105. *
  106. *** ro u, ro v
  107. *
  108.  
  109. GN_1  = 'KCHT'  $DOM1    'VECT' 'CENTRE' (0.0 0.0);
  110. GN_2  = 'KCHT'  $DOM2    'VECT' 'CENTRE' (0.0 0.0);
  111. GN    = 'KCHT'  $DOMTOT  'VECT' 'CENTRE' (GN_1 'ET' GN_2);
  112.  
  113. *
  114. *** ro e 
  115. *
  116.  
  117.  
  118. RE_1  = 'KCHT'  $DOM1   'SCAL' 'CENTRE'  reg ;
  119. RE_2  = 'KCHT'  $DOM2   'SCAL' 'CENTRE'  red ;
  120. RETN = 'KCHT'  $DOMTOT 'SCAL' 'CENTRE' (RE_1 ET RE_2) ;
  121.  
  122. ********************************************************
  123. *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM   ***
  124. ********************************************************
  125.  
  126. MOD1     =  'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE';
  127.  
  128. *
  129. **** Les debits dans le repaire tube
  130. *
  131.  
  132. GNX = 'EXCO' 'UX' GN;
  133. GNY = 'EXCO' 'UY' GN;
  134.  
  135. *
  136. *** GRAPHIQUE DES C.I.
  137. *
  138. 'SI' GRAPH ;
  139. *
  140. *** CREATION DE CHAMELEM
  141. *
  142.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  143.   CHM_GN   =  'KCHA' $DOMTOT 'CHAM' GN ;
  144.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  145.   TRAC CHM_RN   MOD1 'TITR'  ('CHAINE'  'RO   at t=' 0.0);
  146.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' 0.0);
  147.   TRAC CHM_GN   MOD1 'TITR'  ('CHAINE'  'GN   at t=' 0.0);
  148. 'FINSI' ;
  149. *
  150. **** Les variables primitives
  151. *
  152.  VN PN =  'PRIM' 'PERFMONO' 
  153.             RN GN RETN GAMN ;
  154.  
  155. *
  156. **** Les gradients
  157. *
  158.  TITI CHLIM COEFG = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  159.                 ('MOTS' 'SCAL') RN  ;
  160. *                
  161. *************************************************************************
  162. *************************************************************************
  163. *********************** LE CALCUL ***************************************
  164. *************************************************************************
  165. *************************************************************************
  166.  
  167. *
  168. *** Methodes possibles :
  169. *
  170. *   'VANLEER'
  171. *   'VLH    '
  172. *   'HUSVLH '
  173. *   'GODUNOV'
  174. *   'AUSMPLUS'
  175. *   'HLLC    '
  176. *   ...
  177.  
  178.  METO =   'HLLC' ;
  179.  LOGSO = VRAI ;
  180.  LOGRK2 = VRAI ;
  181.  SAFFAC = 0.7 ;
  182.  NITER = 1000 ;
  183.  TFINAL = 0.2 ;
  184. *
  185.  TPS = 0.0 ;
  186.  ZERO = 1.0D-12 ;
  187.  'TEMPS' 'ZERO' ;
  188.  'REPETER' BL1 NITER ;
  189. *
  190. **** Primitive variables
  191. *
  192.     VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
  193.  
  194.     'SI' LOGSO ;
  195.  
  196.        GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  197.           ('MOTS' 'SCAL') RN  'GRADGEO' COEFG ;
  198.  
  199.        GRADP ALP  = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  200.           ('MOTS' 'SCAL') PN  'GRADGEO' COEFG ;
  201.  
  202.        GRADV ALV  = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  203.           ('MOTS' 'UX' 'UY') VN   'GRADGEO'  COEFG ;
  204.  
  205.        ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  2 1
  206.                             $DOMTOT
  207.                             RN    GRADR  ALR
  208.                             VN    GRADV  ALV
  209.                             PN    GRADP  ALP
  210.                             GAMN ;
  211.  
  212.     'SINON' ;
  213.  
  214.        ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  1 1 
  215.                               $DOMTOT 
  216.                               RN   
  217.                               VN   
  218.                               PN   
  219.                               GAMN    ;
  220.  
  221.     'FINSI' ;
  222.  
  223.     RESIDU DELTAT =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
  224.          $DOMTOT  ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')
  225.          ROF VITF PF GAMF ;
  226.  
  227.     DT_CON = SAFFAC '*' DELTAT ;
  228. *
  229. **** The time step linked to tps
  230. *
  231.     DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
  232. *
  233. **** Total time step
  234. *
  235.     DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
  236. *
  237. **** Increment of the variables (convection)
  238. *
  239.     SI LOGRK2 ;
  240.        RESIDU = 0.5 * DTMIN '*' RESIDU ;
  241.  
  242.        DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
  243.        DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU
  244.           ('MOTS' 'UX' 'UY') ;
  245.        DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
  246. *
  247.        RNP   = RN '+' DRN ;
  248.        GNP   = GN '+' DGN ;
  249.        RETNP = RETN '+' DRETN ;
  250.  
  251.        VNP PNP = 'PRIM' 'PERFMONO' RNP GNP RETNP GAMN ;
  252.  
  253.        'SI' LOGSO ;
  254.  
  255.           GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  256.              ('MOTS' 'SCAL') RNP  'GRADGEO' COEFG ;
  257.  
  258.           GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  259.              ('MOTS' 'SCAL') PNP  'GRADGEO' COEFG ;
  260.  
  261.           GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  262.              ('MOTS' 'UX' 'UY') VNP   'GRADGEO'  COEFG ;
  263.  
  264.           ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  2 1
  265.                             $DOMTOT
  266.                             RNP    GRADR  ALR
  267.                             VNP    GRADV  ALV
  268.                             PNP    GRADP  ALP
  269.                             GAMN ;
  270.  
  271.        'SINON' ;
  272.  
  273.           ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  1 1 
  274.                               $DOMTOT 
  275.                               RNP   
  276.                               VNP   
  277.                               PNP   
  278.                               GAMN    ;
  279.  
  280.        'FINSI' ;
  281.  
  282.        RESIDU DELTAT =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
  283.             $DOMTOT  ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')
  284.             ROF VITF PF GAMF ;
  285.     'FINSI' ;
  286.  
  287.  
  288.     RESIDU = DTMIN '*' RESIDU ;
  289.  
  290.     DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
  291.     DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU
  292.        ('MOTS' 'UX' 'UY') ;
  293.     DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
  294.  
  295.     TPS  = TPS '+' DTMIN ;
  296.     RN   = RN '+' DRN ;
  297.     GN   = GN '+' DGN ;
  298.     RETN = RETN '+' DRETN ;
  299.  
  300.     'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
  301.         'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
  302.     'FINSI' ;
  303.  
  304.     'SI' (TPS > TFINAL) ;
  305.        'QUITTER' BL1 ;
  306.     'FINSI' ;
  307.  
  308.  'FIN' BL1 ;
  309.  'TEMPS' ;
  310.  
  311. *                
  312. *************************************************************************
  313. *************************************************************************
  314. *********************** FIN  CALCUL *************************************
  315. *************************************************************************
  316. *************************************************************************
  317. *
  318.  
  319.  VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
  320.  
  321. *
  322. **** La vitesse dans le repaire tube
  323. *
  324.  
  325. VNN = 'EXCO' 'UX' VN;
  326. VNT = 'EXCO' 'UY' VN;
  327.  
  328. *
  329. *** GRAPHIQUE DES SOLUTIONS
  330. *
  331.  
  332. 'SI' (GRAPH);
  333. *
  334. *** CREATION DE CHAMELEM
  335. *
  336.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  337.   CHM_GN   =  'KCHA' $DOMTOT 'CHAM' GN ;
  338.   CHM_RETN =  'KCHA' $DOMTOT 'CHAM' RETN ;
  339.   TRAC CHM_RN  MOD1 'TITR'  ('CHAINE'  'RO at t=' TPS);
  340.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE' 'ROET at t=' TPS);
  341.   TRAC CHM_GN MOD1 'TITR'  ('CHAINE'  'RUN at t=' TPS);
  342. 'FINSI' ;
  343.  
  344. *
  345. *** Objects evolutions
  346. *
  347. 'SI' LOGSO ;
  348.    IE = 2 ;
  349. 'SINON' ;
  350.    IE = 1 ;
  351. 'FINSI' ;
  352. SI LOGRK2 ;
  353.    IT  = 2 ;
  354. 'SINON' ;
  355.    IT = 1 ;
  356. 'FINSI' ;
  357.  
  358. xx yy = 'COORDONNEE' Courb;
  359. ss =  xx ;
  360. evxx = 'EVOL' 'CHPO' ss Courb ;
  361. lxx  = 'EXTRAIRE'  evxx 'ORDO' ;
  362.  
  363. x0   = 'MINIMUM' lxx;
  364. x1   = 'MAXIMUM' lxx;
  365.  
  366. evro = 'EVOL'  'CHPO' RN Courb ;
  367. lro  = 'EXTRAIRE'  evro 'ORDO' ;
  368. evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro;
  369.  
  370. tit =  ('CHAINE' '1D ' METO  ' IT '  IT
  371.    '  IE ' IE '  TPS ' TPS) ;
  372.  
  373. evu  = 'EVOL' 'CHPO' VNN Courb ;
  374. lu   = 'EXTRAIRE'  evu 'ORDO' ;
  375. evu  = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
  376.  
  377. evv  = 'EVOL' 'CHPO' VNT Courb ;
  378. lv   = 'EXTRAIRE'  evv 'ORDO' ;
  379. evv  = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
  380.  
  381. SI GRAPH;
  382.  'DESSIN' evv 'TITRE' tit 'XBOR' x0 x1;
  383. FINSI;
  384.  
  385. evp  = 'EVOL' 'CHPO' PN Courb ;
  386. lp   = 'EXTRAIRE'  evp 'ORDO' ;
  387. evp  = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
  388.  
  389. ls = lro '**' gamgd;
  390. ls = lp '/' ls;
  391. evs  = 'EVOL' 'MANU' 'x' lxx 's' ls ;
  392. ts   = CHAINE  '1D ' METO ' : s  IT '  IT     '  IE ' IE   
  393.       '  tmps ' TPS  ' elem ' 'CUB8' ;
  394.  
  395. *
  396. *** Solution analytique
  397. *
  398.  
  399. lxxa = 'PROG' -0.49000  -0.47000  -0.45000  -0.43000  -0.41000  -0.390
  400.  -0.37000  -0.35000  -0.33000  -0.31000 -0.29000  -0.27000  -0.25000
  401.  -0.23000  -0.21000  -0.19000  -0.17000  -0.15000  -0.13000  -0.11000 
  402.  -9.00000E-2 -7.00000E-2 -5.00000E-2 -3.00000E-2 -1.00000E-2;
  403. lxxa = lxxa 'ET' ('PROG' 
  404.  1.00000E-02  3.00000E-02  5.00000E-02  7.00000E-02  9.00000E-02 
  405.  0.11000 0.13000 0.15000   0.17000  0.19000  0.21000  0.23000  0.25000
  406.  0.27000  0.29000  0.31000  0.33000  0.35000  0.37000  0.39000  0.41000
  407.  0.43000  0.45000  0.47000) ;
  408. lpa = 'PROG' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000      
  409.  1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99193
  410.  .85405 .79151 .69598 .64297 .57841 .51997 .45356 .41754
  411.  .37307 .33114 .29475 .28482 .28482 .28482 .28482 .28482
  412.  .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 
  413.  .28482 .28482 .28482 .28482 .28482 .28482 .10000 .10000
  414.  .10000 .10000 .10000;
  415. lroa = 'PROG' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000      
  416.  1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99423
  417.  .89343 .84619 .77192 .72945 .67635 .62680 .56851 .53588
  418.  .49447 .45410 .41786 .40776 .40776 .40776 .40776 .40776
  419.  .40776 .40776 .40776 .40776 .40776 .20444 .20444 .20444
  420.  .20444 .20444 .20444 .20444 .20444 .20444 .10000 .10000      
  421.  .10000 .10000 .10000;
  422. lua = 'PROG' 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
  423.  0.0 0.0 6.84663E-03 .13185 .19435 .29851 .36173 .44507 .52768
  424.  .63185 .69395 .77728 .86406 .94740 .97167 .97167 .97167
  425.  .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167
  426.  .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167
  427.  0.0 0.0 0.0 0.0 0.0;
  428.  EPSI=0.26;
  429.  
  430. lsa = lroa '**' gamgd;
  431. lsa = lpa '/' lsa;
  432.  
  433. evroa = 'EVOL' 'MANU' 'x' lxxa 'ro' lroa;
  434. evpa  = 'EVOL' 'MANU' 'x' lxxa 'p'  lpa;
  435. evua  = 'EVOL' 'MANU' 'x' lxxa 'ua' lua;
  436. evsa  = 'EVOL' 'MANU' 'x' lxxa 's'  lsa;
  437.  
  438. *PM  du fait de la précision du SIN et du COS, (de l'ordre de 10-10)
  439. *    l'intervalle lxx est plus petit que lxxa.
  440. * Comme IPOL fait une erreur dans ce cas, on restreint lxxa
  441. 'REMPLACER' lxxa 1 ('EXTR' lxx 1) ;
  442. 'REMPLACER' lxxa ('DIME' lxxa) ('EXTR' lxx ('DIME' lxx)) ;
  443.  
  444. lper  = 'IPOL' lxxa lxx lp ;
  445. luer  = 'IPOL' lxxa lxx lu ;
  446. lroer = 'IPOL' lxxa lxx lro ;
  447.  
  448. dlp = 'MAXIMUM' (ABS ( lper '-' lpa));
  449. dlu = 'MAXIMUM' (ABS ( luer '-' lua));
  450. dlro = 'MAXIMUM' (ABS ( lroer '-' lroa));
  451. dl = (dlp '+' dlu '+' dlro) '/' 3.0D0;
  452.  
  453. *
  454. *** Quelque DESSIN
  455. *
  456.  
  457.  
  458. TAB1=TABLE;
  459. TAB1.'TITRE'= TABLE ;
  460. TAB1.1='MARQ TRIB REGU';
  461. TAB1.'TITRE' . 1  = MOT 'Numerique';
  462. TAB1.2='MARQ CROI REGU';
  463. TAB1.'TITRE' . 2  = MOT 'Exacte';
  464. SI GRAPH;
  465. 'DESSIN' (evro 'ET' evroa) 'XBOR' x0 x1
  466. 'LEGE' TAB1 'TITRE' tit ;
  467. 'DESSIN' (evp 'ET' evpa) 'XBOR' x0 x1
  468. 'LEGE' TAB1 'TITRE' tit ;
  469. 'DESSIN' (evu 'ET' evua) 'XBOR' x0 x1
  470. 'LEGE' TAB1 'TITRE' tit ;
  471. 'DESSIN' (evs 'ET' evsa) 'XBOR' x0 x1
  472. 'LEGE' TAB1 'TITRE' tit ;
  473. FINSI;
  474.  
  475. SI (DL > EPSI);
  476.   MESSAGE 'Erreur calcul du tube a choc';
  477.   MESSAGE dl;
  478.   ERREUR 5;
  479. FINSI;
  480.  
  481. FIN;
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  

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