- Email: [email protected]
Least-squares computation of hypersonic flows
Computer Physics Communications 65 (1991) 57-61 North-Holland
57
Least-squares computation of hypersonic flows Charles-Henri Bruneau Laboratoire d'Analyse Num&ique, Bdtirnent 425, Unioersit~Paris-Sud, 91405 Orsay, France
The steady Euler equations are written in conservative form for the density and the components of the velocity as unknowns. The pressure is eliminated in the equations of conservation of momentum by using Bernoulli's equation for steady flows. The resulting first-order nonlinear hyperbolic system is solved globally by means of Newton linearization and a least-squares embedding. This method is stable for the subsonic as well as the supersonic regime. It is applied to compute hypersonic flows around a circular cylinder or an ellipse in 2D and a sphere in 3D. Depending on the initialization, an entropy corrector is required to capture the bow shock, then a mesh adaptation procedure is used to fit the mesh to the shock front. The solutions are very close to those available in the literature.
1. Introduction W i t h the project of b u i l d i n g a space shuttle, m a n y studies have b e e n done, these last few years in Europe, to c o m p u t e h y p e r s o n i c flows [1]. M o s t of t h e m use u p w i n d techniques to c a p t u r e strong p h e n o m e n a such as b o w shocks occuring in front of a b l u n t b o d y [2-4]. D u e to the artificial viscosity c o n t a i n e d in the scheme, the discontinuities s p r e a d out on several cells a n d a sharp representation is o b t a i n e d b y local refinement. I n this p a p e r we show that the n o n d i s s i p a t i v e m e t h o d introd u c e d in former w o r k s [5] is able to c a p t u r e strong shocks within a cell, b y c o m p u t i n g the flow a r o u n d a circular c y l i n d e r or a n ellipse in two d i m e n s i o n s a n d a r o u n d a sphere in three dimensions. O u r m e t h o d , b a s e d on N e w t o n ' s algorithm, needs to have a g o o d initial guess to converge to the right solution. However, as this is n o t easy to handle, it is sometimes necessary to use a n ent r o p y corrector to drive the iterations to the physical solution. Otherwise, if we start with a b a d initialization, the m e t h o d can converge to a n o n physical solution with negative pressures at s o m e points, specially for high M a c h n u m b e r s at infinity. T h a t is the m a i n difficulty e n c o u n t e r e d a n d we show how we c a n get rid of it b y a d d i n g a n
e n t r o p y p e n a l i z a t i o n in the least-squares functional. I n the next section we give the Euler m o d e l a n d the m a i n lines of the m e t h o d . T h e n we p r e s e n t the e n t r o p y c o r r e c t i o n a n d the a d a p t a t i o n p r o c e d u r e that fits the m e s h to the shock. F i n a l l y the n u m e r i c a l results are discussed.
2. Elder model and outline of the method I n this w o r k we are l o o k i n g for s t e a d y solutions o f Euler equations. Then, we use Bernoulli's theor e m to r e d u c e the e q u a t i o n o f c o n s e r v a t i o n o f energy to its algegraic f o r m a n d get the pressure as a function of the d e n s i t y p a n d the velocity q = (u,
U, w)T: P=
3"-3, l p( H - ½q2) = Y l P -
3'2Pq 2,
where q is the m o d u l u s of q, H is the total e n t h a l p y , 3, is the r a t i o of specific heats a n d 3'1 a n d 3'2 a r e two c o n s t a n t s equal to 6 / 7 a n d 1 / 7 respectively for the air (3' = 1.4), w h e n M = 1 is t a k e n as a reference M a c h n u m b e r ( H = 3.0).
0010-4655/91/$03.50 © 1991 - Elsevier Science Publishers B.V. (North-Holland)
C.-H. Bruneau / Least-squarescomputationof hypersonicflows
58
With this expression of the pressure we get the following Euler equations written in conservative form:
3ylpm 3x
3pv 3pw A = - ~ x + - ~ + --~- = 0, 3pu
3-/10 B = ~
3(1 - ")t2)Pu2 +
3x
3y2PW2
aT-
3T
+
= o,
3pvw + 3z
3ylp D=-3z 3(1
-
+
3puw 3x 0.
Dm=
3"y1~ 3z +
U0 given, u~+, = Wm- 0,
where 0 = 05, fi, 5, ~)T is the solution of the first-order linear hyperbolic system (A.,, B.,, Cm, Din) deduced from (A, B, C, D) at each Newton step:
3pvm
3pwm
3pmi".l 3PmO
a ~ = - ~ - - + - - ~ - + --37- + - a Z - +
+
3pm# Oz
3"}tlp B,. = ~
3PmUm
3x 3(1
+
-
3y
3PmVm 3PmWm 3y 3z = 0 ,
"h)~u 2 3x
3Om~-lUm +
Oy
3~UmVm 0"/215o 2 3x
+ 0-----f--
3T2tSw2 02(1 - yz)p,,um5 3x + 3x
3Pm5Wm +
3"~2PW 2 3y
3272PmU,n['t 3.V
3OmUrnU m 3x
--
y2)PmUm~ 3y
0(1 - 3'2) O.C 3y
3------2---+
3PrnUmU O-----y--
32Y2PmOrng 3OmUml~ 32"Y2PmWmff; 3x + 3z 3x
3"[2PU2 OPUmWrn 3"~2pU 2 aT + 3T 3z
3tSvmw m 3(1 - y2)fiw2 3y + 3z
-
+
0pro 5win 3x
+
-
32"hPmU,,fi 3z
02 T2PmVm~ 3pmfw,,, 3z + 3y
30,~um# 3OmVm# 32(1 -- y2)pmWm# 3X + 3y + 3Z
3~tl~rn 3"Y2Pmu2 3z + 3z 3PmOmW,,, 3y
3PmblmW m +
3x
3"~2[.)mO2m 3z
3(1 - T2)Pmw2 =0. 3z
Then a least-squares embedding is applied to that linear system by minimizing the functional,
1£( A .2+ B 2 + C 2 + D 2 )
J ( O ) = ~3~u,,,w,,, + 3z
2
3PmOrnWm 3"Y2PmW 2 3z + 3y =0,
This first-order nonlinear hyperbolic system (A, B, C, D) for the unknown U = (P, u, v, w) T is linearized by the Newton method as follow:
3PUrn
~
3"[2OHm 3plAmUrn 3(1 - - ",/2) p U 2 3y + 3T + 3y
+ -3"~2OmU2m 3y
3pow - - + - 3z 3y 3"y2pV 2
V2)PW 2
3z
3PmUrnUrn 3y
+ -3OmUW - m OPmUmW 32Y2PmWmff; 3YlPm 3z + 3~ 3y 3y
3y2PW2 3y =0,
3y2PU2 - - + 3z
~[2)Omu2 3x
3p,.Umg 32(1 + 3p.,fiv,,, b 3x + - 3x +
3~/2PU 2
C= m3y
--
3PmUrnW m 3"~2[,)mW2 3z + 3x =0,
m + 3PUmW - 3z
3pUU 3(1 - - " y 2 ) P V 2 3~-y--- + ~ + 3y
3"/10
3"~2PmU2 3x
3~1~ 3y
3Ou.___uu 3"/2pv2 3y 3x
+
3puw
+
3(1
dxdydz,
where ~ is the computational domain around the blunt body. The boundary conditions at the limits of the domain are the usual conditions for the supersonic regime. The whole flow is specified at the entrance section (qo~ "n < 0), nothing is given at the exit section (q~ "n > 0) and the tangency
59
C-H. Bruneau / Least-squares computation of hypersonicflows
condition on the body is imposed by adding to the functioinal J ( U ) the term
and using the conservation of mass (div we finally get, as p is always positive,
1
q'Ws>_O. fbody[(qm--q)
"n] 2 do,
where n always denotes the unit normal vector pointing outside of the domain. At last the approximation is achieved by linear quadrilateral finite elements applied to the groups of conservative variables, for instance,
D i=1
where the subscript h denotes the approximate variables, D is the number of free nodes and (%)i = 1. . . . . D is the basis of the finite dimensional space (see ref. [5]).
3. Entropy corrector and numerical improvements Applying the method above we can select bad solutions. For instance, round a cylinder with a coarse mesh of 5 × 5 cells, the computation for Mach number at infinity Mo~ = 8 gives a converged solution in 14 iterations of the Newton method; but the values of the density and the velocity are such that the pressure is negative at a few points! This is due to the fact that if the initial data U0 is far from the solution, then the Newton linearization selects a nonphysical solution. It is well-known that hyperbolic systems need an entropy inequality to avoid expansion shocks and select the good solution. For the Euler systems this inequality can be written in the form
pqs > 0
where s = 3'p/pV is the entropy we use, because it is more convenient to normalize flow at infinity; we always take soo = 1. For steady flows this inequality reduces to div
This formula completes Euler system; the difficulty is to include such a condition in the method. What is important is to control the entropy in the hypersonic region f~oo before the bow shock to avoid expansion phenomena. In this region the flow must be constant and equal to the flow at infinity; by taking qoo in the x direction, we can reduce the previous inequality to as a-~ > 0 in f~oo,
phu = E
Ops 0--t- + div
pq = O)
pqs >_0
where s is the nonlinear function of U given by
s(U)
(
3"P
u2 + U2-[-W2 )
= -p-7 = (3' - 1)P 1-v H -
2
"
Finally, when the constrain Os/ax >_0 in f~oo is unsatisfied, the equation as/Ox = 0 is linearized by Newton's method and added to the leastsquares functional as a penalization,
1csfu [Y(sm)S~l dx dy dz, where Cs is a positive constant, Sm = S ( U m ) ,
Y(sm ) and
= {1 O,
Sm is
if s,~ < s o, otherwise,
the linearized equation:
(
a(3' - ])2p£y
U2 + Um+ ) 2 Wm 2 2
Sm=--
ax a(3' - 1)p~-'(Umfi + V,.5 + W,.~) aX a(3' - ] ) 0 L - "
(
H-
Um + 2
ax The above entropy corrector is well defined and easy to incorporate in the least-squares functional. The only problem is the determination of the
60
C.-H. Bruneau / Least-squares computation of hypersonicflows
c o n s t a n t C s which is a function of the m e s h size a n d the M a c h n u m b e r at infinity. T h e n u m e r i c a l tests show that this c o n s t a n t m u s t be d e v i d e d b y thirty when the m e s h size is d i v i d e d b y two a n d m u s t be t a k e n smaller when the M a c h n u m b e r at infinity increases. F o r instance, on a coarse m e s h of 5 x 5 cells a r o u n d a c y l i n d e r in 2D, we t a k e C s = 1 0 - 2 for Mo~ = 3 a n d C s = 10 5 for M ~ = 8. In the n u m e r i c a l c o m p u t a t i o n s we use a multigrid-like p r o c e d u r e to i m p r o v e the efficiency of the m e t h o d ; we solve the p r o b l e m on a coarse mesh, then we e x t e n d the solution to a finer grid with a mesh size twice as fine as the p r e v i o u s one to give a g o o d p r e d i c t o r on this finer grid. This p r o c e d u r e is r e p e a t e d until the finest grid is reached. T h e r e m a r k s on the c o n s t a n t s Cs a b o v e show that the e n t r o p y c o r r e c t o r is p a r t i c u l a r l y i m p o r t a n t on the coarsest mesh to select the g o o d solution a n d that its influence vanishes when fine grids are concerned. A t the end, we use a n a d a p t a t i o n p r o c e d u r e to fit the mesh to the b o w shock c a p t u r e d with the help of the e n t r o p y corrector. T h e a i m of this p r o c e d u r e is to s m o o t h the shock front which is l o c a t e d on one cell a n d j u m p s f r o m cells to cells
0.4343 0.3814
P
0.0368
~ stagnation point
Fig. 2. Pressure on the axis of symmetry and on the profile for the MS = 3 flow around the cylinder.
w h e n the m e s h is n o t a d a p t e d . M o r e details can be f o u n d in ref. [6].
4. Numerical results T h e first test p r e s e n t e d in figs. 1 a n d 2 is the flow a r o u n d a circular c y l i n d e r of r a d i u s one at M ~ = 3. B o t h the d i s t a n c e b e t w e e n the b o w shock
,4 // /
Fig. 1. Adapted mesh and pressure contours for the cylinder at M~o= 3. o
x
I
shock
/
Fig. 3. Adapted mesh and pressure contours for the M~ = 3 flow around the ellipse.
C.-H. Bruneau / Least-squares computation of hypersonicflows
Fig. 4. Adapted mesh and pressure contours for the M~ = 8 flow around the ellipse.
61
tively. W e can see the b o w shock tightening to the profile as the M a c h n u m b e r at infinity increases. F i n a l l y a 3D c o m p u t a t i o n over a sphere is p r e s e n t e d for M ~ = 8. H e r e again, the b o w shock is c a p t u r e d within a cell with the g o o d distance A = 0.15 o n the axis of s y m m e t r y a n d the g o o d pressure p = 0.012 at the nose. However, in this 3 D case, no m e s h a d a p t a t i o n p r o c e d u r e is a p p l i e d a n d one c a n see s o m e d i s s i p a t i o n on the isobars p r e s e n t e d in fig. 5 d u e to a b a d shape of the shock far f r o m the axis of s y m m e t r y . W e are still d o i n g some w o r k to i m p r o v e the solutions in three dimensions.
5. Conclusion a n d the b o d y A = 0.68 a n d the shape of the shock front, are in very g o o d a g r e e m e n t with the t h e o r y [7]. Moreover, with only twenty cells, the value of the s t a g n a t i o n pressure b e h i n d the shock is given within 2% of error (exact value Ps = 0.4439); so are the R a n k i n e - H u g o n i o t conditions. T h e n the flow a r o u n d an ellipse in 2 D is comp u t e d and again the shock is r e p r e s e n t e d within a cell after the a d a p t a t i o n procedure. F i g u r e s 3 a n d 4 show the flow at M ~ = 3 a n d M ~ = 8, respec-
The least-squares e m b e d d i n g m e t h o d p r e s e n t e d a b o v e can b e a p p l i e d to h y p e r s o n i c flows a n d is able to c a p t u r e strong b o w shocks w i t h o u t the help of artificial viscosity. Both an e n t r o p y corrector a n d an a d a p t a t i o n p r o c e d u r e yield discontinuities which are r e p r e s e n t e d within a cell. T h e purp o s e of the e n t r o p y p e n a l i z a t i o n is to drive N e w ton iterates to the p h y s i c a l solution, b u t it does n o t a d d a n y viscosity in the scheme. O n the o t h e r hand, the m e s h a d a p t a t i o n i m p r o v e s the shape of the shock b y fitting the m e s h to the shock front.
References
Fig. 5. Pressure contours for the M~ = 8 flow around the sphere.
[1] Proc. Workshop on Hypersonic Flows for Reentry Problems, Antibes, 1990, INRIA-GAMNI-SMAI eds. [2] L. Fezoui and B. Stoufflet, A class of implicit upwind schemes for Euler simulations with unstructured meshes, J. Comput. Phys. 84 (1989) 174. [3] B. Van Leer, J.L. Thomas, P.L. Roe and R.W. Newsome, A comparison of numerical flux formulas for the Euler and Navier-Stokes equations, AIAA 8th Comput. Fluid Dynamics Conf., Hawai, 1987, AIAA paper 87-1104. [4] S.R. Chakravarthy, The versatility and reliability of Euler solvers based on high-accuracy TVD formulations, AIAA paper 86-0243 (1986). [5] C.H. Bruneau, J. Laminie and J.J. Chattot, Computation of 3D vortex flows past a flat plate at incidence through a variational approach of the full steady Euler equations, Int. J. Numer. Methods Fluids 9 (1989) 305. [6] C.H. Bruneau, Computation of hypersonic flows round a blunt body in 2D with Euler model, to appear in Comput. Fluids. [7] P. Diringer, Calcul d'ondes de choc d6tach6es, La Recherche A6ronautique 88 (1962) 3.
57
Least-squares computation of hypersonic flows Charles-Henri Bruneau Laboratoire d'Analyse Num&ique, Bdtirnent 425, Unioersit~Paris-Sud, 91405 Orsay, France
The steady Euler equations are written in conservative form for the density and the components of the velocity as unknowns. The pressure is eliminated in the equations of conservation of momentum by using Bernoulli's equation for steady flows. The resulting first-order nonlinear hyperbolic system is solved globally by means of Newton linearization and a least-squares embedding. This method is stable for the subsonic as well as the supersonic regime. It is applied to compute hypersonic flows around a circular cylinder or an ellipse in 2D and a sphere in 3D. Depending on the initialization, an entropy corrector is required to capture the bow shock, then a mesh adaptation procedure is used to fit the mesh to the shock front. The solutions are very close to those available in the literature.
1. Introduction W i t h the project of b u i l d i n g a space shuttle, m a n y studies have b e e n done, these last few years in Europe, to c o m p u t e h y p e r s o n i c flows [1]. M o s t of t h e m use u p w i n d techniques to c a p t u r e strong p h e n o m e n a such as b o w shocks occuring in front of a b l u n t b o d y [2-4]. D u e to the artificial viscosity c o n t a i n e d in the scheme, the discontinuities s p r e a d out on several cells a n d a sharp representation is o b t a i n e d b y local refinement. I n this p a p e r we show that the n o n d i s s i p a t i v e m e t h o d introd u c e d in former w o r k s [5] is able to c a p t u r e strong shocks within a cell, b y c o m p u t i n g the flow a r o u n d a circular c y l i n d e r or a n ellipse in two d i m e n s i o n s a n d a r o u n d a sphere in three dimensions. O u r m e t h o d , b a s e d on N e w t o n ' s algorithm, needs to have a g o o d initial guess to converge to the right solution. However, as this is n o t easy to handle, it is sometimes necessary to use a n ent r o p y corrector to drive the iterations to the physical solution. Otherwise, if we start with a b a d initialization, the m e t h o d can converge to a n o n physical solution with negative pressures at s o m e points, specially for high M a c h n u m b e r s at infinity. T h a t is the m a i n difficulty e n c o u n t e r e d a n d we show how we c a n get rid of it b y a d d i n g a n
e n t r o p y p e n a l i z a t i o n in the least-squares functional. I n the next section we give the Euler m o d e l a n d the m a i n lines of the m e t h o d . T h e n we p r e s e n t the e n t r o p y c o r r e c t i o n a n d the a d a p t a t i o n p r o c e d u r e that fits the m e s h to the shock. F i n a l l y the n u m e r i c a l results are discussed.
2. Elder model and outline of the method I n this w o r k we are l o o k i n g for s t e a d y solutions o f Euler equations. Then, we use Bernoulli's theor e m to r e d u c e the e q u a t i o n o f c o n s e r v a t i o n o f energy to its algegraic f o r m a n d get the pressure as a function of the d e n s i t y p a n d the velocity q = (u,
U, w)T: P=
3"-3, l p( H - ½q2) = Y l P -
3'2Pq 2,
where q is the m o d u l u s of q, H is the total e n t h a l p y , 3, is the r a t i o of specific heats a n d 3'1 a n d 3'2 a r e two c o n s t a n t s equal to 6 / 7 a n d 1 / 7 respectively for the air (3' = 1.4), w h e n M = 1 is t a k e n as a reference M a c h n u m b e r ( H = 3.0).
0010-4655/91/$03.50 © 1991 - Elsevier Science Publishers B.V. (North-Holland)
C.-H. Bruneau / Least-squarescomputationof hypersonicflows
58
With this expression of the pressure we get the following Euler equations written in conservative form:
3ylpm 3x
3pv 3pw A = - ~ x + - ~ + --~- = 0, 3pu
3-/10 B = ~
3(1 - ")t2)Pu2 +
3x
3y2PW2
aT-
3T
+
= o,
3pvw + 3z
3ylp D=-3z 3(1
-
+
3puw 3x 0.
Dm=
3"y1~ 3z +
U0 given, u~+, = Wm- 0,
where 0 = 05, fi, 5, ~)T is the solution of the first-order linear hyperbolic system (A.,, B.,, Cm, Din) deduced from (A, B, C, D) at each Newton step:
3pvm
3pwm
3pmi".l 3PmO
a ~ = - ~ - - + - - ~ - + --37- + - a Z - +
+
3pm# Oz
3"}tlp B,. = ~
3PmUm
3x 3(1
+
-
3y
3PmVm 3PmWm 3y 3z = 0 ,
"h)~u 2 3x
3Om~-lUm +
Oy
3~UmVm 0"/215o 2 3x
+ 0-----f--
3T2tSw2 02(1 - yz)p,,um5 3x + 3x
3Pm5Wm +
3"~2PW 2 3y
3272PmU,n['t 3.V
3OmUrnU m 3x
--
y2)PmUm~ 3y
0(1 - 3'2) O.C 3y
3------2---+
3PrnUmU O-----y--
32Y2PmOrng 3OmUml~ 32"Y2PmWmff; 3x + 3z 3x
3"[2PU2 OPUmWrn 3"~2pU 2 aT + 3T 3z
3tSvmw m 3(1 - y2)fiw2 3y + 3z
-
+
0pro 5win 3x
+
-
32"hPmU,,fi 3z
02 T2PmVm~ 3pmfw,,, 3z + 3y
30,~um# 3OmVm# 32(1 -- y2)pmWm# 3X + 3y + 3Z
3~tl~rn 3"Y2Pmu2 3z + 3z 3PmOmW,,, 3y
3PmblmW m +
3x
3"~2[.)mO2m 3z
3(1 - T2)Pmw2 =0. 3z
Then a least-squares embedding is applied to that linear system by minimizing the functional,
1£( A .2+ B 2 + C 2 + D 2 )
J ( O ) = ~3~u,,,w,,, + 3z
2
3PmOrnWm 3"Y2PmW 2 3z + 3y =0,
This first-order nonlinear hyperbolic system (A, B, C, D) for the unknown U = (P, u, v, w) T is linearized by the Newton method as follow:
3PUrn
~
3"[2OHm 3plAmUrn 3(1 - - ",/2) p U 2 3y + 3T + 3y
+ -3"~2OmU2m 3y
3pow - - + - 3z 3y 3"y2pV 2
V2)PW 2
3z
3PmUrnUrn 3y
+ -3OmUW - m OPmUmW 32Y2PmWmff; 3YlPm 3z + 3~ 3y 3y
3y2PW2 3y =0,
3y2PU2 - - + 3z
~[2)Omu2 3x
3p,.Umg 32(1 + 3p.,fiv,,, b 3x + - 3x +
3~/2PU 2
C= m3y
--
3PmUrnW m 3"~2[,)mW2 3z + 3x =0,
m + 3PUmW - 3z
3pUU 3(1 - - " y 2 ) P V 2 3~-y--- + ~ + 3y
3"/10
3"~2PmU2 3x
3~1~ 3y
3Ou.___uu 3"/2pv2 3y 3x
+
3puw
+
3(1
dxdydz,
where ~ is the computational domain around the blunt body. The boundary conditions at the limits of the domain are the usual conditions for the supersonic regime. The whole flow is specified at the entrance section (qo~ "n < 0), nothing is given at the exit section (q~ "n > 0) and the tangency
59
C-H. Bruneau / Least-squares computation of hypersonicflows
condition on the body is imposed by adding to the functioinal J ( U ) the term
and using the conservation of mass (div we finally get, as p is always positive,
1
q'Ws>_O. fbody[(qm--q)
"n] 2 do,
where n always denotes the unit normal vector pointing outside of the domain. At last the approximation is achieved by linear quadrilateral finite elements applied to the groups of conservative variables, for instance,
D i=1
where the subscript h denotes the approximate variables, D is the number of free nodes and (%)i = 1. . . . . D is the basis of the finite dimensional space (see ref. [5]).
3. Entropy corrector and numerical improvements Applying the method above we can select bad solutions. For instance, round a cylinder with a coarse mesh of 5 × 5 cells, the computation for Mach number at infinity Mo~ = 8 gives a converged solution in 14 iterations of the Newton method; but the values of the density and the velocity are such that the pressure is negative at a few points! This is due to the fact that if the initial data U0 is far from the solution, then the Newton linearization selects a nonphysical solution. It is well-known that hyperbolic systems need an entropy inequality to avoid expansion shocks and select the good solution. For the Euler systems this inequality can be written in the form
pqs > 0
where s = 3'p/pV is the entropy we use, because it is more convenient to normalize flow at infinity; we always take soo = 1. For steady flows this inequality reduces to div
This formula completes Euler system; the difficulty is to include such a condition in the method. What is important is to control the entropy in the hypersonic region f~oo before the bow shock to avoid expansion phenomena. In this region the flow must be constant and equal to the flow at infinity; by taking qoo in the x direction, we can reduce the previous inequality to as a-~ > 0 in f~oo,
phu = E
Ops 0--t- + div
pq = O)
pqs >_0
where s is the nonlinear function of U given by
s(U)
(
3"P
u2 + U2-[-W2 )
= -p-7 = (3' - 1)P 1-v H -
2
"
Finally, when the constrain Os/ax >_0 in f~oo is unsatisfied, the equation as/Ox = 0 is linearized by Newton's method and added to the leastsquares functional as a penalization,
1csfu [Y(sm)S~l dx dy dz, where Cs is a positive constant, Sm = S ( U m ) ,
Y(sm ) and
= {1 O,
Sm is
if s,~ < s o, otherwise,
the linearized equation:
(
a(3' - ])2p£y
U2 + Um+ ) 2 Wm 2 2
Sm=--
ax a(3' - 1)p~-'(Umfi + V,.5 + W,.~) aX a(3' - ] ) 0 L - "
(
H-
Um + 2
ax The above entropy corrector is well defined and easy to incorporate in the least-squares functional. The only problem is the determination of the
60
C.-H. Bruneau / Least-squares computation of hypersonicflows
c o n s t a n t C s which is a function of the m e s h size a n d the M a c h n u m b e r at infinity. T h e n u m e r i c a l tests show that this c o n s t a n t m u s t be d e v i d e d b y thirty when the m e s h size is d i v i d e d b y two a n d m u s t be t a k e n smaller when the M a c h n u m b e r at infinity increases. F o r instance, on a coarse m e s h of 5 x 5 cells a r o u n d a c y l i n d e r in 2D, we t a k e C s = 1 0 - 2 for Mo~ = 3 a n d C s = 10 5 for M ~ = 8. In the n u m e r i c a l c o m p u t a t i o n s we use a multigrid-like p r o c e d u r e to i m p r o v e the efficiency of the m e t h o d ; we solve the p r o b l e m on a coarse mesh, then we e x t e n d the solution to a finer grid with a mesh size twice as fine as the p r e v i o u s one to give a g o o d p r e d i c t o r on this finer grid. This p r o c e d u r e is r e p e a t e d until the finest grid is reached. T h e r e m a r k s on the c o n s t a n t s Cs a b o v e show that the e n t r o p y c o r r e c t o r is p a r t i c u l a r l y i m p o r t a n t on the coarsest mesh to select the g o o d solution a n d that its influence vanishes when fine grids are concerned. A t the end, we use a n a d a p t a t i o n p r o c e d u r e to fit the mesh to the b o w shock c a p t u r e d with the help of the e n t r o p y corrector. T h e a i m of this p r o c e d u r e is to s m o o t h the shock front which is l o c a t e d on one cell a n d j u m p s f r o m cells to cells
0.4343 0.3814
P
0.0368
~ stagnation point
Fig. 2. Pressure on the axis of symmetry and on the profile for the MS = 3 flow around the cylinder.
w h e n the m e s h is n o t a d a p t e d . M o r e details can be f o u n d in ref. [6].
4. Numerical results T h e first test p r e s e n t e d in figs. 1 a n d 2 is the flow a r o u n d a circular c y l i n d e r of r a d i u s one at M ~ = 3. B o t h the d i s t a n c e b e t w e e n the b o w shock
,4 // /
Fig. 1. Adapted mesh and pressure contours for the cylinder at M~o= 3. o
x
I
shock
/
Fig. 3. Adapted mesh and pressure contours for the M~ = 3 flow around the ellipse.
C.-H. Bruneau / Least-squares computation of hypersonicflows
Fig. 4. Adapted mesh and pressure contours for the M~ = 8 flow around the ellipse.
61
tively. W e can see the b o w shock tightening to the profile as the M a c h n u m b e r at infinity increases. F i n a l l y a 3D c o m p u t a t i o n over a sphere is p r e s e n t e d for M ~ = 8. H e r e again, the b o w shock is c a p t u r e d within a cell with the g o o d distance A = 0.15 o n the axis of s y m m e t r y a n d the g o o d pressure p = 0.012 at the nose. However, in this 3 D case, no m e s h a d a p t a t i o n p r o c e d u r e is a p p l i e d a n d one c a n see s o m e d i s s i p a t i o n on the isobars p r e s e n t e d in fig. 5 d u e to a b a d shape of the shock far f r o m the axis of s y m m e t r y . W e are still d o i n g some w o r k to i m p r o v e the solutions in three dimensions.
5. Conclusion a n d the b o d y A = 0.68 a n d the shape of the shock front, are in very g o o d a g r e e m e n t with the t h e o r y [7]. Moreover, with only twenty cells, the value of the s t a g n a t i o n pressure b e h i n d the shock is given within 2% of error (exact value Ps = 0.4439); so are the R a n k i n e - H u g o n i o t conditions. T h e n the flow a r o u n d an ellipse in 2 D is comp u t e d and again the shock is r e p r e s e n t e d within a cell after the a d a p t a t i o n procedure. F i g u r e s 3 a n d 4 show the flow at M ~ = 3 a n d M ~ = 8, respec-
The least-squares e m b e d d i n g m e t h o d p r e s e n t e d a b o v e can b e a p p l i e d to h y p e r s o n i c flows a n d is able to c a p t u r e strong b o w shocks w i t h o u t the help of artificial viscosity. Both an e n t r o p y corrector a n d an a d a p t a t i o n p r o c e d u r e yield discontinuities which are r e p r e s e n t e d within a cell. T h e purp o s e of the e n t r o p y p e n a l i z a t i o n is to drive N e w ton iterates to the p h y s i c a l solution, b u t it does n o t a d d a n y viscosity in the scheme. O n the o t h e r hand, the m e s h a d a p t a t i o n i m p r o v e s the shape of the shock b y fitting the m e s h to the shock front.
References
Fig. 5. Pressure contours for the M~ = 8 flow around the sphere.
[1] Proc. Workshop on Hypersonic Flows for Reentry Problems, Antibes, 1990, INRIA-GAMNI-SMAI eds. [2] L. Fezoui and B. Stoufflet, A class of implicit upwind schemes for Euler simulations with unstructured meshes, J. Comput. Phys. 84 (1989) 174. [3] B. Van Leer, J.L. Thomas, P.L. Roe and R.W. Newsome, A comparison of numerical flux formulas for the Euler and Navier-Stokes equations, AIAA 8th Comput. Fluid Dynamics Conf., Hawai, 1987, AIAA paper 87-1104. [4] S.R. Chakravarthy, The versatility and reliability of Euler solvers based on high-accuracy TVD formulations, AIAA paper 86-0243 (1986). [5] C.H. Bruneau, J. Laminie and J.J. Chattot, Computation of 3D vortex flows past a flat plate at incidence through a variational approach of the full steady Euler equations, Int. J. Numer. Methods Fluids 9 (1989) 305. [6] C.H. Bruneau, Computation of hypersonic flows round a blunt body in 2D with Euler model, to appear in Comput. Fluids. [7] P. Diringer, Calcul d'ondes de choc d6tach6es, La Recherche A6ronautique 88 (1962) 3.