Télécharger la plaquette Cast3M


fr en

Télécharger tubeaxi.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  tubeaxi.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ******************************************************************
  5. * CALCUL DU TUBE A CHOC                                          *
  6. * FORMULATION VF COMPRESSIBLE EXPLICITE                          *
  7. * DIFFERENTS SOLVEURS                                            *
  8. * 'MODE' 'AXIS'                                                  *
  9. *                                                                *
  10. * A. BECCANTINI TTMF  FEVRIER 2004                               *
  11. ******************************************************************
  12.  
  13. 'OPTION'  'DIME' 2 'ELEM' 'QUA4' 'ISOV' 'SULI'
  14.   'ECHO' 1 'TRAC' 'X' 'MODE' 'AXIS' ;
  15.  
  16. * GRAPH  = VRAI ;
  17.  GRAPH  = FAUX ;
  18.  
  19. *
  20. *** Methodes possibles :
  21. *
  22. *   'VANLEER'
  23. *   'VLH    '
  24. *   'HUSVLH '
  25. *   'GODUNOV'
  26. *   'AUSMPLUS'
  27.  
  28.  METO =   'VLH' ;
  29.  LOGSO = VRAI ;
  30.  SAFFAC = 0.5 ;
  31.  NITER = 100000 ;
  32.  TFINAL = 0.2 ;
  33.  
  34. ************ 
  35. * MAILLAGE *
  36. ************
  37.  
  38.  RAF = 1 ;
  39.  NX= 50 '*' RAF ;
  40.  
  41.  L  = 1.0D0  ;
  42.  DX = 0.5 '*' L '/' NX ;
  43.  ZERO = DX '/' 1000. ;
  44.  AM1 = ZERO (-0.5 '*' L)  ;
  45.  A0  = ZERO 0.0  ;
  46.  A1  = ZERO (0.5 '*' L)  ;
  47.  BM1 = DX (-0.5 '*' L)  ;
  48.  B0  = DX 0.0  ;
  49.  B1  = DX (0.5 '*' L) ;
  50. *
  51.  A0B0   = A0  'DROIT' 1 B0  ;
  52.  A1B1   = A1  'DROIT' 1 B1  ;
  53.  AM1BM1 = AM1 'DROIT' 1 BM1 ;
  54. *
  55.  DOM1 = A0B0 'REGLER' NX AM1BM1 ;
  56.  DOM2 = A0B0 'REGLER' NX A1B1 ;
  57. *
  58.  DOMTOT = DOM1 'ET' DOM2 ;
  59.  'ELIMINATION' DOMTOT (DX '/' 100.) ;
  60.  
  61.  $DOMTOT = 'MODE' DOMTOT 'EULER'  ;
  62.  $DOM1   = 'MODE' DOM1   'EULER'  ;
  63.  $DOM2   = 'MODE' DOM2   'EULER'  ;
  64.  
  65.  TDOMTOT = 'DOMA' $DOMTOT 'VF'  ;
  66.  TDOM1   = 'DOMA' $DOM1   'VF'  ;
  67.  TDOM2   = 'DOMA' $DOM2   'VF'  ;
  68.  
  69.  MDOMTOT = TDOMTOT . 'QUAF' ;
  70.  MDOM1   = TDOM1   . 'QUAF' ;
  71.  MDOM2   = TDOM2   . 'QUAF' ;
  72.  'ELIM' (MDOMTOT 'ET' MDOM1 'ET' MDOM2) (DX '/' 100.) ;
  73.  
  74. *
  75. ******* Creation de la ligne Utilisée pour le Post-Traitement
  76. *       reliant les points centres
  77. *
  78.  CENCOR = 'COORDONNEE' 2 (TDOMTOT . 'CENTRE') ;
  79.  PINI   = 'POIN' (CENCOR 'POIN' 'MINI') 1 ;
  80.  PFIN   = 'POIN' (CENCOR 'POIN' 'MAXI') 1 ;
  81.  IAUX = (2 * NX) - 1 ;
  82.  COURB = PINI 'DROIT' IAUX PFIN;
  83.  'ELIMINATION'  (1.0D-6) Courb ('DOMA' $DOMTOT 'CENTRE');
  84.  'SI' GRAPH ;
  85.     'TRACER' (DOMTOT 'ET' (COURB 'COULEUR' ROUG))
  86.        'TITR' '2D domain' ;
  87.   'FINSI' ;
  88.  
  89. ************************************************************
  90. * CONDITIONS INITIALES ET LIMITES.                         *
  91. ************************************************************
  92.  gamgd = 1.4D0;
  93. *
  94. *** Etat gauche 
  95. *
  96.  rog   = 1.0 ;
  97.  pg    = 1.0 ;
  98.  retg  = (pg '/' (gamgd '-' 1.0)) ;
  99. *
  100. *** Etat droite
  101. *
  102.  rod   =  1.0D-1 ;
  103.  pd    =  1.0D-1;
  104.  retd  = (pd '/' (gamgd '-' 1.0)) ;
  105. *
  106. *** gamma
  107. *
  108.  GAMN =  'KCHT'  $DOMTOT 'SCAL' 'CENTRE' gamgd;
  109. *
  110. *** ro
  111. *
  112.  RN =  ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' rog
  113.    'NATU' 'DISCRET') 'ET' ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1
  114.    'SCAL' rod 'NATU' 'DISCRET') ;
  115. *
  116. *** ro u, ro v
  117. *
  118.  GN =  ('MANUEL' 'CHPO' (TDOMTOT . 'CENTRE') 2 'UX' 0.0 'UY' 0.0
  119.         'NATU' 'DISCRET') ;
  120. *
  121. *** ro e 
  122. *
  123.  RETN =  ('MANUEL' 'CHPO' (TDOM1 . 'CENTRE') 1 'SCAL' retg
  124.    'NATU' 'DISCRET') 'ET' ('MANUEL' 'CHPO' (TDOM2 . 'CENTRE') 1
  125.    'SCAL' retd 'NATU' 'DISCRET') ;
  126. ********************************************************
  127. *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM   ***
  128. ********************************************************
  129.  MOD1     =  'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE';
  130.  VN PN =  'PRIM' 'PERFMONO' RN GN RETN GAMN ;
  131. *
  132.  'SI' GRAPH ;
  133. *
  134. *** CREATION DE CHAMELEM
  135. *
  136.    CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN ;
  137.    CHM_PN   =  'KCHA' $DOMTOT 'CHAM' PN ;
  138.    CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN ;
  139.    'TRACER' CHM_RN   MOD1 'TITR'  ('CHAINE'  'RO at t=' 0.0);
  140.    'TRACER' CHM_PN   MOD1 'TITR'  ('CHAINE'  'P  at t=' 0.0);
  141.    'TRACER' CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' 0.0);
  142.  'FINSI' ;
  143. *
  144. **** Les gradients
  145. *
  146.  GRADP ALP   COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  147.                 ('MOTS' 'SCAL') RN  ;
  148.  GRADV ALV   COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  149.                 ('MOTS' 'UX' 'UY') VN ;
  150. * Initialisation of ALP to 0 (for first order computations)
  151.  ALP = ALP '*' 0.0 ;
  152. *
  153. **** Le calcul
  154. *
  155.  TPS = 0.0 ;
  156.  ZERO = 1.0D-12 ;
  157.  DTMIN = 0.0 ;
  158.  FGEO = (TDOMTOT . 'XXSUR2D') '/' (TDOMTOT . 'XXVOLUM') ;
  159.  'TEMPS' 'ZERO' ;
  160. *
  161.  'REPETER' BL1 NITER ;
  162. *
  163. **** Primitive variables
  164. *
  165.     VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
  166.  
  167.     'SI' LOGSO ;
  168.  
  169.        GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  170.           ('MOTS' 'SCAL') RN  'GRADGEO' COEFSCAL ;
  171.  
  172.        GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
  173.           ('MOTS' 'SCAL') PN  'GRADGEO' COEFSCAL ;
  174.  
  175.        GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
  176.           ('MOTS' 'UX' 'UY') VN   'GRADGEO'  COEFVECT ;
  177.  
  178.        ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  2 1
  179.                             $DOMTOT
  180.                             RN    GRADR  ALR
  181.                             VN    GRADV  ALV
  182.                             PN    GRADP  ALP
  183.                             GAMN ;
  184.  
  185. *       ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  2 2
  186. *                            $DOMTOT
  187. *                            RN    GRADR  ALR
  188. *                            VN    GRADV  ALV
  189. *                            PN    GRADP  ALP
  190. *                            GAMN (DTMIN '/' 2) ;
  191.     'SINON' ;
  192.  
  193.        ROF VITF PF GAMF  = 'PRET' 'PERFMONO'  1 1 
  194.                               $DOMTOT 
  195.                               RN   
  196.                               VN   
  197.                               PN   
  198.                               GAMN    ;
  199.  
  200.     'FINSI' ;
  201.  
  202.     RESIDU DELTAT =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
  203.          $DOMTOT  ('MOTS' 'RN' 'RUX' 'RUY' 'RETN')
  204.          ROF VITF PF GAMF ;
  205.  
  206.     DT_CON = SAFFAC '*' DELTAT ;
  207. *
  208. **** The time step linked to tps
  209. *
  210.     DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
  211. *
  212. **** Total time step
  213. *
  214.     DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS) ;
  215. *
  216. *** AXI
  217. *   Contribution of pressure
  218. *   Warning: at first order, ALP initialised to 0
  219. *
  220. *   RESIP = ('NOMC' 'RUX' (PN '*' FGEO)) ;
  221. *
  222.     RESIP = 'FIMP' 'VF' 'AXIS' 'RESI' $DOMTOT
  223.        ('MOTS' 'RN' 'RUX' 'RUY' 'RETN') PN GRADP ALP ;
  224. *
  225. **** Increment of the variables (convection)
  226. *
  227.     RESIDU = DTMIN '*' (RESIDU '+' RESIP) ;
  228.  
  229.     DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
  230.     DGN = 'EXCO' ('MOTS' 'RUX' 'RUY') RESIDU
  231.        ('MOTS' 'UX' 'UY') ;       
  232.     DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
  233. *    
  234.     TPS  = TPS '+' DTMIN ;
  235.     RN   = RN '+' DRN ;
  236.     GN   = GN '+' DGN ;
  237.     RETN = RETN '+' DRETN ;
  238. *
  239.     'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
  240.         'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
  241.     'FINSI' ;
  242.  
  243.     'SI' (TPS '>EG' TFINAL) ;
  244.        'QUITTER' BL1 ;
  245.     'FINSI' ;
  246.  
  247.  'FIN' BL1 ;
  248.  'TEMPS' ;
  249. *
  250. ***** On calcule les variables primitive
  251. *
  252.  VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ;
  253. *
  254. **** La vitesse dans le repaire tube
  255. *
  256.  VNX = 'EXCO' 'UX' VN;
  257.  VNY = 'EXCO' 'UY' VN;
  258.  VNN = VNY ;
  259.  VNT = VNX ;
  260.  GNN = VNN * RN;
  261.  GNT = VNT * RN;
  262. *
  263. *** GRAPHIQUE DES SOLUTIONS
  264. *
  265. 'SI' (GRAPH);
  266. *
  267. *** CREATION DE CHAMELEM
  268. *
  269.   CHM_RN   =  'KCHA' $DOMTOT 'CHAM' RN;
  270.   CHM_GNN  =  'KCHA' $DOMTOT 'CHAM' GNN ;
  271.   CHM_GNT  =  'KCHA' $DOMTOT 'CHAM' GNT ;
  272.   CHM_RETN  =  'KCHA' $DOMTOT 'CHAM' RETN;
  273.   TRAC CHM_RN  MOD1 'TITR'  ('CHAINE'  'RO at t=' TPS);
  274.   TRAC CHM_RETN MOD1 'TITR'  ('CHAINE'  'ROET at t=' TPS);
  275.   TRAC CHM_GNN MOD1 'TITR'  ('CHAINE'  'ROUN at t=' TPS);
  276.   TRAC CHM_GNT MOD1 'TITR'  ('CHAINE'  'ROUT at t=' TPS);
  277. 'FINSI' ;
  278.  
  279. *
  280. *** Objects evolutions
  281. *
  282. 'SI' LOGSO ;
  283.    IE = 2 ;
  284. 'SINON' ;
  285.    IE = 1 ;
  286. 'FINSI' ;
  287. IT = 1 ;
  288. SS = 'COORDONNEE' 2 Courb;
  289.  
  290. evxx = 'EVOL' 'CHPO' ss Courb ;
  291. lxx = 'EXTRAIRE'  evxx 'ORDO' ;
  292.  
  293. x0 = 'MINIMUM' lxx;
  294. x1 = 'MAXIMUM' lxx;
  295.  
  296. evro = 'EVOL'  'CHPO' RN Courb ;
  297. lro  = 'EXTRAIRE'  evro 'ORDO' ;
  298. evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro;
  299. tro  = CHAINE  '1D ' METO  ' : RO   IT '  IT     '  IE ' IE   
  300.       '  tmps ' TPS  ' elem ' 'CUB8' ;
  301. evu  = 'EVOL' 'CHPO' VNN Courb ;
  302. lu   = 'EXTRAIRE'  evu 'ORDO' ;
  303. evu  = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
  304. tu   = CHAINE  '1D ' METO ' : U  IT '  IT     '  IE ' IE   
  305.       '  tmps ' TPS  ' elem ' 'CUB8' ;
  306. evv  = 'EVOL' 'CHPO' VNT Courb ;
  307. lv   = 'EXTRAIRE'  evv 'ORDO' ;
  308. evv  = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
  309. tv   = CHAINE  '1D ' METO ' : V  IT '  IT     '  IE ' IE   
  310.       '  tmps ' TPS  ' elem ' 'CUB8' ;
  311. SI GRAPH;
  312.  'DESSIN' evv 'TITRE' tv 'XBOR' x0 x1;
  313. FINSI;
  314.  
  315. evp  = 'EVOL' 'CHPO' PN Courb ;
  316. lp   = 'EXTRAIRE'  evp 'ORDO' ;
  317. evp  = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
  318. tp   = CHAINE  '1D ' METO ' : P  IT '  IT     '  IE ' IE   
  319.       '  tmps ' TPS  ' elem ' 'CUB8' ;
  320.  
  321. ls = lro '**' gamgd;
  322. ls = lp '/' ls;
  323. evs  = 'EVOL' 'MANU' 'x' lxx 's' ls ;
  324. ts   = CHAINE  '1D ' METO ' : s  IT '  IT     '  IE ' IE   
  325.       '  tmps ' TPS  ' elem ' 'CUB8' ;
  326.  
  327. *
  328. *** Solution analytique
  329. *
  330.  
  331. lxxa = 'PROG' -0.49000  -0.47000  -0.45000  -0.43000  -0.41000  -0.390
  332.  -0.37000  -0.35000  -0.33000  -0.31000 -0.29000  -0.27000  -0.25000
  333.  -0.23000  -0.21000  -0.19000  -0.17000  -0.15000  -0.13000  -0.11000 
  334.  -9.00000E-2 -7.00000E-2 -5.00000E-2 -3.00000E-2 -1.00000E-2;
  335. lxxa = lxxa 'ET' ('PROG' 
  336.  1.00000E-02  3.00000E-02  5.00000E-02  7.00000E-02  9.00000E-02 
  337.  0.11000 0.13000 0.15000   0.17000  0.19000  0.21000  0.23000  0.25000
  338.  0.27000  0.29000  0.31000  0.33000  0.35000  0.37000  0.39000  0.41000
  339.  0.43000  0.45000  0.47000) ;
  340. lpa = 'PROG' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000      
  341.  1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99193
  342.  .85405 .79151 .69598 .64297 .57841 .51997 .45356 .41754
  343.  .37307 .33114 .29475 .28482 .28482 .28482 .28482 .28482
  344.  .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 
  345.  .28482 .28482 .28482 .28482 .28482 .28482 .10000 .10000
  346.  .10000 .10000 .10000;
  347. lroa = 'PROG' 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000      
  348.  1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99423
  349.  .89343 .84619 .77192 .72945 .67635 .62680 .56851 .53588
  350.  .49447 .45410 .41786 .40776 .40776 .40776 .40776 .40776
  351.  .40776 .40776 .40776 .40776 .40776 .20444 .20444 .20444
  352.  .20444 .20444 .20444 .20444 .20444 .20444 .10000 .10000      
  353.  .10000 .10000 .10000;
  354. lua = 'PROG' 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
  355.  0.0 0.0 6.84663E-03 .13185 .19435 .29851 .36173 .44507 .52768
  356.  .63185 .69395 .77728 .86406 .94740 .97167 .97167 .97167
  357.  .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167
  358.  .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167
  359.  0.0 0.0 0.0 0.0 0.0;
  360.  
  361. lsa = lroa '**' gamgd;
  362. lsa = lpa '/' lsa;
  363.  
  364. evroa = 'EVOL' 'MANU' 'x' lxxa 'ro' lroa;
  365. evpa  = 'EVOL' 'MANU' 'x' lxxa 'p'  lpa;
  366. evua  = 'EVOL' 'MANU' 'x' lxxa 'ua' lua;
  367. evsa  = 'EVOL' 'MANU' 'x' lxxa 's'  lsa;
  368.  
  369. lper = 'IPOL' lxxa lxx lp ;
  370. luer = 'IPOL' lxxa lxx lu ;
  371. lroer = 'IPOL' lxxa lxx lro ;
  372.  
  373. dlp = 'MAXIMUM' (ABS ( lper '-' lpa));
  374. dlu = 'MAXIMUM' (ABS ( luer '-' lua));
  375. dlro = 'MAXIMUM' (ABS ( lroer '-' lroa));
  376. dl = (dlp '+' dlu '+' dlro) '/' 3.0D0;
  377. *
  378. *** Quelque DESSIN
  379. *
  380.  
  381.  
  382. TAB1=TABLE;
  383. TAB1.'TITRE'= TABLE ;
  384. TAB1.1='MARQ TRIB REGU';
  385. TAB1.'TITRE' . 1  = MOT 'Numerique';
  386. TAB1.2='MARQ CROI REGU';
  387. TAB1.'TITRE' . 2  = MOT 'Exacte';
  388. SI GRAPH;
  389. 'DESSIN' (evro 'ET' evroa) 'TITRE' tro 'XBOR' x0 x1
  390. 'LEGE' TAB1;
  391. 'DESSIN' (evp 'ET' evpa) 'TITRE' tp 'XBOR' x0 x1
  392. 'LEGE' TAB1;
  393. 'DESSIN' (evu 'ET' evua) 'TITRE' tu 'XBOR' x0 x1
  394. 'LEGE' TAB1;
  395. 'DESSIN' (evs 'ET' evsa) 'TITRE' ts 'XBOR' x0 x1
  396. 'LEGE' TAB1;
  397. FINSI;
  398.  
  399. SI (DLP > 1.0D-1);
  400.   MESSAGE 'Erreur calcul du tube a choc';
  401.   MESSAGE dl;
  402.   ERREUR 5;
  403. FINSI;
  404.  
  405. FIN;
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  

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