* fichier : tube3D.dgibi ************************************************************************ ************************************************************************ ****************************************************************** * CALCUL DU TUBE A CHOC * * FORMULATION VF COMPRESSIBLE EXPLICITE * * DIFFERENTS SOLVEURS * * * * A. BECCANTINI TTMF MARS 1998 * * Remise à jours : JUIILETT 2001 * * Remise à jours : SEPTEMBRE 2002 (syntaxe de 'KONV' changé) * ****************************************************************** 'OPTION' 'ELEM' 'CUB8' ; 'OPTION' 'ISOV' 'SULI' ; 'OPTION' 'ECHO' 1 ; 'OPTION' 'TRAC' 'X' ; GRAPH = VRAI ; GRAPH = FAUX ; * *** Methodes possibles : * * 'VANLEER' * 'VLH ' * 'HUSVLH ' * 'GODUNOV' * 'AUSMPLUS' METO = 'VLH' ; LOGSO = VRAI ; SAFFAC = 0.7 ; NITER = 1000 ; TFINAL = 0.2 ; ************ * MAILLAGE * ************ RAF = 1 ; NX=25 '*' RAF ; NY=2 '*' RAF ; NZ=2 '*' RAF ; L = 1.0D0 ; DX = L/NX/2.0D0; H = NY*DX ; P = NZ*DX ; xD = 0.5D0*L ; xG = -1.0D0*xD ; yH = 0.5D0*H ; yB = -1.0D0*yH ; zV = 0.5D0*P ; zR = -1.0D0*zV ; A1 = xG yB zR ; A2 = 0.0D0 yB zR; A3 = xD yB zR; A4 = xD yH zR; A5 = 0.0D0 yH zR; A6 = xG yH zR; A7 = xG yB zV ; A8 = 0.0D0 yB zV; A9 = xD yB zV; A10 = xD yH zV; A11 = 0.0D0 yH zV; A12 = xG yH zV; VECTG = XG 0.0D0 zV; VECTD = XD 0.0D0 zV; xBG = xG-DX; XBD = xD+DX; B1 = xBG yB zR; B2 = xBG yH zR; B3 = xBD yB zR; B4 = xBD yH zR; B5 = xBG yB zV; B6 = xBG yH zV; B7 = xBD yB zV; B8 = xBD yH zV; * *** Rotation * ANGLE = 30.0D0; ORIG1 = 0.0D0 0.0D0 0.0D0; ORIG2 = 0.0D0 0.0D0 1.0D0; 'MESSAGE'; 'MESSAGE' (CHAIN 'ANGLE = ' ANGLE); 'MESSAGE'; DOM1 = DOM1 'TOURNER' ANGLE ORIG1 ORIG2; DOM2 = DOM2 'TOURNER' ANGLE ORIG1 ORIG2; FRONTG = FRONTG 'TOURNER' ANGLE ORIG1 ORIG2; FRONTD = FRONTD 'TOURNER' ANGLE ORIG1 ORIG2; DOMINT = DOM1 'ET' DOM2 ; 'ELIMINATION' DOMINT 1D-6; FRONT = FRONTG 'ET' FRONTD ; 'ELIMINATION' FRONT 1D-6; FRONT = FRONT 'COULEUR' 'ROUGE'; DOMTOT = DOMINT 'ET' FRONT; 'ELIMINATION' DOMTOT 1D-6; MDOMTOT = TDOMTOT . 'QUAF' ; MDOMINT = TDOMINT . 'QUAF' ; MDOM1 = TDOM1 . 'QUAF' ; MDOM2 = TDOM2 . 'QUAF' ; MFRONTG = TFRONTG . 'QUAF' ; MFRONTD = TFRONTD . 'QUAF' ; MFRONT = TFRONT . 'QUAF' ; MFRONTG 'ET' MFRONTD 'ET' MFRONT) 1.E-5 ; * ******* Creation de la ligne Utilisée pour le Post-Traitement * reliant les points centres * XINIT = XG + (0.5D0*DX); YINIT = YB + (0.5D0*DX); ZINIT = ZV - (0.5D0*DX); XFIN = XD - (0.5D0*DX) ; YFIN = YINIT ; ZFIN = ZINIT ; PINI = XINIT YINIT ZINIT; PFIN = XFIN YFIN ZFIN; IAUX = (2 * NX) - 1 ; COURB = PINI 'DROIT' IAUX PFIN; COURB = COURB 'TOURNER' ANGLE ORIG1 ORIG2; 'REPETER' BLOC1 (IAUX-1); SI (&BLOC1 'EGA' 1); SINON; FINSI; P1 = P2; 'FIN' BLOC1; ************************************************************ * CONDITIONS INITIALES ET LIMITES. * ************************************************************ gamgd = 1.4D0; * *** Etat gauche * rog = 1.0 ; ung = 0.0 ; utg = 0.0 ; uvg = 0.0 ; pg = 1.0 ; rouxg = ((ung '*' ('COS' ANGLE)) '-' (utg '*' ('SIN' ANGLE))) '*' rog ; rouyg = ((ung '*' ('SIN' ANGLE)) '+' (utg '*' ('COS' ANGLE))) '*' rog; rouzg = uvg '*' rog; recing = 0.5D0 '*' rog '*' ((ung '*' ung) '+' (utg '*' utg) '+' (uvg '*' uvg)); retg = (pg '/' (gamgd '-' 1.0)) '+' recing; * *** Etat droite * rod = 1.0D-1 ; und = 0.0D0 ; utd = 0.0; uvd = 0.0; pd = 1.0D-1; rouxd = ((und '*' ('COS' ANGLE)) '-' (utd '*' ('SIN' ANGLE))) '*' rod; rouyd = ((und '*' ('SIN' ANGLE)) '+' (utd '*' ('COS' ANGLE))) '*' rod; rouzd = uvd '*' rod; recind = 0.5D0 '*' rod '*' ((und '*' und) '+' (utd '*' utd)); retd = (pd '/' (gamgd '-' 1.0)) '+' recind; * *** gamma * * *** ro * * *** ro u, ro v * * *** ro e * ******************************************************** *** CREATION DE 'MODE' POUR GRAPHIQUER LE CHAMELEM *** ******************************************************** * **** Les debits dans le repaire tube * GNN = (GNX * ('COS' ANGLE)) '+' (GNY * ('SIN' ANGLE)); GNT = (GNY * ('COS' ANGLE)) '-' (GNX * ('SIN' ANGLE)); GNV = GNZ; * *** GRAPHIQUE DES C.I. * 'SI' GRAPH ; * *** CREATION DE CHAMELEM * 'FINSI' ; * **** Les variables primitives * VN PN = 'PRIM' 'PERFMONO' RN GN RETN GAMN ; * **** Les gradients * CACCA CHLIM COEFSCAL = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE' CACCA CHLIM COEFVECT = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE' * **** Le calcul * TPS = 0.0 ; 'TEMPS' 'ZERO' ; 'REPETER' BL1 NITER ; * **** Primitive variables * 'SI' LOGSO ; GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR' GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR' $DOMTOT RN GRADR ALR VN GRADV ALV PN GRADP ALP GAMN ; 'SINON' ; $DOMTOT RN VN PN GAMN ; 'FINSI' ; ROF VITF PF GAMF ; DT_CON = SAFFAC '*' DELTAT ; * **** The time step linked to tps * * **** Total time step * * **** Increment of the variables (convection) * RESIDU = DTMIN '*' RESIDU ; TPS = TPS '+' DTMIN ; RN = RN '+' DRN ; GN = GN '+' DGN ; RETN = RETN '+' DRETN ; 'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ; 'MESSAGE' ('CHAINE' 'ITER =' &BL1 ' TPS =' TPS) ; 'FINSI' ; 'SI' (TPS > TFINAL) ; 'QUITTER' BL1 ; 'FINSI' ; 'FIN' BL1 ; 'TEMPS' ; * ***** On calcule les variables primitive * * **** La vitesse dans le repaire tube * VNN = (VNX * ('COS' ANGLE)) '+' (VNY * ('SIN' ANGLE)); VNT = (VNY * ('COS' ANGLE)) '-' (VNX * ('SIN' ANGLE)); VNV = VNZ; GNN = VNN * RN; GNT = VNT * RN; GNV = VNV * RN; * *** GRAPHIQUE DES SOLUTIONS * 'SI' (GRAPH); * *** CREATION DE CHAMELEM * 'FINSI' ; * *** Objects evolutions * 'SI' LOGSO ; IE = 2 ; 'SINON' ; IE = 1 ; 'FINSI' ; IT = 1 ; xx yy zz= 'COORDONNEE' Courb; x0 = 'MINIMUM' lxx; x1 = 'MAXIMUM' lxx; tro = CHAINE '1D ' METO ' : RO IT ' IT ' IE ' IE ' tmps ' TPS ' elem ' 'CUB8' ; tu = CHAINE '1D ' METO ' : U IT ' IT ' IE ' IE ' tmps ' TPS ' elem ' 'CUB8' ; tv = CHAINE '1D ' METO ' : V IT ' IT ' IE ' IE ' tmps ' TPS ' elem ' 'CUB8' ; SI GRAPH; 'DESSIN' evv 'TITRE' tv 'XBOR' x0 x1; FINSI; tw = CHAINE '1D ' METO ' : W IT ' IT ' IE ' IE ' tmps ' TPS ' elem ' 'CUB8' ; SI GRAPH; 'DESSIN' evw 'TITRE' tw 'XBOR' x0 x1; FINSI; tp = CHAINE '1D ' METO ' : P IT ' IT ' IE ' IE ' tmps ' TPS ' elem ' 'CUB8' ; ls = lro '**' gamgd; ls = lp '/' ls; ts = CHAINE '1D ' METO ' : s IT ' IT ' IE ' IE ' tmps ' TPS ' elem ' 'CUB8' ; * *** Solution analytique * -0.37000 -0.35000 -0.33000 -0.31000 -0.29000 -0.27000 -0.25000 -0.23000 -0.21000 -0.19000 -0.17000 -0.15000 -0.13000 -0.11000 -9.00000E-2 -7.00000E-2 -5.00000E-2 -3.00000E-2 -1.00000E-2; 1.00000E-02 3.00000E-02 5.00000E-02 7.00000E-02 9.00000E-02 0.11000 0.13000 0.15000 0.17000 0.19000 0.21000 0.23000 0.25000 0.27000 0.29000 0.31000 0.33000 0.35000 0.37000 0.39000 0.41000 0.43000 0.45000 0.47000) ; 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99193 .85405 .79151 .69598 .64297 .57841 .51997 .45356 .41754 .37307 .33114 .29475 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .28482 .10000 .10000 .10000 .10000 .10000; 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 .99423 .89343 .84619 .77192 .72945 .67635 .62680 .56851 .53588 .49447 .45410 .41786 .40776 .40776 .40776 .40776 .40776 .40776 .40776 .40776 .40776 .40776 .20444 .20444 .20444 .20444 .20444 .20444 .20444 .20444 .20444 .10000 .10000 .10000 .10000 .10000; 0.0 0.0 6.84663E-03 .13185 .19435 .29851 .36173 .44507 .52768 .63185 .69395 .77728 .86406 .94740 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 .97167 0.0 0.0 0.0 0.0 0.0; lsa = lroa '**' gamgd; lsa = lpa '/' lsa; *PM du fait de la précision du SIN et du COS, (de l'ordre de 10-10) * l'intervalle lxx est plus petit que lxxa. * Comme IPOL fait une erreur dans ce cas, on restreint lxxa dlp = 'MAXIMUM' (ABS ( lper '-' lpa)); dlu = 'MAXIMUM' (ABS ( luer '-' lua)); dlro = 'MAXIMUM' (ABS ( lroer '-' lroa)); dl = (dlp '+' dlu '+' dlro) '/' 3.0D0; * *** Quelque DESSIN * TAB1=TABLE; TAB1.'TITRE'= TABLE ; TAB1.1='MARQ TRIB REGU'; TAB1.2='MARQ CROI REGU'; SI GRAPH; 'DESSIN' (evro 'ET' evroa) 'TITRE' tro 'XBOR' x0 x1 'LEGE' TAB1; 'DESSIN' (evp 'ET' evpa) 'TITRE' tp 'XBOR' x0 x1 'LEGE' TAB1; 'DESSIN' (evu 'ET' evua) 'TITRE' tu 'XBOR' x0 x1 'LEGE' TAB1; 'DESSIN' (evs 'ET' evsa) 'TITRE' ts 'XBOR' x0 x1 'LEGE' TAB1; FINSI; MESSAGE 'Erreur calcul du tube a choc'; MESSAGE dl; ERREUR 5; FINSI; ERRO1 = 'MAXIMUM' (lv 'ET' lw) 'ABS' ; 'SI' (ERRO1 '>' 1.0D-6); MESSAGE 'Erreur calcul du tube a choc'; MESSAGE dl; ERREUR 5; 'FINSI' ; FIN;
© Cast3M 2003 - Tous droits réservés.
Mentions légales