• Author / Uploaded
• Barbara Schlott

Citation preview

Hypersonic Aerodynamics Barbara Schlottfeldt Maia Engineering Central College of Engineering Swansea University Swansea UK [email protected] Abstract— There is a currently increasing interest in the design and manufacture of civilian hypersonic vehicles, driving the research of methods to employ in the design process. In the present work, a MATLAB program was developed using an approximate and inverse method known as Maslen’s [1] method to solve the entire inviscid hypersonic flow field of different shock waves examples. The results are discussed and compared in this paper along with the exact solution for one of these examples. This work aims to validate the approach of Maslen’s method [1] for both two dimensional and axisymmetric flow in various conditions by using non dimensional variables in equations. Keywords—hypersonic vehicles; shock wave; inverse method;

I.

Figure 1.

INTRODUCTION

SKYLON space plane

The increased research of technologies in aerodynamics field lead to the design of civilian hypersonic vehicles such as the SKYLON (Figure 1) [2] and Virgin Galactic’s SpaceShipTwo (Figure 2) [3]. SKYLON is a design for a single-stage-to-orbit space plane(STTO) by the British company Reaction Engines Limited(REL), able to achieve orbit without staging. The vehicle design is for a hydrogen-fuelled aircraft that would take off from a purpose-built runway, and accelerate to Mach 5.4 at 26 kilometers (16 mi) altitude using the atmosphere's oxygen before switching the engines to use the internal liquid oxygen (LOX) supply to take it into orbit. Once in orbit it would release its payload (of up to 15 tonnes). The vehicle will be unpiloted, but also be certified to carry passengers. SpaceShipTwo is the second stage of a two stage to orbit (TSTO) reusable, winged spacecraft, designed to repeatedly

carry as many as eight people(including two pilots) into space. SpaceShipTwo is powered by a hybrid rocket motor and this spacecraft uses a unique system to safely re-enter the Earth’s atmosphere. With this system, SpaceShipTwo can mimic the performance of a capsule or of a winged vehicle at the appropriate parts of its trajectory. SpaceShipTwo is currently undergoing testing before becoming available for commercial use. A present dilemma is encountered in the aerodynamic design of these hypersonic vehicles i.e. there is a limited number of large working section wind tunnels worldwide for testing the design, besides the ones available are very expensive to operate. CFD (Computational Fluid Dynamics) provides an alternative design approach, however the high accuracy techniques applied in CFD are themselves expensive as well and are generally employed only for final detailed design. In the initial design phase of hypersonic vehicles, cheaper, low accuracy methods are more desirable considering these methods allow a rapid sweeping of parameter spaces. Maslen’s method [1] for example is among this class of methods and produces satisfying results and simple approach when solving hypersonic flow past smooth symmetric bodies in both axisymmetric or planar cases.

Figure 2.

SpaceShipTwo by Virgin Galactic

II.

LITERATURE SURVEY

The method chosen in the present work was Maslen’s [1] and it consists in a simple inverse problem method considering the shock shape is assumed and the calculations start right behind

the shock, followed by the computation of the flow field and determination of the corresponding body shape. The method relies on the assumption that the layer between the shock and the body is very thin, creating a thin shock layer in which the flow is approximately parallel to the shock, with exception of the stagnation region near the stagnation point. Such approximations are the basis of the thin shock layer theory and this method relies strongly on this analysis. III.

METHOD

Maslen’s method [1] solves an inverse problem. A shock wave shape is assumed and the body shape that supports this shock is calculated along with the flow field between the shock and the body. The steps taken are derived from Anderson’s book[4] and Maslen’s scientific article[1].

Figure 3.

1u1   2u2 ( 0 ) p1  1u12  p2   2u2 2 ( 0 )

A. Nomenclature



1 1 h1  u12  h2  u2 2 2 2 (0)

Subscripts :  freestream value

 , x = partial derivatives with respect to  or x

Continuity equation ahead of shock is defined as in Eq.(4) :

  (  u1 ) + (  v1 )  0 x y

Dimensional variables - marked with an over bar: x , y = distance along and normal(inward) to the shock p = pressure  = density u, v = velocities in x, y direction  = stream function h

 enthalpy

L

= characteristic length

V∞ 

= freestream velocity = shock curvature

Dimensionless variables : 1 h h / ( V2 ) 2 X X/L

Flow variables through shock

The stream function across a normal shock can be found with Equations (5) and (6) :

1v1   Since

 a 0 x (0)

1,u1 are constant, the stream function across a normal shock is defined as :

 a  1u1 y ( 0 ) Behind the shock :

 2u2 

Y Y /L U  U / V V  V / V

 2 v2  

   / pV L

B. Normal shock relations – Conservation of mass, momentum and energy

 a y ( 0 )

1u1 

p  p/ ( pV2 )

  / L M = Mach Number

(0)

Since

 b y ( 0 )

 b 0 x (0)

 2 ,u2 are constant, the stream function is defined as :  b   2u2 y( 0 )

Knowing that across the shock

u11  u2  2 :

 b   a( 0 ) The stream function is continuous across the shock.

CNPq – Conselho Nacional de Desenvolvimento Científico e Tecnológico is the funding agency that provided the author with a scholarship inside the Science Without Borders program.

   u y (0)

C. Oblique Shock Relations

Ahead of the shock    VY ( 0 ) Hence, in dimensionless form:

 Y

The body surface is defined by

ψ=0 .

Tangential velocity is continuous i.e. .

u1 cos   u2 cos   v2 sin  ( 0 )

p2 2  1 (M12 sin 2   1) p1 (  1) (0) Conditions with substitution are conditions at ∞ - hence : p2 2  1 (M  2 sin 2   1) p ( 1) (0) and  2 (  1) (M  2 sin 2  1)    (  1) (M  2 sin 2   2) ( 0 ) D. Change Coordinate System From (X,Y) to (x,y). Let κ (x) denote the curvature of the bow shock.

a.

Detail for the analysis by Maslen

Since the stream value ψ = constant across the shock, taking Figure 4 as an example to go through the method, the stream value at point 2 is:

 2   2'  Y2' ( 0)

The shock curvature can be calculated by the formula (16) from M. M. Schulreich and D. Breitschwerdt’s article [5] :

Y 1 is known so is ψ 1 . ψ 2 ' is known so is Y 2 ' .

From Eq.(21), if

:



If

d2 r dz2 dr 2   1 ( dz)  

3/2

(0)

The mass conservation equation then becomes:   ( u)   (1   y) v  0 x y (0) in dimensionless form. Stream function in (x,y) is then such that :   (1  y) v x (0)

Since 2-2’ is a streamline then from isentropic relations and knowing that the flow is isentropic along any streamline i.e. s 2=s2 ' :

1 1 h2  (U 2 2  V2 2 )  h  V 2 2 2 (0) 1 2 1 1 1 V h2  V 2 (U 2 2  V2 2 )  V 2 h  V 2 2 2 2 2 i.e. (0) It can be demonstrated that:

1 2 1 2 2 1 2 2 V  V V2 ≪ V u2 2 2 2 (0)

 C  A 

1   AuA  CuC   yC  yA  2 (0)

Then 2 C  A  yC  yA   AuA   CuC ( 0 )

Hence, Eq.(24) can be approximated as : (0)

Giving an approximate value for

u2 : , (0)

E. Von Mises Transformation A Von Mises Transformation of governing equations is done to change the independent variables from (x,y) to (x, ψ ). Then, approximating for a thin shock layer gives:  p  u   x (0)

Turning equation (27) into a non dimensional form, gives :  p  u   x (0)

It can be approximated as :  p  u1 1   x (0)

This Eq.(29) gives the formula to calculate the pressure at point 2 :

p2  p1  u1 1 ( 2   1 ) ( 0 )

To find the physical coordinate Y, which corresponds to a particular value of ψ an integration of the definition of the stream function must be done : d  u dy x (0) then s d y  u 

(0)

This also locates the body coordinate where : s d yb   u ( 0 ) 0 From a point A to another point C , both at the y direction normal to the shock, eq. (33) can be approximated as:

IV.

PARABOLOIDAL SHOCK EXAMPLE

F.

Approximate equations for axisymmetric case The same procedure of the method can be done in the axisymmetric case, only requiring small modifications to some of the equations. The shock is defined as :

r 2  2Lz with

z  L z and r  L r shock equation in a non dimensional form is:

r 2  2z

The stream function from Eq. (20) becomes:

   V

r2 2 (0)

2 since    V L 

then, at the shock:

s 

rs2  zs 2 (0)

Equation (27) , at the shock in axisymmetric becomes :  p us s   rs ( 0 ) Turning Eq. (38) into non dimensional form :  p us s   rs ( 0 ) A change is also made in Eq.(30) to find the pressure: u p  ps  s s (   s ) rs (0) Hence, Eq.(32) turns into:  1 s d y  rs   u (0)

or Eq.(42) below can be used directly :

s

r1w  rs1w  (w 1)us  

d u

     1   2z     /  s   1  p ps  1        4 2z 1   (0)

(0)

For an axisymmetric case w=1, hence: s d r 2  rs2  2us  u 

Where  1  1   p 1 2     1  2 

(0)

Eq.(42) can be used as well in a two dimensional case(w=0).



4 2 p  u   2     1  1  

G. Exact Solution The equation of the shock is rs  (2z) . Consider a perfect gas with γ =1.4 and M ∞=∞ . From this, at the shock, the exact solution can be found with the equation used in Maslen’s [1], using Rankine-Hugoniot conditions : 2  1  ps      1  2z 1  ( 0 )



s  

s

r   2z 2us   

1/2

d   u 

1/

(0)

1/2

(0)

1/2

(0)

The result of this exact solution example in axisymmetric fow is shown in Figure 8 at section VI of this paper.

 1   1( 0 ) 1/2

2z  us    1 2z  (0)

 s  z( 0 ) V.

PROGRAM DEVELOPED

The MATLAB code presented below is for the axisymmetric example with a shock shape assumed of r= √ 2 z at free stream Mach Number and γ =1.4 The code was developed in a way that makes it easy to be adjusted for different shock shape examples, because it only needs the sub-functions to be edited. %========================================================================== function maslen_axi %========================================================================== clear all;  close all;  clc; %========================================================================== % flow data %========================================================================== gamma= 1.4;             % Cp/Cv Mach=5*10^10;           % free stream Mach number  w=1;                    % w=1  : axisymmetric flow                         % w=0  : plane flow %========================================================================== % define number of solution lines, ns, and number of points on each line,np %==========================================================================

ns=100; np=200;  nsh=100;                % number of shock points %========================================================================== %initialise vectors ­ non dimensional variables %==========================================================================   u=zeros(ns,np);     % tangential velocity p=zeros(ns,np);     % pressure rho=zeros(ns,np);   % density psi=zeros(ns,np);   % stream function beta=zeros(1,ns);   % angle beta  sin2b=zeros(1,ns);  % squared sine of beta cos2b=zeros(1,ns);  % squared cosine of beta zs=zeros(1,nsh);    % shock shape z coordinate  rs=zeros(1,nsh);    % shock shape r coordinate  k=zeros(1,ns);      % shock curvature y=zeros(ns,np);     %  y2=zeros(ns,np);    % body shape y coordinate  x2=zeros(ns,np);    % body shape x coordinate dn=zeros(ns,np);    %  beta1=zeros(ns,np); % angle beta at point 1 p1=zeros(ns,np);    % pressure at point 1 rho1=zeros(ns,np);  % density at point 1   %========================================================================== % compute shock points for plotting %========================================================================== X0=0.01; Xf=10;   for is=1:nsh     zs(is)=X0+(Xf­X0)*(is­1)/(nsh­1);     rs(is)=shock_eq(zs(is),1); end %========================================================================== % compute solution along each normal %========================================================================== for i=1:ns     %======================================================================     % conditions immediately behind the shock     %======================================================================     [beta(i)]=slope(zs(i));     cos2b(i)=(cos(beta(i)))^2;     sin2b(i)=(sin(beta(i)))^2;     u(i,1)=sqrt(cos2b(i));     r1=(gamma+1)*(Mach^2)*sin2b(i);     r2=(((gamma­1)*(Mach^2)*sin2b(i))+2);     rho(i,1)=r1/r2;     sp=(2*gamma)*((Mach^2)*sin2b(i)­1)/(gamma+1);     p(i,1)=(1+sp)/(gamma*Mach^2);     psi(i,1)=(rs(i)^(1+w))/(1+w);   % stream function normalized     k(i)=curvat(zs(i));     y(i,1)=0;     x2(i,1)=zs(i);

    y2(i,1)=rs(i);     %======================================================================     % compute solution at each other point     %======================================================================     for j=2:np         if j==np             psi(i,j)=0;         else         psi(i,j)=psi(i,1)­(j­1)*(psi(i,1)/(np­1));         end         p(i,j)=p(i,j­1)+ (u(i,1)*k(i)*(psi(i,j)­psi(i,j­1)))/((rs(i))^w);                x1=shock_eq(psi(i,j),0);         beta1(i,1)=beta(i);         [beta1(i,j)]=slope(x1);         s2b=(sin(beta1(i,j)))^2;         sp1=(2*gamma)*((Mach^2)*s2b­1)/(gamma+1);         p1(i,1)=p(i,1);         p1(i,j)=(1+sp1)/(gamma*Mach^2);         rho1(i,1)=rho(i,1);         rho1(i,j)=((gamma+1)*(Mach^2)*s2b)/(((gamma­1)*(Mach^2)*s2b)+2);         %h1=2*gamma*p1/((gamma­1)*rho1);         %u1=sqrt(c2b);         h0=(2/((gamma­1)*Mach^2))+1;  % total enthalpy         rho(i,j)= rho1(i,j)*(p(i,j)/p1(i,j))^(1/gamma);         hp=(2*gamma/(gamma­1))*(p(i,j)/rho(i,j));         u(i,j)=((h0­hp))^(1/2);                          dn(i,j)=(2*(psi(i,j)­psi(i,j­1)))/(rho(i,j)*u(i,j)+rho(i,j­1)*u(i,j­1));        y(i,j)=y(i,j­1)­dn(i,j);        y2(i,j)=((rs(i))^2­(2*u(i,1)*y(i,j)))^(0.5);        x2(i,j)=zs(i)+(rs(i)­y2(i,j))*tan(beta(i));             end end   %========================================================================== % plot results %========================================================================== if w==1 %========================================================================== % axisymmetric flow %==========================================================================  figure(1)  hold on  contourf(x2,y2,p,100,'EdgeColor','none')  plot(zs,rs,'b',x2(:,np),y2(:,np),'r');  grid on  xlabel(' Z coordinate ','FontSize',16);  ylabel(' R coordinate ','FontSize',16 );  legend('pressure profile','shock shape','body'); %========================================================================== % pressure plot over body %==========================================================================

 figure(2)  plot(zs,p(:,np),'b');  set(gca,'FontSize',16)  xlabel(' Z AXIS ' , 'FontSize',16);  ylabel(' Pressure at the body','FontSize',16 )   end   if w==0 %========================================================================== % two­dimensional(plane) flow %==========================================================================  figure(1)  plot(zs,rs,'b',x2(:,np),y2(:,np),'r');  grid on  xlabel(' X coordinate ');  ylabel(' Y coordinate ' );  legend('shock wave','body'); %========================================================================== % pressure plot over body %==========================================================================  figure(2)  plot(zs,p(:,np),'b');  xlabel(' X AXIS ' );  ylabel(' Pressure at the body' );   end %========================================================================== %SUB­FUNCTIONS function Y=shock_eq(X,imark) %========================================================================== % shock equation Y^2=X %========================================================================== if imark==1     Y=sqrt(2*X); end if imark==­1     Y=(X^2)/2; end if imark==0     Y=X; end return %========================================================================== % tan(beta) = dy/dx %========================================================================== function [beta]=slope(X)   beta=atan((2*X)^(­0.5));   return %========================================================================== % shock curvature k %========================================================================== function k=curvat(X)

dydx=(2*X)^(­0.5);          % first derivative of shock equation d2ydx2=(­1/((2*X)^(1.5)));  % second derivative of shock equation k=(abs(d2ydx2))/((1+(dydx)^2)^(1.5)); %Rs=k^­1; return

VI.

RESULTS

The results are shown for several tests made to solve the flow and body shape of different shock wave shapes stated on the code in both two dimensional and axisymmetric flow. The comparison of results is made, validating the accuracy of this Maslen’s approximate approach.

Figure 6.

Approximate solution for paraboloidal shock of

equation

r= √ 2 z

at

M ∞=∞

using the

MATLAB ® program in axisymmetric case. The dimensionless pressure profile is shown with red the high and blue the low pressure.

Figure 4.

Body associated with shock of equation

y=√ x

(

M ∞=5.8, γ =1.4 ¿

. Example

solved and plotted with the MATLAB ® program in a 2D case.

Figure 7. Figure 5. equation

Pressure on the body associated with the shock

y=√ x

( M ∞=5.8, γ =1.4)

.

Exact solution for the paraboloidal shock of equation r= 2 z at infinite free stream Mach number in axisymmetric case. The dimensionless pressure profile is shown

√

with red the high and blue the low pressure.

Figure 8.

MATLAB ® program in axisymmetric case.

Pressure on the body associated with the shock of

equation equation

r= √ 2 z

at

M ∞=∞

using the Figure 9. Comparison between the numerical and analytical solution of a bow shock wave generated by a hypersonic flow (M∞ = 500) hitting the

solid unit sphere. A perfect monoatomic gas is assumed (γ = 5/3). The analytical body prediction is represented by the red line. The bow shock wave is parametrized by rˆ = (2.620 · zˆ)0.471 with w= 1. The

VII.

CONCLUSIONS

The MATLAB ® code was validated by comparison with the exact solution for the case of infinite free stream Mach number. Maslen’s method is currently an indirect method and for future applications it should be extended to create a direct method in which the body shape is prescribed and the problem is iteratively solved by comparing a computed body shape with the defined shape, adjusting the shock There is also a uniformly valid extension that wasn’t used in this project and should be pursued to naturally provide the REFERENCES [1]

[2]

[3]

Maslen S. H. Inviscid hypersonic flow past smooth symmetric bodies (AIAA). AIAA Journal. 1964 [cited 23 August 2016];2(6):1055-1059. Available from: http://dx.doi.org/10.2514/3.2494 Reaction Engines Ltd - Space Access: SKYLON [Internet]. Reactionengines.co.uk. 2016 [cited 25 August 2016]. Available from: http://www.reactionengines.co.uk/space_skylon.html Our vehicles and spaceships [Internet]. Virgin Galactic. 2016 [cited 25 August 2016]. Available from: http://www.virgingalactic.com/human-spaceflight/our-vehicles/

(dimensionless) density profile shown is colour-coded, with red the high and blue the low density. Example taken from [5].

solution in the region behind the shock that is near the axis and the stagnation point. location appropriately and repeating this process to convergence.

ACKNOWLEDGMENTS

B. Schlottfeldt Maia thanks Kenneth Morgan for taking part as a supervisor in this research project and providing assistance from its start to finish.

[4]

[5]

[6]

Anderson J. D. Hypersonic and high-temperature gas dynamics. Reston, Va.: American Institute of Aeronautics and Astronautics; 2006. (pp. 36, 167-172) Schulreich, M. M., Breitschwerdt, D. Astrophysical bow shocks: an analytical solution for the hypersonic blunt body problem in the intergalactic medium. Astronomy & Astrophysics. 2011;531:A13. Sokhanvar F. , Khoshnevis A. B. Approximate Method of Calculation of Inviscid Hypersonic Flow. International Journal of Mechanical, Aerospace, Industrial, Mechatronic and Manufacturing Engineering. 2009;3(10):1272-1277.