Control-oriented modeling and real-time simulation method for a dual-mode scramjet combustor

Accepted Manuscript Control-oriented modeling and real-time simulation method for a dual-mode scramjet combustor Jicheng Ma, Juntao Chang, Junlong Zhang, Wen Bao, Daren Yu PII:

S0094-5765(18)31018-X

DOI:

10.1016/j.actaastro.2018.10.002

Reference:

AA 7123

To appear in:

Acta Astronautica

Received Date: 15 June 2018 Revised Date:

18 August 2018

Accepted Date: 1 October 2018

Please cite this article as: J. Ma, J. Chang, J. Zhang, W. Bao, D. Yu, Control-oriented modeling and real-time simulation method for a dual-mode scramjet combustor, Acta Astronautica (2018), doi: https:// doi.org/10.1016/j.actaastro.2018.10.002. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT 1 2 3 4

Control-Oriented Modeling and Real-Time Simulation Method for a Dual-Mode Scramjet Combustor Jicheng Ma1, Juntao Chang2*, Junlong Zhang3, Wen Bao4, Daren Yu5 Harbin Institute of Technology, 150001 Heilongjiang, People’s Republic of China Abstract: The control-oriented modeling and real-time simulation method for a dual-mode scramjet

6

combustor has been conducted in this paper. The 1-D unsteady model coupled with isolator shock train

7

model and oblique shock wave modification can treat variable area, fuel addition, combustion heat

8

release, variable specific heat, inflow air vitiation component, wall friction and mixing efficiency.

9

Combustion is simplified into mass and heat addition process. So this paper discusses heat release

10

distribution obtaining method which can be applied in scramjet combustor 1-D flowfield simulation in

11

detail. It is fit for certain scramjet combustor configuration and acceptable over a wide range of inflow

12

conditions. By mixed programming of C and Matlab, the code of model solving is capsuled into

13

S-Function which can be used in SIMULINK. The control-oriented model is built on SIMULINK

14

platform. When the model runs in a single-core 3.6GHz Intel Core i7-4790 processor, computation

15

time-consuming is found to require 5-10ms within a control period. The computing speed indicates

16

very promising, because model is compact enough to run in real time, and it can be used as an

17

embedded model in control system research. This model can also be used to analyze unsteady and

18

steady flowfield characteristic. The safety boundary codetermined by unstart and over-temperature is

19

conducted using this model.

SC

M AN U

TE D

20

RI PT

5

Keywords: Control-oriented modeling; Real-time simulation; One-dimensional model; Dual-mode

22

combustor; Heat release distribution

24 25 26 27 28 29

AC C

23

EP

21

1 PhD Candidate, School of Energy Science and Engineering, [email protected] 2 Professor, Academy of Fundamental and Interdisciplinary Sciences, [email protected]. 3 PhD Candidate, School of Energy Science and Engineering, [email protected]. 4 Professor, School of Energy Science and Engineering, [email protected]. 5 Professor, School of Energy Science and Engineering, [email protected]. 1

ACCEPTED MANUSCRIPT Nomenclature

30 =

cross-sectional area

x

=

axes coordinate

Cw

=

wetted perimeter

ρ

=

density

Et

=

total energy

=

combustion efficiency

f

=

friction coefficient

ϕ

=

fuel equivalence ratio

H

=

isolator height

θ

=

momentum thickness

Hs

=

total enthalpy

Ma =

Mach number

Subscripts

=

fuel addition ratio

comb =

combustion

p

=

static pressure

i

=

mesh serial number

Re

=

Reynold number

s

=

t

=

time

w

=

SC

RI PT

A

strut

M AN U

wall

Ttotal =

total temperature

u

=

velocity

=

fuel injection velocity along x axis

Supscripts n

=

time step serial number

1. Introduction

31

Scramjet is regarded as an ideal propulsion device during hypersonic flight within the atmosphere,

33

which owns many advantages such as long range, high specific impulse, high Mach cruise and direct

34

access to Earth orbit [1],[2]. Due to broad application prospects in space transportation and weapon

35

equipment, it aroused much interest in many countries. From starting acceleration, cruise to terminal

36

guidance, scramjet works over a broad range of hypersonic flight conditions. Hence, control system

37

must rapidly adjust working condition to match complex flight task. For hardware-in-loop simulation

38

system [3], scramjet combustor working process is replaced by control-oriented model, and working

39

states estimation does not rely on real experiment measurement directly. Correspondingly, combustor

40

model used for control algorithm evaluation is needed to estimate quickly performances such as thrust

41

to simulate the situation where flight Mach number and angle of attack change or any perturbations to

42

the design trajectory occurs. Control-oriented model of combustor requires a reduced-order model

43

with about 10ms of computational time. Based on above mentioned reasons, control-oriented

44

modeling and real-time simulation method for a dual-mode scramjet combustor are needed to develop

45

and estimate control system.

AC C

EP

TE D

32

46

A control-oriented model requires advantages both relatively small amount of computational time

47

and relatively accurate solution for the combustor flow. Multi-dimensional models can provide

48

detailed flowfield structure in scramjet combustor. The kinetic model and chemistry, which are the 2

ACCEPTED MANUSCRIPT most typical for the combustion mode in the engine, were applied in numerical simulation and show

50

satisfactory result [4]. Huang [5]-[9] developed many meaningful investigations about shock wave

51

transition, combustion and ram-scram transition mechanism in scramjet. They have instructive

52

significances to the scramjet numerical simulation. But computing burden of multi-dimensional

53

models is huge so that they cannot run in real time to meet the requirements of control design and

54

evaluation. A large number of table look-ups were also used to reduce execution time [10],[11]. It

55

means that run conditions were restricted between the bounds of pretabulated cases. Solving model

56

was reduced into an interpolation or a regression scheme. When pretabulated cases are sufficient, it

57

may have an acceptable precision for interpolation. But there exist uncertainties for extrapolation.

58

Since the conservation laws are not solved directly, it’s in the absence of physics content and accuracy.

59

Furthermore, there is phenomenon called mode transition accompanied with sudden change such as

60

thrust during working process of combustor. When fuel equivalence ratio changes from small to large

61

or from large to small, there exists different mode transition fuel equivalence ratio. Meanwhile angle

62

of attack in same flight Mach may arouse mode transition in certain case. A great deal of simulation

63

cases are needed if mode transition wants to be captured exactly. In addition, control research not only

64

needs thrust, but also pays more attention to changes of unsteady flowfield. As a consequence,

65

control-oriented model relies heavily on zero-or-one dimensional simulation. In previous works,

66

Lumped Parameter Method were applied to analysis overall performance of engine [12]-[14]. This

67

simple method can realize fast simulation. But it is possible that control system obtains unreliable

68

information so that scramjet combustor cannot fulfill a given control task. To solve this problem, a

69

method of distributed parameter control for dual-model scramjet combustor was proposed [15],[16].

70

One-dimensional (1-D) model can reflect distributed parameter characteristics of scramjet combustor

71

flowfield. It has been developed to support hypersonic vehicle design studies, flight dynamics, control

72

system design and analysis because of its moderate precision [17],[18].

EP

TE D

M AN U

SC

RI PT

49

At present, there are two kinds of 1-D scramjet combustor model, which is based on Partial

74

Differential Equations (PDEs) or Ordinary Differential Equations (ODEs). 1-D scramjet combustor

75

model based on PDEs was presented by Bussing [19] in 1980’s. The Navier-Stokes equations were

76

reduced to a set of one dimensional unsteady inviscid flow equations, namely Euler equations. The

77

governing equations were solved by MacCormack method. After that, this method was improved to

78

handle the effects of heat release source terms addition, area variation, friction, mass addition, heat

79

transfer to combustor wall, and variable specific heat, and it was used to analyze unsteady process and

80

supersonic combustor performance [20],[21]. Evaporation of liquid fuels and the oblique shock train in

81

the isolator were taken consideration into 1-D unsteady model so that the model could simulate various

82

fuels and show dual-mode characteristic [22]. Fewer models based on PDEs were used in scramjet

83

combustor because solving model needs complex finite-difference method such as MacCormack or the

84

advection upstream splitting method (AUSM). 1-D scramjet combustor model based on ODEs can be

AC C

73

3

ACCEPTED MANUSCRIPT referred in literature [23]. This method has been widely applied in performance analysis, combustor

86

initial design and control research. Previous models made utmost efforts to develop the accuracy of

87

simulation by taking consideration of the effects of fuel injection, mixing, chemical production rates,

88

heat transfer, and viscous losses [24]-[26]. The computation time-consuming of engine flowfield is on

89

the order of 1-2s. This model has been applied to compute the sensitivity of the thrust to the vehicle

90

angle of attack, as well as the poles and zeros of the transfer functions that relate control inputs of fuel

91

setting and elevator [27],[28]. In addition, many other researchers made considerable efforts to develop

92

fast 1-D simulations that can run in short time periods [29],[30] and evaluate dual-mode characteristic

93

[31]. However, solving ODEs process exists a mathematical singularity in the equations when the flow

94

velocity approaches the speed of sound [32],[33]. This singularity problem becomes significantly

95

important when considering subsonic combustion cases. Furthermore, dealing with combustion

96

chemistry in 1-D model is a key issue both in 1-D model based on PDEs and ODEs. Firstly, there exists

97

stiffness problem during solving chemical governing equations due to time scales mismatching

98

between chemistry and flowing. Thus, the solution procedure is computationally expensive. Meanwhile

99

the governing equations of fluid flow with detailed chemistry may generate a highly nonlinear problem

100

during integrating these equations. Secondly, even though the detailed chemistry mechanism is

101

accounted, the accuracy of computation would not change substantially. Lastly, if combustion is

102

simplified into combination of heat release and mass addition, there is a problem in heat release

103

distribution law and its universality over a wide range of conditions. It has apparent impact on

104

computation accuracy of flowfield. Moreover, heat release distribution has not yet coincident empirical

105

modification method and relies heavily on experimental data.

TE D

M AN U

SC

RI PT

85

In the field of scramjet combustor, related researchers have also developed many reduced-order

107

models, but the computation time-consuming is mostly measured on the order of seconds. It is obvious

108

that these models cannot accomplish real-time optimization of scramjet combustor control system. By

109

pre-setting controller parameters to tolerate all kinds of uncertainty and perturbations, performance of

110

scramjet combustor is not optimal. As a whole, control research is in sore need of 1-D control-oriented

111

model which has advantages both relatively small amount of computational time and relatively

112

accurate solution for the combustor flow. Hence, this paper proposes control-oriented modeling and

113

real-time simulation method for a dual-mode scramjet combustor. In this paper, a 1-D unsteady model

114

based on PDEs is discussed in detail. This method can conduct transient analysis without assumption of

115

the shock location. Flowfield characteristics distributions are calculated at each time step of a

116

simulation, instead of being stored in look-up tables. Then flowfield characteristics in this time step are

117

given to next time step initial state of the combustor. The model coupled with isolator shock train

118

model and oblique shock wave modification can treat variable area, fuel addition, variable specific heat,

119

inflow air vitiation, wall friction, combustion and mixing efficiency. To shorten computation

120

time-consuming, combustion is simplified into mass and heat addition process without the process of

AC C

EP

106

4

ACCEPTED MANUSCRIPT 121

spraying and mixing. Meanwhile, mixed programming of C and Matlab is applied. Using these

122

methods mentioned above, this model can meet the requirement of real-time simulation. Additionally,

123

model presented in this paper can also been used to analyze unsteady and steady flowfield

124

characteristic. The safety boundary decided by unstart and over-temperature is conducted in the end of

125

this paper.

127

RI PT

2. One-dimensional modeling method

126

2.1. Governing equations and numerical solution method

Flow in a duct having a slowly varying cross-section where the duct height is small compared to the

129

radius of curvature of the axis of the duct is termed a quasi-one-dimensional flow [34]. A truly 1-D

130

unsteady flow is one where the flow properties depend only on space coordinate x and time t. The

131

pathlines of the fluid particles are straight lines, and the flow properties are uniform across all surfaces

132

perpendicular to the main flow direction. The governing equations for unsteady 1-D flow are derived in

133

this section. If the effects of work and body forces are neglected, the flow equations are inviscid

134

adiabatic, so an equivalent wall friction model has been incorporated into the 1-D model to give more

135

realistic drag estimates. In that case, the flow area is only a function of x and t. Based on unsteady 1-D

136

Euler equations, the model accounted for variable area, fuel addition, combustion heat release, variable

137

specific heat, inflow air vitiation component, wall friction and mixing efficiency. Mixture gas

138

parameters can be obtained by mass-weighted average. Constant-pressure specific heat Cp and enthalpy

139

of single component are function of temperature that can be referred at NIST. Governing equations can

140

be expressed as follows:

141 142

=

In equation (1)：

144

,

,

=

We want to evolve

/

!"# !

,

+

!$ !

；

−

&'(

,

= +

(1)

；

) ∙ +, +

!"# !

,- ∙

from time 0 = 0 / to the new value

/12

∆0 = 0 /12 −0 / . The splitting can be derived as follows [35]:

∙

!"# !

. ；

at 0 = 0 /12 in a time step

148

AC C

145

=

,

+

EP

143

TE D

M AN U

SC

128

149

For solving the advection (2), AUSM flux vector splitting is applied for the approximation of the

146 147

150 151

152

+

!

5, 0 / =

!

= 8 /12

=4 /

6 ⇒ 8 /12

9⇒

(2)

/12

(3)

convective flux function. Flux in the duct with varied area can be expressed as follows: :12⁄

=

:12⁄

<

+

=

:12⁄

+ >:1 <

0 0

1 = + >:12 < :

0 0

=

:12

(4)

For source terms that are independent of time or involve spatial derivatives, one can use the Euler 5

ACCEPTED MANUSCRIPT 153

method resulting from evaluating the integral at the old time. The Euler method is explicit and

154

first-order accurate. For solving equation (3), it can be expressed as follows: /12 = 8 /12 + ∆0 ∙ 8 /12

155 156 157

Where ∆0 = 0 /12 −0 / is the time step and

/

approximately equals

2.2. Normal shock wave treatment and isolator train model

(5)

0/ .

In 1-D model, there are substantial limitations, such as no oblique shock wave but normal shock

159

wave, which is unreasonable. Therefore, proper modification is necessary. After model is solved, if

160

back-to-front static pressure ratio along axial detected is greater than a specified value, it is thought that

161

there exists normal shock wave. Then the calculation results are modified according to the law of

162

quadratic function as shown in Fig. 1. By this modification method, the model can simulate transition

163

between ram-mode and scram-mode operation which previous models did not include.

p/MPa

0.3

0.25

0.2 0.25

M AN U

Calculation Modification

0.35

164

SC

RI PT

158

0.3

0.35

0.4

0.45

x/m

Fig. 1. Pressure modification schematic at normal shock wave

TE D

165 166

As a performance adjustable bumper and engine regulating actuator, isolator can expand scramjet

167

combustor operating domain. Isolator shock train has obvious and complex non-1-D characteristic. To

168

catch this special phenomenon accurately in 1-D model, Billig’s empirical method [36] is used in this

169

work. Firstly, selecting isolator exit pressure as backpressure

170

A BCD( E2 IFG@H

171

Pressure distribution is modified from shock train leading edge point to isolator exit. In the formula, L

172

is shock train length in the isolator, if it exceeds isolator length, unstart will generate and computation

173

is terminated; Mai is isolator inlet Mach number; Reθ is Reynold number which based on inlet

174

momentum thickness, generally it is 10000; H is isolator height; θ is inlet boundary layer momentum

175

thickness; θ/H is equal to 0.02 generally[1],[37]. It is just a method in 1-D simulation when

176

experimental data is lacked. Precision will fit for control research as long as they are defined as

177

relatively reasonable values..

178

2.3. Mixed programming method of C and Matlab

= M50

OP OD

− 1. + 170

OP OD

and then equation

− 1. S is applied to solve isolator shock train length.

AC C

√KL

EP

(

@,

179

S-function is system function in Simulink which is used to describe module macro function. When

180

S-function modules provided by Simulink default don’t meet the requirements, users can make

181

modules by themselves to realize user-defined algorithm and desired action. S-function can be 6

ACCEPTED MANUSCRIPT 182

programmed in several kinds of languages, including M, FORTRAN and C. In reference [38], a

183

separate FORTRAN code interfaces with MATLAB as a MEX file in the research of reacting flows

184

with detailed chemistry and transport. C language has flexibility in programming, enormous computation function and a faster running

186

speed. Similarly, S-function can also be programmed into C Mex file using C language. Compared

187

with running interpreted S-function which is programmed using M language, MATLAB interpreter is

188

not called repeatedly during simulation, so the method has obvious advantages in runtime and

189

calculation efficiency. In addition, S-function can be applied to hardware-in-loop simulation platform

190

conveniently. Hence, this paper presents that the model is solved using C programming language, and

191

then the code is encapsulated into S-function.

193

According to program structure sequence, C Mex S-function file is composed of such sections as

SC

192

RI PT

185

follows [39]:

(1) The S-function begins with #define statements for the S-function name and level;

195

(2) Header files include parts of C Mex S-function core data structure and other library files.

196

Moreover it also includes other header files as needed. For example, math.h is included if code goes on

197

arithmetical operation utilizing C standard library; fixedpoint.h is included if code uses built-in

198

fixedpoint data type in SIMULINK. Every C Mex S-function has its own simstruct data type plant to

199

ensure each S-function internal data independent.

M AN U

194

(3) The definition of parameter dialog visiting macros function: it includes definition of macros

201

variable TRUE/FALSE. By these macros, C Mex S-function module can obtain the value of each

202

control in parameter dialog conveniently, so that it is applied in each submethod.

203 204

TE D

200

(4) After that, in C Mex S-function each submethod is defined, such as Initialize, InitializeSampleTimes, updating state variables, Output, mdlRTW, terminate and so on. (5) At the end, S-function must include necessary source and header files according using condition,

206

so that it connects SIMULINK or SIMULINK Coder products. If this section is ignored, the whole

207

compilation will fail.

AC C

EP

205

208

After code structure is built as listed above, through a special compilation process, code can be

209

compiled into machine language files, such as mexw32 or mexw64 file. Only after code is compiled

210

into executable file can it be called by SIMULINK. Code using compiled language directly instead of

211

compiling repeatedly can ensure high efficiency operation. This is also the reason why users are

212

recommended the use of C Mex S-function when fast simulation is desired.

213

3. Heat release distribution obtaining and real-time simulation

214

3.1. Heat release distribution obtaining

215

3.1.1. Experimental setup descriptions

216

To validate the accuracy of scramjet combustor 1-D control-oriented model numerical results, parts

7

ACCEPTED MANUSCRIPT 218

Technology in China [40],[42]. In the scramjet combustor, there are a constant width of 100mm and

219

single side divergence lower wall. The total length of the combustor is 1.620m with a total divergent

220

ratio of 2.5. Some pressure sensors are distributed on the center line of side wall. Schematic

221

configuration of combustor is shown in Fig. 2. The detailed geometry of strut can be found in reference

222

[43]. Strut fuel injector employed as the first stage fuel injector is located at x=0.36m. The second stage

223

fuel injected from wall fuel injector is located at x=1.08m. Hydrocarbon fuel at room temperature

224

(about 300K) is injected in the direction normal to the airflow direction, with fuel total equivalence

225

ratio ranging from 0.1 to 1. The hydrocarbon fuel used in experiments is China No.3 aviation kerosene.

226

In computation, kerosene surrogate molecular formula is C10H22, whose heating value of combustion is

227

42.5MJ/Kg. PJ ignitor

100

M AN U

228

Pressure sensor

84

60

Wall fuel injector

100

Strut fuel injector O2

SC

RI PT

of experiments has been conducted at the supersonic combustion research facility of Harbin Institute of

40

217

260

400

229 230

100

320

420

100

Fig. 2. Schematic configuration of combustor (mm) 3.1.2. Heat release distribution obtaining method

Modifications and simplifications must be made for the mathematical model to allow the application

232

of a given tasks [44]. When detailed chemistry mechanism is accounted, accuracy of code can be

233

developed in a degree, but computation time-consuming is huge. In the model presented in this paper,

234

combustion process is simplified into combination of mass and heat addition. In calculation, heat

235

addition process can be described as heat release distribution function. By tailoring the heat release

236

profile, combustion operation can be modified to meet desired requirements. In 1-D model presented

237

by Jiang [45], mixing between air and fuel was taken into account, and fuel addition ratio along axial

238

direction was described using Rayleigh probability density method. Heat release was associated with

239

fuel mass distribution, so they had the same distribution law. Zhang [46] summarized existing

240

experimental results where OH spontaneous radiation plane images were observed in combustor when

241

fuel equivalence ratio was low, presented 1-D heat release distribution peicewise linear fitting model.

242

But this model can only be used in specific combustor configuration and fuel component. For

243

multi-stage injection and hydrocarbon fuel, some corresponding experiments were needed to analyze

244

image characteristic and adjust model parameters. Furthermore, this model is only fit for high speed

245

inflow condition and ramjet mode. In addition, parabolic fuel distribution and exponent mixing

246

distribution were also adopted in some literatures [47],[48].

AC C

EP

TE D

231

247

The above mentioned studies about the heat release distribution function mainly rely on existing

248

experimental data or specific parameter setup in their own model. However, it is not known if the heat 8

ACCEPTED MANUSCRIPT release function is fit for different scramjet combustor configuration, different fuel, different injection

250

form and multiple inflow conditions. For the moment, researchers lack unified, normative description

251

method and measurement means with regard to exploration and description of heat release distribution.

252

Hence, this paper presents a method that can be applied in scramjet combustor 1-D flowfield

253

simulation. Using this method and utilizing several existing experimental data, the heat release

254

distribution which is fit for certain scramjet combustor configuration and is acceptable over a wide

255

range of inflow conditions can be found. Heat release distribution modification and optimization flow

256

chart is shown in Fig. 3. Detailed procedures are listed from step 1 to 5 as follows:

RI PT

249

Step 1:

258

In the supersonic combustion research facility introduced in Section 3.1.1, cold flow experiments

259

were conducted for the scramjet combustor, and isolator inflow conditions, including total temperature,

260

static pressure, Mach number, and inflow air vitiation component, were measured. They are listed in

261

case 1, case2, case 3 and case 4 as shown in Table 1. In computation, boundary conditions are same as

262

Table 1. Wall is treated as adiabatic surface. Then density, constant-pressure specific heat, specific heat

263

ratio and other corresponding parameters can be obtained.

TE D

M AN U

SC

257

265

AC C

EP

264

Fig. 3. Heat release distribution modification and optimization flow chart

266

Table 1 Isolator incoming flow conditions in experiments Mole fractions ratio O ∶ N ∶ H O ∶ CO

Y + Y, allocation proportion

FZ[K

I

Case

Ma

Ttotal (K)

p (MPa)

1

2.0

965

0.153

1.000:3.190:0.333:0.238

30%+70%

√\1757.5

2

2.5

1027

0.089

1.000:3.143:0.381:0.238

48%+52%

1837.5

3

2.7

1480

0.059

1.000:2.810:0.571:0.381

65%+35%

1780.0

4

2.7

1865

0.063

1.000:2.657:0.662:0.443

100%+0%

1711.3

9

(

ACCEPTED MANUSCRIPT 267

Step 2:

268

Combustion experiments in case 1 at ϕ=0.77, in case 2 at ϕ=1.0 and in case 3 at ϕ=1.0 separately

269

were conducted. Pressure distributions along flow direction were obtained. Fuel allocation proportion

270

injected from strut and side wall of combustor has been shown in Table 1. Step 3:

272

Two categories of heat release distribution function which are exponent and polynomial form are

273

selected as the examples in this method as follows:

274 ]

275

=

276

279

^

^

−5̅ + 25̅

(6)

8 2.2bcd ̅ e.f(g(h ∙ie.I(fj(k

2.lm2n ̅ 1c.d mnn

5̅ =

(7)

E opq

Arstu

(8)

^

is fuel consumption mass flow rate which reacts chemically;

is fuel total mass flow; x

SC

278

]

Where

=

is combustor position along axial; 5vwx is fuel injection position; Lcomb is the length of combustion section.

M AN U

277

]

RI PT

271

Step 4:

281

Equation (6) is applied to 1-D model firstly. If calculation results appear relatively large deviations

282

comparing with experimental results, Lcomb should be modified and optimized until achieving the

283

satisfied results. Otherwise return to step 3 and select equation (7). In this paper, when upstream fuel

284

injection Lcomb is between 0.3m and 0.5m, downstream fuel injection Lcomb is 0.3m in heat release

285

distributions laws based on Equation (6), normalized wall pressure distributions along with axial

286

direction can meet the accuracy requirements against with experimental results. It can be seen in Fig. 4,

287

Fig. 5 and Fig. 6 separately. Based on pressure values, Table 2 gives root-mean-square deviation

288

between experimental and numerical results with different Lcomb. Root-mean-square deviation between

289

experimental and numerical results can be expressed as follows:

TE D

280

∑D‰Š} yz = {

EP

ˆ

290 :

v ‹Œ• v w

/

(9)

: i Ži]v iw

is pressure value of simulation result in position i;

AC C

291

(

D D D O#ot~•€•osp EO‚ƒ„‚…ot‚p• .†O‚ƒ„‚…ot‚p• ‡

is pressure value of

292

experiment in position i; n is total number of pressure sensors. As we can see from Table 2, in heat

293

release distribution law based on Equation (6), upstream fuel injection Lcomb equals to 0.4m is the most

294

appropriate. 4

L

comb

=0.30m

L

3.5

comb

=0.40m

L

comb

=0.50m

φ=0.77experiment

5.5

area

5

70%

30%

Lcomb=0.30m

Lcomb=0.40m

Lcomb=0.50m

48%

φ=1.0experiment

area

52%

4.5

2

3.5

area/m2

pressure ratio

4

2.5

area/m2

pressure ratio

3

3 2.5 2

1.5

1.5

1

0.01

0.5 0

0.004 1.6

0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.01

1 0.5 0

10

0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

ACCEPTED MANUSCRIPT Fig. 4. Experimental and numerical normalized wall pressure distributions in case 1 at ϕ=0.77 7

L

=0.30m

comb

L

6

comb

=0.40m

Fig. 5. Experimental and numerical normalized wall pressure distributions in case 2 at ϕ=1.0 L

comb

=0.50m

φ=1.0experiment

area

35%

65%

4

area/m2

pressure ratio

5

3 2 1 0 0

0.2

0.4

0.6

0.8 x/m

1

1.2

RI PT

0.01 0.004 1.6

1.4

Fig. 6. Experimental and numerical normalized wall pressure distributions in case 3 at ϕ=1.0

Table 2 Root-mean-square deviation between experimental and numerical results with different Lcomb Case Fuel equivalent ratio Lcomb=0.3m Lcomb=0.4m Lcomb=0.5m 1 0.77 12.8% 6.2% 6.4% 2

1.0

21.7%

3

1.0

13.9%

SC

295

10.1%

11.0%

8.7%

9.9%

Step 5:

297

Combustion experiments in case 3, 4 at different fuel equivalence ratio and fuel allocation proportion

298

separately were carried out for this scramjet combustor geometry. They are used to validate the

299

accuracy of upstream Lcomb=0.4m in other fuel equivalent ratios, fuel allocation proportion and inflow

300

conditions. Numerical results are compared with experimental results as shown in Fig. 7-Fig. 9. If

301

Lcomb=0.4m cannot meet the needed accuracy, return to step 4 and select different Lcomb until satisfied

302

value, otherwise return to step 3 and select equation (7).

3

1

0.2

0.4

EP

2

0.6

0.8 x/m

1

1.2

φ=1.00

φ=1.00experiment

area

100%

5 4 3 2

1.4

0.01

1

0.004 1.6

0 0

0.01 0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

Fig. 8. Experimental and numerical normalized wall pressure distributions in case 4 at ϕ=0.74 and 1.0

φ=0.48experiment

φ=1.2 (40%+60%)

100% 40%

5

φ=1.2experiment

area

0% 60%

area/m2

4 pressure ratio

AC C

φ=0.74experiment

6

Fig. 7. Experimental and numerical normalized wall pressure distributions in case 3 at ϕ=0.74 φ=0.48 (100%+0%)

φ=0.74

7

area/m2

4

0 0

8

area

35%

65%

5

pressure ratio

φ=0.74experiment

area/m2

φ=0.74

pressure ratio

6

TE D

M AN U

296

3

2

1

0 0

0.01 0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

Fig. 9. Experimental and numerical normalized wall pressure distributions in case 3 at ϕ=0.48 and ϕ=1.2 303

Root-mean-square deviations between experimental and numerical result are listed in Table 3. It is 11

ACCEPTED MANUSCRIPT 304

obvious that the simulation results have a deviation within 10%, it can meet the need of control

305

research. Hence, numerical results of Lcomb=0.4 in other fuel equivalent ratios, fuel allocation

306

proportion and inflow conditions also meet the requirement of accuracy for 1-D simulation. With this,

307

heat release distribution is confirmed. Table 3 Root-mean-square deviation between experimental and numerical results in different case Case

Fuel equivalent ratio

Y + Y, fuel allocation

3

0.74

65%+35%

4

0.74

100%+0%

4

1.00

100%+0%

3

0.48

100%+0%

3

1.20

40%+60%

Root-mean-square deviation

RI PT

proportion

6.8%

5.9% 9.9% 8.9% 9.7%

SC

308

In addition, 3-D CFD simulation is necessary to show detailed flowfield structure and provides

310

comparison with 1-D simulation. In 3-D CFD simulation, accuracy of numerical scheme is

311

second-order upwind, which has two-order accuracy. The method of error estimation is from Smirnov

312

[49],[50]. The numerical calculation error estimates of 3-D simulation are shown in Table 4. The

313

40*50*800 grids are adopted as the calculation mesh at last. Computational fluid dynamics model

314

referred to reference [40]. Allowable error

Table 4 Numerical calculation error estimates of 3-D simulation Allowable Grid Time step Number of Accumulated number of time resolution size(µs) time steps error steps

5%

20*50*800

5%

40*50*800

5%

60*80*800

Reliability Rs=nmax/n

10

50000

1.330·10-4

1.413·105

2.8271

10

50000

2.363·10-5

4.479·106

89.583

10

50000

6.583·10-6

5.769·107

1153.9

TE D

315

M AN U

309

Fig. 10 shows static pressure contour in case 1 at ϕ=0.77. Meanwhile mass-weighted average of

317

static pressure is given in Fig. 11. 1-D simulation can have similar distribution along x-axis with 3-D

318

mass-weighted average of static pressure even though 1-D simulation cannot show detailed flowfield

319

structure. Nevertheless, 1-D simulation method is enough to control design and evaluation.

80000 Pa

1-D simulation

φ=0.77experiment

3-D CFD simulation

area

70%

30%

3 2.5

area/m2

Static pressure contour

4 3.5

pressure ratio

AC C

EP

316

2 1.5

400000 Pa

1 0.5 0

Fig. 10. Static pressure contour in case 1 at ϕ =0.77

0.01 0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

Fig. 11. Comparison of mass-weighted average of static pressure distribution, 1-D simulation and experiment

320

This 1-D model has been proved valid within the experimental bounds by above analysis. During the

321

scramjet vehicle control, the angle of attack will result in different total and static pressure incoming 12

ACCEPTED MANUSCRIPT 322

flow condition. In order to validate 1-D model in such flight conditions which is outside the

323

experimental bounds, a 3-D CFD is performed. Bear in mind, for the situation where the vehicle works

324

in some an angle of attack, incoming flow will still embody uniform equivalent conditions in 1-D

325

simulation. In the previous study, some inlet designs and 3-D CFD simulations had been conducted

326

according to the reference[41]. Flight conditions have been listed in Table 5, and the equivalent

327

incoming flow conditions of combustor has been listed in Table 6. Table 5 Flight conditions of scramjet vehicle Angle of Altitude(km) Static pressure(Pa) attack(degree) 2 21 4729

Flight Mach 5

RI PT

328

329

Table 6 Equivalent incoming flow conditions of combustor Incoming Mach

Static pressure(MPa)

2.15

0.1459

Static temperature(K) 678.4

By 3-D CFD simulation, static pressure contour at Y + Y, =0.3+0.5 is shown in Fig. 12.

M AN U

331

217.6

SC

330

Static temperature(K)

332

Meanwhile comparison of mass-weighted average of static pressure distribution and 1-D simulation is

333

given in Fig. 13. 1-D simulation can have similar distribution along x-axis with 3-D mass-weighted

334

average values. Hence, applicability of the 1-D simulation is checked for the situation where the

335

combustor works in angle of attack.

4

1-D simulation

3.5

3-D CFD simulation

area 63%

37%

Static pressure contour

100000 Pa

400000 Pa

area/m2

TE D

pressure ratio

3

2.5 2 1.5 1 0.5 0

Fig. 12. Static pressure contour at Y + Y, =0.3+0.5

0.01 0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

EP

Fig. 13. Comparison of mass-weighted average of static pressure distribution and 1-D simulation

It is necessary to validate 1-D model near unstart conditions. In case 1 at ϕ=0.77 and case 4 at ϕ=1.0,

337

the combustor works in such conditions. In case 1 at ϕ=0.77, comparison of 1-D, 3-D CFD simulation

338

and experiment have been conducted which are shown in Fig. 10 and Fig. 11. In case 4 at ϕ=1.0,

339

comparison of 1-D and experiment is shown in Fig. 8. By above researches, 1-D simulation shows

340

relatively satisfactory results near unstart conditions.

AC C

336

341

As a whole, numerical results obtained from 1-D model are in good agreement with experimental

342

and 3-D CFD results no matter the incoming flow is in or outside experimental bounds. It does

343

successfully reproduce the pressure rise trend along flow direction. So by above steps, heat release

344

distribution function which is fit for 1-D model can be obtained aiming at this scramjet combustor

345

configuration over a wide range of inflow conditions.

346

3.2. Real-time simulation method and validation

347

The goal of developing a real-time simulation is not to replace high-fidelity and 13

ACCEPTED MANUSCRIPT high-spatial-resolution CFD solutions or real-world experiments. Instead, we seek to match numerical

349

and experimental results to an adequate degree of accuracy for control design and evaluation. The

350

computation time-consuming for full flow path kept on the order of 10ms can satisfy this goal.

351

Handling methods are listed as followed: (1) Combustion is simplified into mass and heat addition

352

process, heat release distribution function used in 1-D model can be obtained from Section 3.1.2. (2)

353

The relationships between inflow air specific heat at constant pressure, enthalpy and temperature are

354

tabulated to look up at all times in later calculations. (3) 1-D unsteady Euler equation is solved using C

355

programming language, then it is compiled into machine language in MATLAB. Lastly it is

356

encapsulated into S-function referred in Section 2.3. On SIMULINK platform, 1-D real-time model is

357

built, where memory module can return the results computed by S-function as flowfield initial value at

358

the next time step. In this way, the model can realize the maximum reduction of the computation

359

time-consuming. Flow chart of real-time model is shown in Fig. 14. (4) In premise that accuracy meets

360

the requirement of control on the whole, and calculation grid number should be minimized.

TE D

M AN U

SC

RI PT

348

362

AC C

EP

361

Fig. 14. Scramjet combustor flowfield 1-D real-time simulation block diagram

363

Based on the order of the numerical scheme accuracy and grids, the different grids, including 50, 70

364

and 140, are used as the error estimates of numerical calculation. In this paper, accuracy of numerical

365

scheme is one-order upwind, which has one-order accuracy. CFL is selected to 0.7. According to the

366

method for accumulation of errors estimations in the reference [49],[50], several simulations have been

367

conducted. The numerical calculation error estimates of 1-D simulations are shown in Table 7. Taking

368

into account of computation accuracy and time-consuming, grid total numbers are terminally selected

369

to 70.

370

This paper develops 1-D model which can run in real time. On hardware-in-loop simulation platform

371

the model is available as a module. Bear in mind that although some parameters in the model do not 14

ACCEPTED MANUSCRIPT 372

attempt to match the combustor flowfield parameters exactly, they should be functionally accurate. For

373

example, if this model is used to fuel control research, it is sufficient since fuel control is not

374

instantaneous but rather a stepwise process. This model can predict an averaged flowfield which is

375

useful for fuel calculations. Allowable error 5%

Table 7 Numerical calculation error estimates of 1-D simulations Grid Time step Number of Accumulated Allowable number resolution size(µs) time steps error of time steps 50

2.68

2500

4.000·10-4 -4

5%

70

1.34

5000

2.041·10

5%

140

0.67

10000

5.102·10-5

Reliability Rs=nmax/n

1.563·104

6.25

RI PT

376

4

12.01

9.604·105

96.04

6.003·10

To verify the real-time characteristic of this model, a large number of calculations in multiple flow

378

conditions should be conducted. In this paper, fuel supply signal is given as a slope and inflow

379

conditions is case 1, 2, 3 and 4. Computation conditions have been shown in Table 1. The real-time

380

model which runs in a single-core 3.6GHz Intel Core i7-4790 processor can acquire the full combustor

381

1-D flowfield structure. In control system research, one control period is generally 20ms. In Table 8, it

382

can be seen that all computation time-consuming are less than 10ms which is enough to satisfy on-line

383

control, so the control-oriented model can realize real-time simulation of scramjet combustor flowfield

384

at different fuel equivalence ratio over a wide range of inflow conditions.

385

Table 8 Computation time-consuming at different fuel equivalence ratio in each case Fuel equivalence ratio Computation time-consuming (ms) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

M AN U

SC

377

Case 2 Case 3 Case 4

387

6

5

6

5

6

7

7

—

—

6

6

5

6

5

6

6

6

7

5

7

6

6

7

6

7

7

6

7

8

8

8

7

7

7

8

9

9

9

9

4. Flowfield characteristic analysis

EP

386

5

TE D

Case 1

1.0

4.1. Unsteady analysis

The process of this model being solved is based on time iteration, so it can show unsteady flowfield

389

characteristic. Steady flowfield at ϕ=0.2 are assigned to the whole flowfield as the initial values.

390

Numerical normalized wall pressure and Ma distributions at different time in case 1 at ϕ=0.77 are

391

shown in Fig. 15 and Fig. 16.

AC C

388

392

Overall, the example calculation of heat addition from the initial flowfield (at t =0) to an established

393

diverging duct flowfield containing a shock is considered. In the rear part of the combustor, flowfield

394

tends to be stable choke state quickly and generates high backpressure. The flow can be only adjusted

395

via a shock moving forward. In this way, the flowfield tends to stable at a slower rate. The flowfield

396

and shock movement requires a time of 8ms to reach equilibrium. Note the theoretical times only give

397

an order of magnitude estimate of the overall time scales. From the example, we could see that the flow

398

is not stationary and shocks travel upstream to adjust inflow conditions. The 1-D model is able to 15

ACCEPTED MANUSCRIPT 399

simulate transient behavior of choked flow. 4

t=2ms

t=4ms

3.5

t=6ms

t=8ms

t=10ms

area

t=2ms

t=4ms

2

70%

30%

t=6ms

t=8ms

30%

t=10ms

area

70%

Ma

2

area/m2

1.5

area/m2

pressure ratio

3 2.5

1

1.5 1 0.5

0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

0

Fig. 15. Numerical normalized wall pressure distributions at different time in case 1 at ϕ=0.77

0.2

0.4

RI PT

0 0

0.01

0.01

0.5

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

Fig. 16. Numerical Ma distributions at different time in case 1 at ϕ =0.77

Firstly, pressure field is analyzed. Fuel allocation proportion injected from the rear region is larger

401

than the front region in this case. Hence in the initial time, higher pressure peak firstly appears in the

402

rear region. In the same position Ma is equal to 0.4 that is far less than 1.0, so a strong thermal

403

chocking is inferred. In other word shock leads to pressure peak. Then the shock moves forward to

404

small area region, as a consequence of pressure lifting point moving forward which is located in front

405

of heat release region. Pressure peak value continue rising over time, simultaneously pressure lifting

406

point moves forward and is connected with the first stage heat release region. After t=6ms, flowfield

407

shows a slower growing trend. The backpressure at x=0.36m has risen less and maintain a high level

408

due to isolator. Ultimately, approximate constant pressure combustion is observed in the first stage heat

409

release region.

M AN U

SC

400

Secondly, from the perspective of Mach number, the second stage combustion happening in constant

411

area region generates shock waves and leads to subsonic region. Nevertheless, in the initial flowfield

412

establishing process due to the second stage combustion has not yet established enough backpressure,

413

strong thermal choking impact is not obvious. So flowfield in the first stage combustion region which

414

happens in diverging area, transforms between subsonic and supersonic. After that, since the steady

415

stage of the flow relies on the propagation of information from the choking location in the combustor

416

up to the upstream supersonic inflow. Hence, subsonic region caused by the second stage combustion is

417

enlarged and combines with the first stage combustion region. Finally, the subsonic region generates

418

from the location of the isolator entrance, and the whole combustor is in a subsonic combustion mode.

AC C

EP

TE D

410

419

Eventually, flowfield is steady after t=8ms. Pressure rises from x=0.2m in the isolator due to

420

pre-combustion shock train. Then combustion heat release in the diverging region is balanced out by

421

area change, so pressure platform and Ma platform occurs. When heat release process is over, flowfield

422

is regarded as the simple flowing process of the high temperature and high pressure gas. In second

423

stage combustion region, flowfield velocity rises and pressure declines by the heating process, and heat

424

release process accelerates the flow towards the sonic limit till supersonic in nozzle region.

425

4.2. Steady analysis

426

Steady pressure distribution characteristics at different fuel equivalence ratio are analyzed below. Fig.

16

ACCEPTED MANUSCRIPT 17 and Fig. 18 show numerical normalized wall pressure and Ma steady distributions at different fuel

428

equivalent ratios in case 2. Pressure along with the flow direction rises from vicinity of isolator

429

entrance gradually as fuel equivalence ratio increases until to 1.0, at which scramjet combustor is still

430

not unstart and pressure reaches the peak in the first stage heat release region. Heat is added in the

431

diverging area regions and the constant area, leading to the appearance of thermal choking. The

432

position of heat addition has influence on the shock positions and heat addition causes the shock to

433

move forward to a region of smaller area. Fuel allocation proportions injected from the two stages are

434

approximately the same in this case. But in the rear region pressure peak and subsonic flowfield arise

435

when the whole fuel equivalence ratio is smaller. It is explained that constant area regions can contain

436

less combustion heat release than diverging area regions. φ=0.2

φ=0.5

φ=0.7

φ=0.85

48%

φ=1.0

3

area

52%

φ=0.0

2

φ=0.5

1.5

M AN U

2

Ma

2

area/m

pressure ratio

4 3

φ=0.2

φ=0.7

φ=0.85

48%

2.5

φ=1.0

area

52%

area/m2

φ=0.0

SC

5

RI PT

427

1

1

0.5

0.01

0 0

0.2

0.4

0.6

0.8 x/m

1

1.2

1.4

0.004 1.6

0 0

0.2

0.4

0.6

0.01 0.8 x/m

1

1.2

1.4

0.004 1.6

Fig. 18. Numerical Ma steady distributions at different fuel equivalent ratio in case 2

Fig. 17. Numerical normalized wall pressure steady distributions at different fuel equivalent ratio in case 2

Several key points should be identified firstly. Thermal choking is the process in which heat is added

438

to the flowfield until the flow Mach number becomes 1 no matter if initial Mach number greater or less

439

than 1.0. When thermal choking occurs along with back pressures, and the flow can only be adjusted

440

via a shock. The heat release and area change are coupled to decelerate or accelerate the flow. For

441

supersonic flow, diverging area makes pressure decline, and heat release makes pressure rise. When

442

fuel equivalent ratio is less than 0.3, pressure is in the condition of declining in the first stage heat

443

release region, so the effect of area change dominates at the moment. On the contrary the heat release

444

dominates. When fuel equivalent ratio is between 0.3 and 0.5, enormous heat release causes a pressure

445

rise and Mach number decrease. Here, Mach number has not yet met the sonic limit caused by thermal

446

choking. When fuel equivalent ratio is equal to 0.5, it is noted that the flow has choked at the second

447

stage heat addition region. Subsonic flow can be found in the rear part of the combustor, but in the

448

front part flow is supersonic, so combustion is in transition mode. As fuel equivalent ratio continues to

449

increase, pressure peak rises due to stronger thermal choking, and shock wave travels upstream to

450

decelerate the flow. Until fuel equivalent ratio reaches 1.0, the whole flowfield is subsonic, and

451

combustion is in subsonic combustion mode. As a whole, combustion realizes transition from scram

452

mode to ram mode with fuel equivalent ratio increasing. Behind the second stage heat release region,

453

heat release process is over, the flow is accelerated in the constant area and diverging area of the

454

combustor, leading to pressure decreasing and Mach number increasing. The area change shifts the

AC C

EP

TE D

437

17

ACCEPTED MANUSCRIPT 455

flow towards the supersonic in nozzle region. Isolator can tolerate a certain amount of backpressure, so in the flow direction the static pressure

457

decreases and Mach number increases to match combustor entrance condition. Isolator shock train

458

model using Billig’s empirical method is accounted in this paper. In case 2, in the engine isolator

459

entrance Mach number is fixed at 2.5 and it is equal with no accelerating of engine. Therefore, as more

460

heat is continuously added to a combustor, the pre-combustion shock train becomes stronger, and

461

pressure ratio can reach to 3.2.

462

4.3. Control law analysis

RI PT

456

For dual-mode scramjet combustor, there exist at least two kinds of boundaries for control research.

464

Unstart boundary and over-temperature boundary are important safety control boundaries of them.

465

During the ascent process of vehicle, combustor always works close to safety boundary so as to obtain

466

a high performance. The safety boundary codetermined by unstart and over-temperature will be

467

conducted in this section.

SC

463

As is known that the divergent structure is beneficial to accommodate more heat release and improve

469

the combustor performance. In the case of same total fuel equivalence ratio, the larger the first stage

470

fuel equivalence ratio is, the higher the performance of combustor is, the easier unstart phenomenon

471

happens [42]. Hence, this paper only investigates the first stage fuel maximum equivalence ratio when

472 473

M AN U

468

Y, is fixed at 0.5. Incoming flow conditions of combustor are from case 1, 2 and 3 listed in Table 1. Others are obtained by interpolation. In addition, accounting for thermal protection system in passive

cooling combustor, we define the maximum total temperature allowed of gas in combustor outlet as

475

2100K.

TE D

474

At a fixed incoming flow Mach, with the increasing of the first fuel equivalence ratio, working state

477

of combustor gradually tends to unstart. The first stage fuel maximum equivalence ratio allowed which

478

can be found in y-axis is unstart boundary. Like this, with the increasing of the first fuel equivalence

479

ratio, the total temperature of gas in combustor outlet will gradually rise. When gas in combustor outlet

480

is 2100K, the first fuel equivalence ratio is maximum under the over-temperature restrict. It is called

481

over-temperature boundary. For a same incoming flow Mach, the less value of the first stage fuel

482

maximum equivalence ratio is chosen under the limit of unstart boundary and over-temperature

483

boundary. Using this way, the safety boundary is obtained under unstart boundary and

484

over-temperature boundary as shown in Fig. 19.

AC C

EP

476

18

ACCEPTED MANUSCRIPT 0.6 0.5

Over-temperature boundary

Unstart boundary

φ1

0.4 0.3 0.2

2.1

2.3

Ma

2.5

2.7

RI PT

0.1

Fig. 19. Safety working region of combustor under unstart boundary and over-temperature boundary

When Mach is less than 2.3, combustor is mainly limited by unstart boundary. With the increasing of

486

Mach, the first stage fuel maximum equivalence ratio which can be injected in combustor increases

487

gradually. When Mach is between 2.3 and 2.5, the main safety risk of combustor working state transits

488

from unstart to over-temperature. When Mach is greater than 2.5, the first stage fuel maximum

489

equivalence ratio decreases gradually due to the limit of over-temperature boundary.

M AN U

SC

485

490

In conclusion, different fuel equivalence ratio has effect on the distance between current working

491

state and safety boundary. The larger the fuel equivalence ratio is, the nearer the distance away from

492

safety boundary is. But safety boundary of combustor does not change. Only in this way which

493

combustor adopts variable geometry or active cooling, will safety boundary rise up, will safety working

494

region of combustor enlarge.

5. Conclusions

496

In order to develop a control-oriented model to meet the real-time requirements of control design and

497

evaluation, we have followed an approach that restricts the numbers of the physical problem at hand,

498

rather than simply compute a large table of flowfield structure data over a wide range of operating

499

conditions. This paper limits the dimension of the physical problem as much as possible, presents a

500

real-time simulation method. The model is solved at each time step. It means that upon use for control

501

design and evaluation the boundary condition is not limited by pre-specified bounds, which exist in

502

performance table. Flowfield parameter distributions can be computed at a desired operating condition

503

directly. The conclusions are summarized as follows:

AC C

EP

TE D

495

504

(1) A 1-D unsteady flow model is developed and combustion is simplified into mass and heat

505

addition process. This paper presents a heat release distribution obtaining method that can be applied in

506

scramjet combustor 1-D flowfield simulation. Using this method and utilizing several existing

507

experimental data, the heat release distribution which is fit for certain scramjet combustor

508

configuration and is acceptable over a wide range of inflow conditions can be found. This model can

509

simulate transition process between ram-mode and scram-mode operation，and numerical results can

510

agree well with experiments.

511

(2) By mixed programming of C and Matlab, the control-oriented model is built on SIMULINK 19

ACCEPTED MANUSCRIPT platform. When the model runs in a single-core 3.6GHz Intel Core i7-4790 processor, computation

513

time-consuming is found to require 5-10ms within a control period. The computing speed indicates

514

very promising, because model is compact enough to run in real time and it can be used as an

515

embedded model in control system research. Utilizing this model, unsteady and steady flowfield

516

characteristic can be analyzed in the scramjet combustor. The safety boundary codetermined by unstart

517

and over-temperature is conducted in the end of this paper.

Acknowledgments

518

TE D

M AN U

SC

51722601, 91741123).

EP

520

This research work is supported by National Natural Science Foundation of China (Grants No.

AC C

519

RI PT

512

20

ACCEPTED MANUSCRIPT References [1] W.H. Heiser, D.T. Pratt. Hypersonic Airbreathing Propulsion, AIAA Education Series, AIAA, Washington, DC 1994. [2] J.J. Bertin, R.M. Cummings. Fifty years of hypersonics: where we've been, where we're going, Prog. Aerosp. Sci. 39(6) (2003) 511-536.

RI PT

[3] I.D. de Souza, S.N. Silva, R.M. Teles, M.A.C. Fernandes, Platform for real-time simulation of dynamic systems and hardware-in-the-loop for control algorithms, Sensors 14(10) (2014) 19176-19199.

[4] N.N. Smirnov, V.B. Betelin, V.F. Nikitin, Yu.G. Phylippov, Jaye Koo, Detonation engine fed by

SC

acetylene-oxygen mixture. Acta Astronaut. 104 (2014) 134-146.

[5] W. Huang, Z.G. Wang, M. Pourkashanian, L. Ma, D.B. Ingham, S.B. Luo, J. Lei, J. Liu, Numerical investigation on the shock wave transition in a three-dimensional scramjet isolator,

M AN U

Acta Astronaut. 68(11-12) (2011) 1669-1675.

[6] W. Huang, L.Q. Li, L. Yan, L. Liao, Numerical exploration of mixing and combustion in a dual-mode combustor with backward-facing steps, Acta Astronaut. 127 (2016) 572-578. [7] W. Huang, L. Yan, J.G. Tan, Survey on the mode transition technique in combined cycle propulsion systems, Aerosp. Sci. Technol. 39 (2014) 685-691.

[8] W. Huang, L. Yan, Numerical investigation on the ram–scram transition mechanism in a

TE D

strut-based dual-mode scramjet combustor, Int. J. Hydrog. Energy, 41(8) (2016) 4799-4807. [9] W. Huang, L. Ma, M. Pourkashanian, D.B. Ingham, Flow-field analysis of a typical hydrogen-fueled dual-mode scramjet combustor, J. Aerosp. Eng. 25(3) (2012) 336-346. [10] R.W. Weeks, J.J. Moskwa. Automotive engine modeling for real-time control using

EP

MATLAB/SIMULINK, in: International Congress and Exposition, Detroit, Michigan, 1995-0417. [11] N.K. Gupta, B.K. Gupta, N. Ananthkrishnan, G.R. Shevare, I.S. Park, H.G. Yoon, Integrated modeling and simulation of an air-breathing combustion system dynamics, in: AIAA Modeling

AC C

and Simulation Technologies Conference and Exhibit, 2007-6374.

[12] M.L. Liu, H.B. Zhang, Simulation of hypersonic scramjet engine real-time model, J. Aerosp. Power 32(6) (2017) 1447-1455. (In Chinese)

[13] J.A. Roux, L.S. Tiruveedula, Constant velocity combustion parametric ideal scramjet cycle analysis, J. Thermophys. Heat Transf. 30(3) (2016) 697-703.

[14] T. Mitani, T. Hiraiwa, S. Sato, S. Tomioka, T. Kanda, K. Tani, Comparison of scramjet engine performance in Mach 6 vitiated and storage-heated air, J. Propul. Power, 13(5) (1997) 635-642. [15] T. Cui, D.R. Yu, W. Bao. Distributed parameter control arithmetic for an axisymmetrical dual-mode scramjet, Aeronaut. J. 112(1135) (2008) 557-565. [16] T. Cui. Simplified procedure for controlling pressure distribution of a scramjet combustor, Chin. J.

21

ACCEPTED MANUSCRIPT Aeronaut. 27(5) (2014) 1137-1141. [17] J. Rajasekaran, A.Maity, A. Lal, A. Lal, R. Padhi. Nonlinear control design for an air-breathing engine with state estimation, in: American Control Conference, IEEE Press, 2009, 913-918. [18] K.P.B. Chandra., N.K. Gupta, N.Ananthkrishnan, I.S. Park, H.G. Yoon, Modeling, simulation, and controller design for an air-breathing combustion system, J. Propul. Power, 26(3) (2010) 562-574.

RI PT

[19] T.R.A. Bussing, E.M. Murman. A one-dimensional unsteady model of dual mode scramjet operation, in: AlAA 21st Aerospace Sciences Meeting, Reno, Nevada, 1983-0422.

[20] J.H. Liu, W.H. Ling, X.Z. Liu, L. Liu, Z. Zhang, A quasi-one dimensional unsteady numerical analysis of supersonic combustor performance, J. Propul. Techo. 19(1) (1998) 1-6. (In Chinese)

SC

[21] L. Wang, J.W. Xing, Z.H. Zheng, J.L. Le, One-dimensional evaluation of the scramjet flowpath performance, J. Propul. Techo. 29(6) (2008) 641-645. (In Chinese)

[22] D.S. Niu, L.Y. Hou., P.F. Pan. Quasi-one-dimensional model of supersonic a combustor with

M AN U

various fuels, J. Tsinghua Univ.: Nat. Sci. Ed. 53(4) (2013) 567-572. (In Chinese) [23] A.H. Shapiro. Dynamics and Thermodynamics of Compressible Fluid Flow, Ronald Press, New York, 1953.

[24] R.P. Starkey, M.J. Lewis. Sensitivity of hydrocarbon combustion modeling for hypersonic missile design, J. Propul. Power, 19(1) (1971) 89-97.

[25] T.F. O’Brien, R.P. Starkey, M.J.Lewis. Quasi-one-dimensional high-speed engine model with

TE D

finite-rate chemistry, J. Propuls. Power, 17(6) (2001) 1366-1374. [26] S.M. Torrez, J.F. Driscoll, M. Ihme, M.L. Fotia, Reduced-order modeling of turbulent reacting flows with application to ramjets and scramjets, J. Propuls. Power, 27(2) (2011) 371-382. [27] S.M. Torrez, J.F. Driscoll, M. A. Bolender, M.W. Oppenheimer, D.B. Doman, Effects of

EP

improved propulsion modelling on the flight dynamics of hypersonic vehicles, in: AIAA Atmospheric Flight Mechanics Conference and Exhibit, Honolulu, Hawaii, 2008-6386.

AC C

[28] S.M. Torrez, J.F. Driscoll, D. J. Dalle, M.A. Bolender, D.B. Doman, Hypersonic vehicle thrust sensitivity to angle of attack and Mach number, in: AIAA Atmospheric Flight Mechanics Conference, Chicago, Illinois, 2009-6152.

[29] D.J. Dalle, S.G.V. Frendreis, J.F. Driscoll, C.E.S. Cesnik, Hypersonic vehicle flight dynamics with coupled aerodynamics and reduced-order propulsive models, in: AIAA Atmospheric Flight Mechanics Conference, Toronto, Ontario Canada, 2010-7930. [30] C. Birzer, C.J. Doolan. Quasi-one-dimensional modeling of hydrogen fueled scramjet combustors, J. Propul. Power, 25(6) (2007) 1220-1225. [31] T. Lu, L.H Chen, Q. Chen, F. Li, X.Y Chang, Quasi-one-dimensional multimodes analysis for dual-mode scramjet, J. Propul. Power, 30(6) (2014) 1559-1567. [32] S.M. Torrez, D.J. Dalle, J.F. Driscoll. New method for computing performance of choked reacting 22

ACCEPTED MANUSCRIPT flows and ram-to-scram transition, J. Propul. Power, 29(2) (2015) 433-445. [33] R.F. Cao, T. Cui, D.R. Yu, J.T. Chang, W. Bao, Z.Q Wang, New method for solving one-dimensional transonic reacting flows of a scramjet combustor, J. Propul. Power, 32(6) (2016) 1-10. [34] M.J. Zucrow, J.D. Hoffman. Gas Dynamics, Wiley, New York 1976.

Springer Science and Business Media, New York 2013.

RI PT

[35] E.F.Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction,

[36] P.J. Waltrup, F.S. Billig. Prediction of precombustion wall pressure distributions in scramjet engines. J. Spacecr. Rockets 10(9) (1973) 620-622.

[37] Q.C. Yang, J.T. Chang, W. Bao, Thermodynamic analysis on specific thrust of the hydrocarbon

SC

fueled scramjet, Energy, 76 (2014) 552-558.

[38] A.C. Zambon, H.K. Chelliah. Explicit reduced reaction models for ignition, flame propagation, and extinction of C2H4/CH4/H2, and air systems, Combust. Flame 150(1-2) (2007) 71-91.

M AN U

[39] The Math Works Inc. Simulink Coder Target Language Compiler, The Math Works Inc. Natick, MA 2012.

[40] C.L. Zhang, J.T. Chang, Y.S. Zhang, Y.Y. Wang, W. Bao, Flow field characteristics analysis and combustion modes classification for a strut/cavity dual-mode combustor, Acta Astronaut. 137 (2017) 44-51.

[41] Z.G. Jin, K.Y. Zhang, Y. Liu, Design and performance investigation of sidewall compression

Chinese)

TE D

scramjet inlets operating from Mach 4 to Mach 7, J. Aerosp. Power 26(6) (2011) 1201-1208. (In

[42] J.C. Hu, J.T. Chang, W. Bao, Q.C. Yang, J. Wen, Experimental study of a flush wall scramjet combustor equipped with strut/wall fuel injection, Acta Astronaut. 104(1) (2014) 84-90.

EP

[43] J.L. Zhang, J.T. Chang, J.C. Ma, C.L. Zhang, W. Bao, Investigation of flame establishment and stabilization mechanism in a kerosene fueled supersonic combustor equipped with a thin strut,

AC C

Aerosp. Sci. Technol. 70 (2017) 152-160. [44] J.T. Parker, A. Serrani, S. Yurkovich, M.A. Bolender, D.B. Doman, Control-oriented modeling of an air-breathing hypersonic vehicle, J. Guid. Control Dyn. 30(3) (2015) 856-869.

[45] J. Jiang, M. Chu, X. Xu. A Quasi-one-dimensional method for prediction of dual mode scram-jet combustor performance, J. Propul. Techo. 34(6) (2013) 802-808. (In Chinese)

[46] Y.X. Zhang, Z.G. Wang, M.B. Sun, Y.X. Yang, H.B. Wang, An investigation on one-dimensional model of heat release distribution in scramjet combustor, J. Propul. Techo. 36(12) (2015) 1852-1858. (In Chinese) [47] K.W. Liu, Y.J. Zhu, J.M. Yang, X.B. Mao, Y.C. Wu, Comparison of two typical flow-parameter-matching schemes in a combustion wind tunnel, J. Propul. Techo. 38(6) (2017) 1226-1234. (In Chinese) 23

ACCEPTED MANUSCRIPT [48] D. Zhang, Y. Feng, S.L. Zhang, J. Qin, K.L Cheng, W. Bao, D.R. Yu, Quasi-one-dimensional model of scramjet combustor coupled with regenerative cooling, J. Propul. Power 32(3) (2016) 1-11. [49] N.N. Smirnov, V.B. Betelin, V.F. Nikitin, L.I. Stamov, D.I. Altoukhov, Accumulation of errors in numerical simulations of chemically reacting gas dynamics. Acta Astronaut. 117 (2015), 338-355. [50] N.N. Smirnov, V.B. Betelin, R.M. Shagaliev, V.F. Nikitin, I.M. Belyakov, Yu.N. Deryuguin, S.V.

RI PT

Aksenov, D.A. Korchazhkin. Hydrogen fuel rocket engines simulation using LOGOS code. Int. J.

AC C

EP

TE D

M AN U

SC

Hydrog. Energy, 39 (2014), 10748-10756.

24

ACCEPTED MANUSCRIPT

Highlights 1.

Control-oriented 1-D unsteady model was presented for dual-mode scramjet combustor.

2.

By mixed programming, real-time model was built on SIMULINK platform.

AC C

EP

TE D

M AN U

SC

RI PT

3. The heat release distribution obtaining method was discussed and explained. 4. The unsteady, steady flowfield process and safety boundary were analyzed.