ESPSS SIMULATION PLATFORM J. Moral(1), R. Pérez Vara(1), J. Steelant(2), M. de Rosa(2) (1) Empresarios Agrupados Internacional. S.A. Magallanes, 3. 28015 Madrid. Spain. (email: [email protected]) (2) Aerothermodynamics and Propulsion Analysis Section (TEC-MPA), ESTEC-ESA, Keplerlaan 1, P.O. Box 299, 2200 AG Noordwijk, The Netherlands (email: [email protected])
Abstract This paper describes the EcosimPro Libraries provided with the ESPSS (European Space Propulsion System Simulation) software, which is concerned with setting up a common European platform for propulsion system simulation. A description of the components and the most important aspects of the physical modelling related with the libraries’ components are included. Three groups of libraries are presented: 1) libraries providing a standard properties database for propellants, pressurants and other fluids and materials, 2) general purpose libraries for two-phase fluid systems including tanks and turbo machinery, and 3) libraries for the simulation of combustion chambers and cooling systems. Hydraulic or pneumatic systems with coupled heat transfer networks and controllers can be easily represented with these libraries thanks to the EcosimPro object oriented graphical capabilities. Cavitation and priming phenomena under two-phase flow (with or without a non-condensable gas travelling in a liquid) can be calculated in pipes or in other components. Models where one or more combustion chambers are present (staged engines) can be evaluated allowing to analyze in great detail the transients during the start up and shut down processes, where the valve sequences are decisive. Three application examples are presented below, one simulating a water hammer case, another related with a pressurisation system and the third one with the simulation of a staged combustion cycle, showing some of the capabilities of the ESPSS libraries.
INTRODUCTION ESPSS is an ESA Software requiring the coordination of multiple disciplines. The combination of an expert in the design and production of software simulation tools (Empresarios Agrupados Int.) together with the main industrial architect and stage integrator of Ariane launchers (EADS Astrium – Space Transportation), a company at the forefront of CFD and Optimization methods (CENAERO) and a company with a large experience in Satellite Propulsion Systems of Chemical and Electrical type (KopooS) provided the best guarantee for the successful development of the project Three main targets were in the base of ESPSS V1 being successfully fulfilled: • To develop a set of EcosimPro Libraries with components and properties functions for the simulation of spacecraft and launch vehicle propulsion systems.
• To share a basic tool (kernel) providing modular model building, a standard interface and a standard database of fluid and material properties. • Finally, to be able to interact with other softwares, especially concerning the following areas: Optimization for Design and Test Fitting and Coupling to External Software. Current version 2 accounts for upgrades and experience acquired under the Industrial Evaluation phase: • Experimental validation of priming cases, two-phase tanks filling and the ESCA engine • New capabilities developed for ESPSS V2: Liquid injection in combustion chambers, upgrade of pressurization models, fluid properties, etc ESPSS EcosimPro libraries provide palettes of components that can be used to build graphically complex systems. EcosimPro basically incorporates the following features: a object-oriented programming language which enables encapsulation, inheritance and aggregation, a powerful DAE solver and a friendly Graphic User Interface (GUI).
ESPSS LIBRARIES A) Fluid_properties Library. The most important features are summarized as follows: • Most of the fluids used for rockets applications are available. Fluids are supported in different categories depending on the type used: 1) Perfect gases (transport and heat capacity properties obtained from CEA, [2], coefficients (temperature dependent)): Cp = a1 ⋅ T − 2 + a2 ⋅ T −1 + a3 + a4 ⋅ T + a5 ⋅ T 2 + a6 ⋅ T 3 + a7 ⋅ T 4 R T T C p (T ) dT − R log(P / P0 ) h = h 0 + ∫ C p (T )dT ; s = s0 + ∫ T T0 T0
2) Perfect gases (transport and heat capacity properties interpolated in files (temperature dependent)): C p / R = Table _ Interp(T )
3) Simplified liquids interpolated in tables (properties not pressure dependent). The state equation is simplified as follows
κ
ρ ( P, T ) =ρ (T )[1 + κ ( P − Pref )]
Where is numerically calculated from the density, sound speed and specific heat interpolated in user data tables depending on T. 4) Real fluids properties are directly interpolated in external 2D properties tables (normally obtained running the REFPROP code, [1]) considering liquid,
superheated, supercritical or two-phase (temperature and pressure dependent)
flow
FLUID_PROPERTIES library. The actual fluid phase will be automatically calculated.
• The homogeneous equilibrium model is used to calculate mixtures of a real fluid with a noncondensable gas (quality, void fraction, etc) in case of two phase flow assuming a relative humidity of one. The following state equations are involved:
• Flow inversion, inertia, gravity forces and high speed phenomena are considered in pipes, volumes and junctions.
ρ liq , u liq = f state ( fluid , P, T ); P = Pnc + Pvap ρ vap , u vap = f sat ( fluid , Pvap ); Pvap = f sat ( fluid , T ) “f” denote the corresponding pure fluid functions. In this system of equations, Pnc and T are unknowns. Assuming that the global density, ρ, the noncondensable mass fraction, xnc, and the mixture energy, u, are known (dynamic variables in the capacitive fluid components), the following closing equations allow an iterative calculation of the homogeneous temperature, the non condensable pressure and the void fraction α:
• Heat transfers at the walls/fluid interface and concentrated (valves) or distributed (pipes) load losses calculation is included in the formulation. • 1D Pipe incorporates an area-varying non-uniform mesh 1-D spatial discretization into n (input data) volumes using some of the most up-to-date, robust and accurate Roe or centred schemes. Cavitation and priming phenomena under two-phase flow can be calculated in pipes or other components. The components of this library are listed below, where the inheritance hierarchy is shown:
ρ ncα = ρ x nc ; ρ nc = f state ( fluid nc , Pnc , T ) u = (1 − x nc )( x u vap + (1 − x)u liq ) + x nc u nc where
α = ( ρ liq − (1 − x nc ) ρ ) /( ρ liq − ρ vap ) x = α ρ vap /( ρ liq − α ( ρ liq − ρ vap ))
• User-defined fluids are available. An EXCEL tool generating Peng Robinson real fluid files (MMH and N2H2 files, see [3] to [5]) is included. Available real fluids are summarized below: Fluid
Crit. Cond. Pbar) T (K)
P/T Range Pr. (bar) Temp. (K)
ParaH2
12.84
32.94
0. – 1200
13.8 – 400
NormalH2
13.03
33.24
0. – 500
18.0 – 700
Dodecane
18.17
658.1
0.– 3500
263.6 – 1000
Methane
45.99
190.6
0 – 10000
90.7 – 625
Propane
42.47
369.8
0 – 1030
85.5 – 623
Butane
37.96
0. – 690
134.9 – 637
MMH
80.64
425.1 567.
0 – 350
220. – 620
Hydrazine
149.4
653.5
0 – 400
182.3 – 710
Oxygen
50.43
154.6
0. – 820
54.4 – 1000
N2O
72.45
309.5
0 – 500
182.3 – 525
N2O4
101.3 61.48
431.15
0 – 300
261.9 – 600
Ethanol
513.9
0. – 2800
250 – 770
Helium
2.27
5.19
2.18 – 1500
Nitrogen
33.96
126.2
0. – 1000 0. – 22000
63.15 – 2000
Xenon
58.42
289.7
0. – 7000
161.4 – 750
Argon
48.63
150.7
0. – 10000
83.8 – 700
Water
220.6
647.1
0. – 10000
273.2 – 1275
se113.3
405.4
0. – 10000
195.5 – 700
R12
41.36
385.1
0. – 2000
116.1 – 578
R22
49.90
369.3
0. – 600
115.7 – 554
R134A
40.59
374.2
0. – 700
169.8 – 561
R141B
42.12
477.5
0. – 4000
169.7 – 716
Ammonia
B) Fluid_flow_1d Library. This library allows 1D transient simulation of two-fluid, two-phase systems. The most important features are the following: • The working fluid(s) can be easily selected from a large collection of fluids included in the
In an ESPSS fluid network, every component is either a resistive component or a capacitive component. A resistive component receives the state variables (pressure, density, velocity, chemical composition and enthalpy) as input and gives back the flow variables (volumetric, mass and enthalpy flows) as output. A capacitive component receives the flow variables as input and gives back the state variables at output. Pipe component behaves as a chain of C-R-C-R …-C elements. Special components as Tees behave as volumes with each one of its ports ending with R elements. This computational scheme prevents the appearance of algebraic loops and high index DAE (Differential Algebraic Equations) in the mathematical model of the pipe network. The conservation equations include under the same formulation gas, liquid and two-phase flow regimes for ideal or real fluids. Below, for illustrative purposes, a summary of the Pipe equations: ⎞ ⎛ ρv ⎛ρ ⎞ ⎟ ⎜ ⎜ nc ⎟ nc ⎟ ⎜ ρ vx ∂ω ∂ f (ω ) ⎜ρx ⎟ + = Ω(ω ); ω = A ⎜ ⎟ ⎟; f (ω ) = A ⎜ 2 ∂t ∂x ρ v ⎜ ρv + P⎟ ⎟ ⎜ ⎜ ρu ⎟ ⎟ ⎜ ρ vh ⎠ ⎝ ⎠ ⎝ ⎛0 ⎞ ⎜ ⎟ ⎜0 ⎟ Ω(ω ) = ⎜ − 0.5(Δξ / Δx) ρ v v A − ρ g A + P(ΔA / Δx) ⎟ ⎜ ⎟ ⎜ ΔQ / Δx − ρ g vA ⎟ ⎝ ⎠
ρ, P, h (=u+P/ρ) are the gas/liquid mixture density, pressure and total enthalpy respectively. A is the flow area. x, xnc are the mass fractions and v the velocity. The different source terms are the following: Source terms • The equivalent distributed friction, Δξ(i) is calculated using empirical correlations of the friction losses, elbows pressures drops, etc [6, 7] • g represents the body forces. It is computed as the scalar product of an extended “gravity” vector (gx, gy, gz) with the direction of the pipe:
extra pressure drop coefficient, ζcrit, in such a way it forces the sonic flow under critical conditions: ζ crit = (G st / G crit ) 2 ; G st = 2 ρ up (( P + 0.5ρv 2 ) 1 − ( P − 0.5ρv 2 ) 2 )
“Gcrit” is the sonic flow per unit of area (ρc)crit calculated exactly by iteration even under two-phase fluxing flow [12] or approximately [13] C) Tanks Library. It is an important complement of the FLUID_FLOW_1D library for the simulation of rocket engine and spacecraft two-phase two fluid tanks. Several models of Tanks are available:
Z2 ΔP =ρ·g ΔZ Z1
• The term ΔQ(i) will permit the exchange of heat with the walls:
⎧ Pinner (tp _ in.T (i ) − Ti ) + ⎫ , ΔQi = hc film ,i Δx i ⎨ ⎬ ⎩ Pouter (tp _ out.T (i ) − Ti ) ⎭ where Pinner and Pouter are the wet perimeters; tp_in, tp_out are thermal ports connected to the tube with the same number of nodes as the fluid vein. The port temperatures tp_..T(i) behave as the wall internal temperatures, to be determined in the connected THERMAL diffusive components. The film coefficient hc is calculated using empirical two-phase correlations [8, 9].
• Wall compressibility. Assuming linear elasticity and fully anchored pipe wall: (dA / dt ) wall = A κ wall (∂P / ∂t ); κ wall =
(
Din ⋅ 1− μ 2 th ⋅ M E
)
where ME is the Young's modulus of elasticity, μ is the Poisson’s ratio and th is the wall thickness. Numerical schemes. Pressures and temperatures are associated with n nodes. The dynamic mass flows are calculated either at the internal junctions (centred staggered scheme, where each junction has associated two half volume inertias), or at the volume center as for the other magnitudes (ROE upwind scheme, [10, 11]): Node 1
Node 2
Nodal
Port 1
Junction 1
volume (i)
Node n
Port 2
Junction 2
Junction (j)
Junction n+1
Pressure, temperature and quality calculation At each discretized volume, the non derivative state variables (pressures, qualities and temperatures) will be calculated using thermodynamic functions and the homogeneous equilibrium model if necessary (see the FLUID_PROPERTIES library). The sonic flow limitation Valves and Junctions components use a momentum equation similar to that of the Pipe, but in this case the mass flow being limited to the sonic flow adding an
• Tank_Single component allows calculating the liquid level evolution assuming that the gas mixture (vapor plus a non-condensable) is separated from the liquid at one equilibrium liquid/gas temperature. • The Tank_Bladder component assumes two different fluid temperatures and pressures, one for the liquid one for the gas mixture. It also considers a bladder separating the liquid from the gas. It can be used as an accumulator or a pressurized tank. • 1D Tank models calculate the temperature variation along a “vertical” axis. Bubbles rising to the gas part (by boiling inside the liquid, Tsat>Tliq) is included in the formulation. Local boiling at the wall liquid interface is only considered for the calculation of the heat transfer coefficient but the vapor produced is supposed to be condensed in the liquid making higher its temperature. The heat and mass transfer at the liquid/gas interface is also calculated. Tank models can simulate the liquid level passing through a fluid port in such a way that the exiting flow can be liquid or gas depending on the liquid level. 1D Tanks conservation equations Inner level of coding corresponds to the 1D (liquid or gas) cavities. Mass, energy and momentum equations are programmed similarly than in Pipes (centred scheme) but using a moving grid. An abstract component is built topologically with two 1D fluid cavities, one simulating the liquid, and one for the gas. Mass & energy exchange equations at the interface are programmed here: assuming a very thin layer at Tsat, the conservation equations allow calculating the evaporation (positive or negative) flow: m vap = A(hc g (Tg − Tsat ) + hc l (Tl − Tsat ))/(h sat_vap − h sat_liq ) q gas_to_liq = A·hc g (Tg − Tsat ) − m vap ·h sat_vap
Tg, Tl are the gas liquid side temperatures. The saturation temperature is calculated at the partial vapor pressure of the gas volume. hcg, hcl are the gas and the liquid side heat exchange coefficients.
Wetted areas calculation: Finally, the calculation of the wetted areas and other geometrical parameters needed to the grid definition is made throughout calls to the Tank_Shapes functions calculating the surface shapes. Cover concept:
matrix
calculation of the heat exchange between a moving fluid grid and a fixed wall grid
Gas Gas Gas Gas Gas Gas
Node Node Node Node Node Node
5 4 4 3 2 1
Upper Dome Node 4 Upper Dome Node 3 Upper Dome Node 2 Upper Dome Node 1
Liquid Node 5 Cylinder Node 2 Liquid Node 4
Cylinder Node 1
Liquid Node 3 Lower Dome Node 4 Liquid Node 2 Liquid Node 1
Lower Dome Node 3 Lower Dome Node 2 Lower Dome Node 1
Several diffusive thermal components with 2D (radial and axial) spatial discretization simulate the walls accordingly with the tank geometry. The “Sphere_ins”, “Cylinder_ins” and “Dome_ins” components can simulate a layer of insulation covering the 1D Tanks. Fluid Separator component: It allows the simulation of a common circuit pressurizing different tanks with different liquids. The component includes a non-return valve avoiding the contact of different kinds of vapors connected by a common pressurization system, see the example below. It is pointed out that the common pressurization circuit can use real properties of gas which can be decisive calculating the behavior of the tank at high pressures and low temperatures.
Programming of the turbo-machinery components is non-dependent on the working fluid type. Components are built topologically with a Volume component receiving/giving the head rise and the mechanical work. A continuous block contains the equations used to calculate the head and the mechanical work (power) based on the performances maps: Turbine/Compressor maps a) User-defined maps are introduced by means of 2D input data tables. Independent variables are: Speed coefficient: N = r ·ω / v son Total pressure ratio: Π = Pin / Pout Dependent variables: Mass flow coefficient: Q + = m ⋅ v son /(r 2 Pin ) Specific torque: ST = T /(r m ⋅ vson ) The turbine model is very general because the parameters are scaling with the mean radius and the speed of sound. Therefore, given turbine characteristics remain valid for geometrically similar turbines, and for different fluids (in first order approximation). b) Generalized maps (no maps available) estimate 1 the efficiency as a function of the reduced speed and pressure ratio parameters: ⎛ 0.5·π ·( N / N ) ⎞ . ⎛ Π −1 ⎞ nom ⎟ ⎟⎟ ·sin⎜ eta = etanom ·⎜⎜ tt ⎜ (Π − 1) /(Π − 1) 0.5 ⎟ ⎝ Π nom − 1 ⎠ tt nom ⎝ ⎠ 0.5
(
)
GasTurbo_efficiency 0.9 0.8
• Pump model provided with user-defined turbo-pump characteristic curves adapted to positive and negative speeds and flow zones. • Turbine and Compressor components provided with user-defined performance maps as a function of the reduced axial speed and pressure ratio. • Special turbo machinery components allowing simple calculation of generalized performance maps as a function of the nominal performances and other significant design data. Two different types of pumps/turbines/compressors are then available: one “generic” model (with a G on the symbol) if the off-design characteristics are unknown, one specific model which can only be used with tables, for well-defined turbo machinery.
PI=12
0.6
Efficiency
D) Turbo-Machinery Library. It is an important complement of the FLUID_FLOW_1D library for the simulation of pumps, turbines and compressors. The most important features are the following:
PI_NOM = 10 N_NOM = 0.8 eta_NOM=0.8
PI=10
0.7
PI=6
0.5 0.4 0.3
PI=3
0.2 0.1 0 0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
N
Depending on the pressure and flow values, it is possible to distinguish four different cases:
• Pout > Pin (Compressor) 1. Positive mass flow: the power absorbed by the fluid is calculated as W = m·( hise ,out − hin ) / eta 2. Negative mass flow: the power taken from the fluid is W = −eta·m·( hout − hise ,in )
• Pout < Pin (Turbine)
1 Following simple functions can be easily changed by the user of a proof knowledge of turbo machinery is available
3. Positive mass flow: the power taken from the fluid is calculated as W = eta·m·( H ise ,out − hin ) 4. Negative mass flow: the power absorbed by the fluid is W = − m·( hout − H ise ,in ) / eta The pressure ratio, dp_rel, is a function of the beta parameter (blade angle), the reduced speed and the inlet Mach number: 2 ⎡ M − M0 ⎤ tan( β )⎥ dp _ relturb = sign( M in − M 0 ) ⎢ in ⎣ N nom ⎦ ⎡ M − M0 ⎤ tan(β )⎥ dp _ relcompr = − k ⋅ sign( M in − M 0 ) ⎢ in ⎣ N nom ⎦
2
In both cases the speed contribution its self is calculated as follows:
⎛ N dp _ rel = dp _ rel − 0.2 * sign( N ) ⎜⎜ ⎝ N nom
⎞ ⎟⎟ ⎠
2
coefficients [2] for an arbitrary mixture or chemical’s reactants. The equilibrium molar fractions at of a mixture of reactants are derived from the Minimum Gibbs energy method [15]
• Non adiabatic 1D combustor components: The equilibrium combustion gases are calculated using previous capabilities and single models for the liquid droplets evaporation • Cooling jacket components including the injection torus and a 3D model of wall temperatures. • Modeling of nozzles provided with a 1D spatial discretisation, solid propellant starters, igniters and thermal coating protection are available using combustor components. The available components are:
Pump Performances a) User-defined maps are introduced by means of input data tables. Independent variables are: Mass flow coefficient:
ϕ + = m /( ρ out ω )
(m3)
Dependent variables:
ψ + = ΔP /( ρ out ω 2 ) C + = T /( ρ out ω 2 ) (m5)
Head rise coefficient: Reduced torque:
(m2)
The pump model described here has the disadvantage of being not dimensionless. Hence, the characteristics differ for each pump, even for geometrically similar pumps. b) Generalized maps. In this case, the pump performances (head and torque) are built in by means of fixed 2D data tables covering all zones of operation and defined as functions of the specific speed Ns and a dimensionless variable θ that preserves homologous relationships in all zones of operation (see [14], Chapter 7.2): θ = π + tan −1 (v / n) ; Ns = rpm Q / n s (TDH / n st ) 0.75 ns is number of suctions, nst the number of stages, Q is the volumetric flow and TDH the total dynamic head of the pump. “v” and “n” are the reduced flow and speed variables: v = Q / Q R n = 30ω / π / rpm R
The COMB_CHAMBERS components can be connected with other ESPSS components for the simulation of a complete rocket engine cycle. Models where one or more chambers are present (staged engines) can be evaluated. This library permits to analyze in great detail the transients during the start up and shut down processes, where the valve sequences are decisive. Pre-burners and main combustion chambers are topologically built by means of a combustor, two injectors and two cavity components. They can work either under liquid, two-phase or gas injection conditions:
Dependent variables are the dimensionless head and torque parameters respectively:
wh =
TDH T ; wt = 2 2 2 TDH R (n + v ) TR ( n + v 2 )
Mechanical Balance The mechanical balance allows the calculation of the axial speed dynamically:
I mech ⋅ ω ' = T shaft , port − T ω is the mechanical speed Imech is the mechanical inertia and T is the torque, calculated as W/ω. E) Comb_Chambers Library. The most important features are the following:
• The properties of the combustion gases (transport and heat capacity) are obtained from the CEA
1D Combustor conservation equations The numeric method used for the resolution of the subsonic sections of a combustor is based on the transient conservation equations. Using this transient approach, the number of implicit equations is reduced (the state variables being dynamic). On the other way
wall heat exchange, pressure drops and vaporization terms are more easily taken into account: a) Gas mixture mass and energy equations:
Vi ρ´i = m jun ,i −1 − m jun ,i + mvap _ red ,i + mvap _ oxy ,i (u´i ρ i + u i ρ´i ) Vi = mh jun ,i −1 − mh jun ,i + q wall ,i + q vap ,i “i” is the volume number, “jun” is the junction between the volume i and i+1. • At injector level, the gas mass and enthalpy flows are calculated by the Injector components taking into account the current quality in Cavities components. Injector and cavity components (see figure above) are Fluid_flow_1D components. • Total enthalpies “hjun” are calculated using the upstream cell conditions: hjun,i = hi-1.= (u + P/ρ)i-1; u = ust + v2/2 • Node velocities “v” are calculated using the adjacent junction mass flows and the mean junction densities: m jun,i ⎞ i=2,n; 1 ⎛ m jun ,i −1 ⎟ vi = ·⎜⎜ + Ai ⎝ ρ i + ρ i −1 ρ i + ρ i +1 ⎟⎠ b) Droplets vaporization model (mvap, qvap terms). Assuming a very thin saturated layer between the droplets and the surrounding gases, the sum of convective heat plus enthalpy mass flow are the same at both sides of the layer: m vap_red,i = A liq_red,i
hc i (Ti − Tsat_red,i ) + hc liq_red (Tliq_red,i − Tsat_red,i )
ignition flag activated. Otherwise the burning factors are set to zero. d) Vapors / non-condensable mass equations: The vapors and non-condensable mass conservation equations take into account the burned gases production and the vaporization terms previously calculated:
( x´red ,i ρ i + x red ,i ρ´i ) Vi = x red ,i −1 m jun ,i −1 − x red ,i m jun ,i + m vap _ red ,i − x red ,i mbu ,i ( x´oxy ,i ρ i + x oxy ,i ρ´i ) Vi = x oxy ,i −1 m jun ,i −1 − x oxy ,i m jun ,i + mvap _ oxy ,i − x oxy ,i mbu ,i ( x´nc ,i ρ i + x nc ,i ρ´i ) Vi = x nc ,i −1 m jun ,i −1 − x nc ,i m jun ,i The mass flow of the injected vapors (mred, moxy, see before) and injected non-condensable gases will be added as source terms to the respective conservation equation of the first chamber volume. In the same way, the mass flow of the solid propellant gases, “starter_m”, will be added as source a term to the solid propellant gases conservation equation of the first chamber volume. e) Liquids conservation equations: The liquid mass and enthalpy flows are calculated assuming that Vel_gas=Vel_liq and neglecting Cp_liq derivatives: M 'liq _ red ,i = M liq _ red , i −1veli −1 / Li −1 − M liq _ red , i veli / Li − mvap _ red ,i ' ( M ⋅ T )liq _ red , i = ( M ⋅ T ) liq _ red , i −1 veli −1 / Li −1 − ( M ⋅ T ) liq _ red , i veli / Li
h vap_red,i − h liq_red,i
q vap_red,i = m vap_red,i h vap_red,i − A liq_red-gas,i hc i (Ti − Tsat_red,i )
hc liq_red = 2λliq _ red ,i / Ddroplet _ red Aliq_red,i = f vap ,i 6 M liq _ red ,i / ρ liq _ red / Ddroplet _ red (Similar set of equations for the oxidizer fluid)
Ddroplet is the mean droplets diameter and Aliq,i is the equivalent exchange area between the droplets and the gas. We point out that the droplet diameter has been modulated by the vaporization factor. This factor is an input data depending on the time and on the volume number. In theory, assuming a known droplet size at the injection plate, the droplet diameter evolution could be determined by “simple” equations relating the evaporated mass flow with the liquid mass conservation equations. Nevertheless, due to the high penetration and breakup of the liquid jets, it seems more realistic to assume a known (adjusted) droplet size at each chamber volume, the number of droplets being determined by the current liquid mass. Combustor models include an automatic limitation of the vaporization rate when the chamber temperature is below the saturation temperature. c) Burning rate. The burned gases mass flow is calculated assuming a global characteristic burning time. It is also supposed that any specie (vapor or burned gas) present in gas mixture contributes to the global reaction rate, so the burning rate will be proportional to the total gas mixture density:
−qvap _ red ,i / Cpliq _ red ,i M 'liq _ oxy ,i = M liq _ oxy ,i −1veli −1 / Li −1 − M liq _ oxy ,i veli / Li − mvap _ oxy ,i ' ( M ⋅ T )liq _ oxy , i = ( M ⋅ T ) liq _ oxy , i −1 veli −1 / Li −1 − ( M ⋅ T ) liq _ oxy , i veli / Li
−qvap _ oxy ,i / Cpliq _ oxy , i
where Li is the length of volume no. i. f) Momentum equations are basically the same as in the Pipe component:
[
0.5( Li −1 + Li )·m 'jun ,i = Ai −1 P + qn + ρ ·v 2 − 0.25ξ ·ρ ·v v
[
− Ai P + qn + ρ ·v + 0.25ξ ·ρ ·v v 2
]
i
]
i −1
− 0.5( Pi −1 + Pi )( Ai −1 − Ai )
The artificial dissipation, qn(i), is calculated as follows: qn(i ) = − Damp
m jun (i + 1) − m jun (i ) A
V sound (i )
Damp is a global input data (0.3 – 1) of the FLUID_FLOW_1D library. Momentum equations are applied to the exit of any volume in which the combustor is discretized with the exception of the last volume that should end with the throat to avoid numerical problems (transition from subsonic to supersonic flow). The outlet mass flow will be calculated by the inherited components (preburners) or by a nozzle connected to the combustor. Combustion gases properties calculation It is supposed that any molar fraction follows a global reaction rate accordingly with the previously mentioned burning time:
mbu ,i = f bu ,i ρ iVi / τ bu
y 'k _ bu ,i = ( y k _ eq ,i − y k _ bu ,i ) / τ bu
where τbu is the characteristic burning time and fbu,i is the burning factor at volume no. i. The burning factors are automatically set to one if the burning conditions are true: mixture ratio within the allowed limits and
For each chamber volume i, the equilibrium composition (needed for the calculation of the actual burned gases composition) is calculated using the Minimum Gibbs energy method [2, 15] as a function of the gas mixture molar fractions, pressure and
enthalpy. The gas mixture molar fractions are calculated as follows: Reducer vapors contributions: y k _ red N k ,i = x red ,i ; MWmix , red = ∑ y k _ red MWk , red MWmix , red k =1, Nchem Oxidizer vapors contributions: N k ,i = N k ,i + x oxy ,i
y k _ oxy MWmix ,oxy
; MWmix , xoy =
∑y
MWk ,oxy
∑y
MWk _ bu
k _ oxy k =1, Nchem
Burned gases contribution: N k ,i = N k ,i + xbu ,i
yk _ bu ,i MWmix,bu
;
MWmix,bu =
k _ bu ,i k =1, Nchem
Once the number of moles of the mixture has been evaluated, and using the “dynamic” enthalpy value obtained from the conservation equations, it is possible to call to the Minimum Gibbs energy method to obtain the equilibrium combustion gases composition:
(yk_eq , Teq)i = fminGibbs(Nk,i , hi–vi2/2, Pi) Two possibilities are foreseen calling the Gibbs function: A/ Equilibrium and B/ frozen flow. In the last case (no ignition) the molar fractions remain constant: yk_eq,i = Nk,i. The effective combustion gas constants (Ri, Cpi, condi, visci) will be derived using the mixture properties equations as a function of yk_bu,i. The pressure is obtained from the perfect gas equation: P i = ρ i ·R i ·T i . Using the Combustor_eq component (no burning delay times), the pressure equation and the Minimum Gibbs function would become an algebraic loop because the actual molar fractions (yk_bu,i) of the mixture are not dynamic variables. Heat exchanged with walls: The term qwall appearing in the energy conservation equation permits the exchange of heat through the walls (that can be represented by THERMAL components or by the Cooling Jacket component): q wall ,i = h film ,i Awet ,i (Ti ,aw − tp.T (i )) + σ ⋅ Awet ,i (Ti 4 − tp.T (i ) 4 )
tp is the name of the thermal port (with n nodes in axial direction) connected to the Combustor. σ is the Stefan-Boltzmann constant = 5.67.10-8 W·m-2·K-4. tp.T(i) behaves as the internal wall temperature, to be determined in the connected wall component. Ti,aw is the adiabatic wall temperature defined as: 0.33 γ i − 1 2 Ti , aw = Ti (1 + Pri , ref M i ); Pri , ref = (Cpi λi / ηi ) ref 2 The reference conditions (ref) are calculated at a temperature halfway between the wall and the free stream static temperature. It’s supposed that the mixture composition do not change between these two temperatures. The film coefficient is calculated using empirical correlations according to Bartz [16]:
either to a complete nozzle or to a nozzle extension. The important point is to know the inlet total conditions calculated from the upstream component:
• If the nozzle is connected to a chamber (case of a complete nozzle), total conditions will be those of the exit of the chamber. • If the nozzle is connected to another nozzle (case of a nozzle extension), total conditions will be those of the exit of the upwind nozzle. a) Throat section calculation. The choked throat conditions are calculated with the following equations assuming known total enthalpy and the static entropy (those of the last combustor volume):
hch ,tot = h( N k ,th , Tth ) + vel th / 2; s ch = s ( N k ,th , Pth , Tth ) 2
vel th = γ th Rth Tth ;
mth = Ath velth Pth /( Rth Tth )
Sub index “ch” denotes the conditions at the exit of the combustor. h, s are the CEA thermodynamic functions. Nk,th is the number of moles of the chemical constituent k at throat. It is calculated from the Minimum Gibbs energy method as for the chamber sections. b) Supersonic sections calculation. It is assumed that the non-adiabatic process takes place into two separated steps: 1) Heat losses. From section i to i+1, the following losses in enthalpy and entropy will take place:
hi +1,tot = hi ,tot − q wall ,i / mt ; s i +1 = s i − q wall ,i / mth / Ti where qwall,i is calculated using Bartz correlations (see above). 2) Expansion. Assuming now that the isentropic relations are valid we have:
hi +1,tot = h( N k ,i +1 , Ti +1 ) + veli +1 / 2 2
si +1 = s( N k ,i +1 , Pi +1 , Ti +1 ); mth = ρ i +1veli +1 Ai +1 Nk,i is the number of moles of the chemical constituent k at the nozzle station i. It is calculated from the Minimum Gibbs energy method. The three equations above with three unknowns (pressure, temperature and speed) are solved iterating in pressure and temperature. Two possibilities are foreseen, equilibrium or frozen flow. In the last case, molar fractions remain constant. c) Thrust calculation. The expressions used to calculate the thrust F and the ISP at section no. i are:
Fi = mth veli + (Pi − Pout )Ai ;
ISPi = Fi / mth
where Pout is the external pressure. Cooling Jackets These components represent Regenerative Circuits of a Chamber including or not inlet/outlet torus. A 3D geometry (built by means of several 3D walls around the channels) will be taken into account:
.2 6 4 h film,i = 0.026 ⋅ η i0,ref ⋅ (λ / η ) i0,.ref ⋅ Cpi0,.ref ⋅ mth0.8 / Ai0.9 ⋅ (πDth / 4 / R curv ) 0.1
1D Nozzle conservation equations A resolution method based on the transient conservation equations would probably had numeric problems (passage from subsonic to supersonic regime, shock waves, etc), so a 1D quasi-steady implicit method has been implemented including nonisentropic effects under frozen or equilibrium conditions. The mathematical model can be applied
Regenerative circuit
Regenerative circuit level
Regenerative circuit section
The cooling jacket is divided into a variable number of sections in axial direction. Each section is made of:
• one fluid node of the “Tube” component (FLUID_FLOW_1D library), which is simulating the cooling channels and, • five slides of the “wall_3D” component simulating the heat conduction in every direction including the axial. They are arranged as follows:
EXTERIOR
w_ch
0.5(
Ext.
π ·D − w _ ch n _ ch
Ext_right nr
nr
ns
th e
ns Right
Channel nr
t ch
nr
th i
ns Int.
• When this backflow is stopped again at the closed end, the “negative” pressures that should be produced are limited to the vapor pressure. Then, a vapor bubble generation takes place, producing a water column separation, the column entering into the tank. • The bubble collapse begins when this liquid column is decelerated by the tank pressure and begins to enter the pipe, collapsing the vapor. A new cycle starts when the liquid column is stopped at the closed end like in a bouncing “ball”. Wall friction progressively smooth water hammer peaks as it is shown in the plots. b) With non-condensable gas. Assuming only a 0.1 % in mass of non-condensable gas (air) not solved in the liquid, the water hammer is greatly reduced:
Int_right nr
ns
ns
CHAMBER
The three wet surfaces in contact with the refrigerant are supposed to be at three different temperatures. The heat flux at the channel walls interfaces is coupled with the “Tube” component simulating the channels, where the heat exchange coefficient is calculated depending on the phase of the refrigerant
APPLICATION EXAMPLES A) Water hammer with cavitation This example shows the ESPSS capabilities concerning water-hammer analysis (hydrodynamic cavitation) in a pipe using real properties fluid (two phase flow). See other more complex validation cases in [17].
The model represents a pipe between two boundaries: The left one imposing P-T (5 bar, 300 K). The right one forces the mass flow circulating through the pipe. The pipe data are: Fluid: water. Diameter: 10 mm. Roughness: 50 μm. Nodes: 25. Length: 2 m. A previous steady state is calculated to determine the pipe conditions at a constant mass flow of 1 kg/s. The transient will be started by a sudden flow stop at the right side at TIME = 0.05. Below main plots obtained: a) Without non condensable gas. The following phenomena can be shown:
B) Pressurisation circuit A model showing a pressurization system by means of electro-valves is presented below. The Tank is simulated with the more complex component capturing vaporization, vertical fluid temperature distribution, etc. The pressurization system consists of a set of three valves. One of them is normally open. The second will open when the Tank pressure goes down a predefined (but variable in time) threshold, and will close when the Tank pressure goes up another predefined threshold. The third one is normally closed.
The tank is pressurized with Hydrogen. The initial condition in the Tank is: 95% full of LH2 with the ullage volume full of Helium. By simplicity, in this example all the thermal ports of the tank are supposed adiabatic. Below main plots obtained:
• The tank pressure evolves (in accordance with the valve activations) between the predefined thresholds.
• Pressure surge is correctly simulated accordingly with the theoretical value of ρcv. Then, the wave trips caused by a sudden stop of the fluid are reflected in the open end as a negative flow.
• The vaporization level produced at the liquid/gas interface is very important in order to know the number of valve activations. It is noted that the initial mixture process of H2 with He at different temperature probably produce the pressure overshot observed at the beginning, just after the first valve actuation.
• In this respect, the volume discretization can have non-negligible influence on the model results: Indeed, the local wall temperature seen by a liquid volume will depend on its size. The consequence is that the liquid evaporation can be affected by the discretization.
While pressurizing a tank with an inflow jet, vaporization results can be better with an only gas volume in which a complete non-condensable gas / vapor mixture is calculated. The evaporation rate is highly dependent on the vapor partial pressure that can be overestimated with a 1D model at the gas volume near to the interface (jet effects are not taken into account). C) Staged engine cycle This model represents a Staged engine cycle type. Input data are fictitious values, the aim this example being just to show the ESPSS capabilities regarding this type of engine. In this model, the default performances maps of the turbines (subsonic) have been used.
nodes and the cooling jackets, it is possible to use Multiplexer/Demultiplexer components allowing to split a chamber thermal port into two or more different CoolingJackets components. The definition of the nozzle and chamber geometry (diameters and dx length as a function of the axial position) needs common data affecting to the Chamber, Nozzle and the Cooling Jacket. ESPSS libraries facilitate this task by using common variables (non dimensional nozzle profiles affecting both the diameter and the axial mesh size). b) Initialization and ignition control. As a general rule, the combustor should be initialized with the external conditions. Under vacuum conditions, the combustor and the regenerative circuit can be initialized with pressures lower than the saturation conditions at the external temperature (a minimum of 100 Pa is advised). If a proper calculation of the start-up is not necessary, use higher values of the initial pressure The valves opening-closing sequences must be treated with care. In principle, if only steady conditions are wanted, it is possible to begin the simulation with all the valves at its final position. Other times it will be necessary to test more physical opening laws. In general, it is advised to:
• Enter convenient values for the minimum/maximum Mixture Ratio and for the liquid vaporization delay ('MR_min. MR_max, tau_v' input data). • Initialize the chamber with the same vapor as that contained in the injector opening first. Combustor chamber can be initialized with only reducer vapor, oxidizer vapor or non-condensable gas. The ignition flag (a time dependent input data) allows the user to coarsely simulate an igniter. By default, its value is one (ignition activated). Nevertheless, it will be more physical to set this variable to zero while no selfignition is foreseen, no ignition order has been commanded or when the MR at injection level is not enough to allow the combustion heat vaporize all the injected liquids. The combustion chamber also includes a solid propellant starter model, the starter gases mass flow and its temperature being time dependant input data. By default, these values are zero (no starter). The molar fractions of the solid propellant gases are input data. The composition of the solid propellant gases must be within a predefined set of chemicals {H2O, CO, CO2, N2, H2, He}, where most of the starter gases are included. Below main plots obtained:
Building the model. The propellant fluids must be defined inserting two “WorkingFLuid” components in the circuits connected to the respective injector ports. (Mono propellant combustors can be modeled plugging one of the injector inlets). a) Geometry arrangement. Preburner and Combustor have been modeled using a 1D discretization. They have been provided with a thermal port which can be connected to a Cooling Jacket or other thermal components. In some cases, where there is not a direct correspondence between the combustor/nozzle
• With the default burning and evaporation times, the chamber temperatures are less than the equilibrium ones, the model being able to simulate the progressive consumption of the liquid oxygen injected. (see differences between section 1 and 2, T[1], T[2] plots) • The model is able to simulate the start-up and the shut-down of the engine, including the relevant phenomena: cavities priming, liquid vaporization, pressure/temperature rise in the chamber during ignition, cavitations in the feeding lines, etc.
• Two steady working points have been simulated by controlling the VPrO valve, the second one leading to an extremely high pump pressure without gaining in thrust, thus the difficult regulation of this kind of cycle engines. • Temperatures can be very low in the chamber before ignition because the feeding tanks are cryogenic. The simplified model for the liquid vaporization inside the combustor prevents the chamber ignition if the injection’s conditions do not have an appropriate mixture ratio, letting the liquid jets leave the chamber without vaporization.
As for the others engine models, the molar fraction of the products entering and exiting the combustion chambers can be found by exploring the results on the monitor tool.
CONCLUSIONS Main targets of the ESPSS V1 project has been successfully obtained during the first phase of the Project. Based on the outcome of the first Phase, two main targets formed the basis of ESPSS-V2 were successfully fulfilled: a) Identify areas of the tool where further refinement is needed. In this respect ESPSS V2 is appreciably more robust and powerful: 1. Simulations of pipe networks are more robust and faster, specially dealing with two-phase (priming), two-fluids systems. 2. Models of tanks (boiling process) and pressurisation systems (real gas properties for the pressuring gas) are improved with faster simulations 3. Combustor models are upgraded in many respects permitting more successful simulations during the ignition and shut down processes b) Validation of the tool against cases for which experimental data exist: 1. Multiple upgrades concerning robustness were needed to validate the code against experimental data (priming cases and HM7B engine). 2. Version 2.0 incorporates specific new user requirements and stabilizes previous versions. After this validation and refinement phase, a third phase is proposed with the following upgrades: 1. Modeling of boiling at internal walls of 1D Tanks 2. Modeling of gas dissolution in Tanks and Pipes, 3. Modeling of surface tension Tanks 4. Possibility of steady and quasi-steady direct calculations dealing with parametric studies.
Nomenclature α
Void fraction
-
β
Blade angle
º
γ
Isentropic coefficient
-
ϕ
Mass flow coefficient
m3
ψ
Head rise coefficient
m2
η
Dynamic viscosity
μ
Poisson’s ratio
λ
Conductivity
Ns m-2 W m-1 K-1
κ
Isothermal compressibility
Pa-1
κ
Wall compressibility
Pa-1
pw red R
Solid propellant gases Reducer propellant rated (nominal)
Π
Pressure ratio
ρ
Density
ξ
Pressure drop coefficient
-
References
τ
Characteristic time
s
1
ω ai A
Angular speed Cp adjusting coefficients Flow area
Cp C+ D eta F G h hc I m ME
specific heat P=cte J kg-1 K-1 Reduced torque m5 Pipe diameter m Efficiency Thrust N Mass flow per unit of area kg s-1m-2 Total specific enthalpy J kg-1 Heat exchange coef. W m-2 K-1 Inertia kg m2 Mass flow rate kg s-1 Young's modulus Pa
MW n N
Molecular weight Pump reduced speed Axial speed coefficient
gr mol-1 -
N Ns P q Q Q+ R r rpm s ST T T TDH u v v, vel V W x
Number of moles Pump specific speed Pressure heat flow Volumetric flow Mass flow coefficient Specific gas constant Radius revolution per minute Entropy Specific torque Torque Temperature Total dynamic head Total specific energy Pump reduced flow Fluid speed Volume Power vapour mass fraction
Pa W m3 s-1 J kg-1 K-1 m
y Index cr, son ch
molar fraction critical conditions Chamber conditions
bu
Burned gases
eq
Equilibrium conditions
th nc oxy
throat non condensable fluid Oxidizer propellant
kg m-3
rad s-1 m2
J kg-1 K-1 Nm K m J kg-1 m s-1 m3 W -
2
3 4
5
6
7 8
9
10
11
12
13
14 15
16
17
REFPROP V7.0. Thermodynamic & Transport Properties. National Institute of Standards & Technology. (NIST) Gordon, S. and B.J. McBride: Computer Program for the Calculation of Complex Chemical Equilibrium Compositions with Applications. NASA Reference Publication 1311 (1994) Reid, R.; Prausnitz,J.M; Poling, B.; The properties of Gases & Liquids, 4th Edition, 1987 C.-A. Schley & J. Kollien. Algorithms for efficient computation of Propulsion Fluid Properties. AIAA 96-0235, January 15-18, 1996 / Reno, NV E. W. Schmidt, Hydrazine and its Derivatives Preparation, Properties, Applications, Second Edition, Vol. 1, J. Wiley & Sons, 2001 Idelchik, I.E.; Handbook of Hydraulic Resistance, 3rd edition, Research Institute for Gas Purification, Moscow, Russia, 1994 Engineering Data Book III. Wolverine Tube. Inc Chen, J.C., 1966, “ A correlation for Boiling Heat Transfer to Saturated Fluid in Convective Flow, “ I.Eng. Chem. Process Des. Dev., 5, pp. 322-329. Dougall, R.S., Rohsenow, W.H.; Film Boiling on the inside of the inside of the Vertical Tubes with Upward Flow of the Fluid at Low Qualities, M.I.T. Report No. 9079.26, (1963) E.F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. 2nd Edition . Springer, Germany, 1999 K.Michalak and C. Ollivier-Gooch. Differentiability of Slope Limiters on Unstructured Grids. Technical report, Department of Mechanical Engineering, University of British. Columbia, Vancouver, BC, V6T 1Z4, Canada, 2006. Fluid Mechanics, Landau and Lifshitz, Pergamon Press, (Chapter VIII Sound, Section 63, Problem 1: Determine the velocity of sound in a nearly homogeneous two-phase system IEC 60534 2 1. Industrial process control valves. Part 2-1 : Flow capacity. Sizing equations for fluid flow under installed conditions. September -1998 E. Benjamin Wylie & Victor Streeter; Fluid Transients in Systems. Prentice Hall, Inc. 1993 http: //courses.washington.edu/mengr524/Handouts/gib bs.pdf Bartz, D.R., Turbulent Boundary-Layer Heat Transfer from Rapidly Accelerating Flow of Rocket Combustion Gases and of Heated Air, Advances in Heat Transfer, pp. 2-108, 1965 EUCASS 2009 paper 273. A Satellite Platform Modelling with EcosimPro: Simulation and Ground Tests Comparison. Versailles 7th June, 2009, France