Télécharger la plaquette Cast3M


fr en

Télécharger tube3Daxi.dgibi

Retour à la liste

Numérotation des lignes : 

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

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