ISSIM Spark Ignition

The imposed stretch spark ignition model (ISSIM) is a Eulerian spark-ignition model for 3D RANS simulations of internal combustion engines.

The spherical kernel equation is used directly in the Eulerian transport equations—allowing a description of all local phenomena. The use of a flame surface density equation allows the simulation of multi-spark ignition and of the flame holder effect. Both flame growth and wrinkling are modelled via the ECFM equation.

Assuming an inductive ignition system, the ISSIM model provides the amount of energy transferred to the gas during the glow phase as well as the evolution of the spark, which is then convected and wrinkled by the flow. At the instant of electrical breakdown, an initial burnt gas profile is created on the 3D CFD mesh. From then on, the reaction rate is directly controlled by the flame surface density (FSD) equation. A number of adjustments are then made to this equation so that it can correctly describe the combustion evolution during the early ignition stage.

Electric Circuit

The energy that is required to successfully initiate combustion—which is initially stored as electrical energy in the ignition system—is converted to thermal energy when it is released into the spark gap. The evolution of the electrical energy

E_{s}

in the secondary circuit is given as:

图 1. EQUATION_DISPLAY

\frac{d E_{s} (t)}{d t} = - R_{s} i_{s}^{2} (t) - V_{s p k} i_{s} (t)

(3963)

where

V_{s p k}

is the inter-electrode voltage,

i_{s}

is the current, and

R_{s}

is the resistance.

Not all of the electrical energy is transferred to the spark—a substantial amount is lost by the Joule effect. The intensity in the secondary circuit is estimated by:

图 2. EQUATION_DISPLAY

i_{s} (t) = \sqrt{\frac{2 E_{s} (t)}{L_{s}}}

(3964)

where

L_{s}

is the circuit inductance.

At breakdown, the spark length is equal to the spark gap

d_{i e}

. The spark is then stretched by convection

l_{s p k}^{m e a n}

and by the turbulent motion

Ξ_{s p k}^{t u r b}

. The total length of the spark is given by:

l_{s p k} = Ξ_{s p k}^{t u r b} l_{s p k}^{m e a n}

(3965)

The spark wrinkling evolution equation of

Ξ_{s p k}^{t u r b}

takes the form:

图 3. EQUATION_DISPLAY

\frac{1}{Ξ_{s p k}^{t u r b}} \frac{d Ξ_{s p k}^{t u r b}}{d t} = \frac{1}{2} (K_{s p k, T} + K_{s p k, M})

(3966)

where

K_{s p k, T}

and

K_{s p k, M}

are the spark strains—that are caused by the turbulent flow, and the mean flow, respectively.

K_{s p k, T}

corresponds to the effect of the turbulent eddies that are greater than the arc thickness, but less than half of the length of the spark. The computation of

K_{s p k, T}

is similar to the ITNFS function (net flame stretch function) in the equation for the flame surface density.

During the arc and glow phases, only a fraction of the spark energy is released to the gas. The energy that is released in a thin region near the electrodes is essentially lost due to a fall in the voltage. The voltage distribution during the arc and glow phases is represented by the following diagram:

The electric potential difference between the electrodes (the spark voltage) is given by:

图 4. EQUATION_DISPLAY

V_{s p k} = V_{c f} + V_{a f} + V_{i o n} + V_{g c}

(3967)

where

V_{c f}

V_{a f}

V_{i o n}

and

V_{g c}

are the cathode voltage drop, anode voltage drop, the ionic potential of the cathode, and the gas column voltage, respectively. The anode voltage drop is similar for both the arc and glow modes.

For Inconel (an alloy based on nickel, chrome, and iron), $V_{i o n}$ is 7.6V. During the arc phase, $V_{c f}$ is included in $V_{i o n}$ , and is equal to 252 V during the glow phase.

The gas column voltage is given by:

图 5. EQUATION_DISPLAY

V_{g c} = 40.46 l_{s p k} i_{s}^{- 0.32} p^{0.51}

(3968)

The breakdown and arc phases are complex due to the extremely short duration, presence of a plasma channel, very high temperatures, and unsteady behaviour. At breakdown, about 60% of breakdown energy (which is a function of the spark gap

d_{i e}

and the breakdown voltage

V_{b d}

) is released to the gas, and provides the ignition energy

E_{i g n}

that is required for the glow phase. The breakdown energy is expressed as:

图 6. EQUATION_DISPLAY

E_{b d} = \frac{C_{b d}}{d_{i e}} {(\frac{V_{b d}}{K})}^{2}

(3969)

where $C_{b d}$ and $K$ are constants.

The breakdown voltage

V_{b d}

is expressed as:

图 7. EQUATION_DISPLAY

V_{b d} = K_{b d} F_{c}^{1 / 2} K (\frac{ρ}{ρ_{0}}) d_{i e}

(3970)

where

ρ_{0}

is the density at standard conditions,

K_{b d}

is a user adjustment parameter for the first breakdown voltage evaluation,

K = 1.5 e^{6} v o l t . m^{- 1}

, and

F_{c}

is a correction function—which is dependant on the thermodynamic conditions and is different for the first breakdown and secondary breakdown (restrike).

First Breakdown

图 8. EQUATION_DISPLAY

F_{c} = \frac{\max (16, \frac{ρ}{ρ_{0}})}{\max (0.1, \frac{ρ}{ρ_{0}})}

(3971)

Secondary Breakdown

No Correction: 图 9. EQUATION_DISPLAY

$F_{c} = 1$

(3972)
Corr. Factor: 图 10. EQUATION_DISPLAY

$F_{c} = K_{b d 2}$

(3973)
Corr. Factor, Rho_g and EGR Mod: 图 11. EQUATION_DISPLAY

$F_{c} = K_{b d 2} \cdot f (Y_{e g r s p k})$

(3974)
Corr. Factor, Rho_plasma Mod: 图 12. EQUATION_DISPLAY

$F_{c} = K_{b d 2} \frac{\max (16, \frac{ρ_{p}}{ρ_{0}})}{\max (0.01, \frac{ρ_{p}}{ρ_{0}}) 1600}$

(3975)

K_{b d 2}

is a user adjustment parameter for the second (restrike) breakdown voltage evaluation,

Y_{e g r s p k}

is the EGR mass fraction near the spark, and

ρ_{p}

is the plasma density which is determined by:

图 13. EQUATION_DISPLAY

ρ_{p} = ρ \frac{T}{T + Δ T}

(3976)

where

Δ T

is the increase in temperature due to the spark energy.

During the glow phase, the voltage drop is localised in the vicinity of the electrodes. It is therefore assumed that the energy that is released within this region is lost to the electrodes. Finally, the energy that is transferred to the gas is deduced from the gas column voltage and current by:

图 14. EQUATION_DISPLAY

\frac{d E_{i g n} (t)}{d t} = V_{g c} (t) i_{s} (t)

(3977)

ISSIM also accounts for the turbulent dilution of the energy release by the spark outside its influence zone with the following equation:

图 15. EQUATION_DISPLAY

E_{i g n} = E_{i g n} (1 - 2 K_{d i l} \frac{d t}{τ})

(3978)

where

K_{d i l}

is a user constant, and

τ

is a spark turbulent mixing timescale.

E_{i g n}

contributes to locally increasing the gas temperature around the spark, and is distributed on the volume mesh according to a Gaussian pdf, which adapts itself in one of two ways:

height: ${S p}_{N R}$ has a fixed value of 2.
width: ${S p}_{N R}$ varies to respect an acceptable increase of temperature.

图 16. EQUATION_DISPLAY

E_{s p k}^{g a s} = E_{i g n} \frac{\exp [- {(\frac{| x_{c} - x_{s p k} |}{0.007 {S p}_{N R}})}^{2}] v o l}{\int \exp [- {(\frac{| x_{c} - x_{s p k} |}{0.007 {S p}_{N R}})}^{2}] v o l}

(3979)

where $x_{c}$ denotes the cell centroid and $x_{s p k}$ is the spark location coordinates.

E_{i g n}

is used to decide if ignition has occurred or not. It is also assumed that a flame kernel forms around the spark, and ignition starts only if

E_{i g n} (t) > E_{c r i t} (t)

within a minimum ignitable volume fraction around the spark

V_{c r l}

E_{c r i t}

is a critical ignition energy that depends on thermodynamic properties in the vicinity of the spark:

图 17. EQUATION_DISPLAY

E_{c r i t} = \frac{γ}{γ - 1} l_{c} 4 π p δ_{L}^{2}

(3980)

where

l_{c}

can be defined as

l_{s p k}

d_{i e}

δ_{L}

is the flame thickness,

p

is the absolute pressure, and

γ

is the ratio of the specific heats.

Under these conditions, an amount of burnt gas

m_{b}^{i g n}

is produced at the spark, equivalent to the mass in a cylinder of radius

2 δ_{L}

and height

l_{s p k}

. There are three options available to evaluate

m_{b}^{i g n}

Minimum Critical Burning Mass: $m_{b}^{ign} = K_{mbign} 〈 ρ_{u} l_{spk} 4 π {δ_{L}}^{2} 〉 spkplg$

(3981)
From Actual Spark Energy Released: $m_{b}^{ign} = K_{mbign} 〈 \frac{ρ_{u} E_{ign} (γ - 1)}{γp} 〉 spkplg$

(3982)
From Fixed Plasma Temperature: $m_{b}^{ign} = K_{mbign} 〈 Max (ρ_{u} l_{spk} 4 π {δ_{L}}^{2}, \frac{E_{ign}}{T_{plasma} C_{p}^{gb}}) 〉 spkplg$

(3983)

where:

$K_{m b i g n}$ is a user specified initial ignition burnt mass coefficient
$ρ_{u}$ is the density of the unburnt gas
$T_{p l a s m a}$ is an assumed temperature of plasma of 10000 K
$C_{p}^{g b}$ is the specific heat of the burnt gases at constant pressure
$< > spkplg$ means average over a sphere that is centered on the spark plug location and weighted by a Gaussian type probability density function (PDF) in the radial direction of the sphere

Flame Kernel Growth Using ISSIM

If the critical energy criterion is satisfied, an initial burnt gas mass is produced. This is evaluated using a target burnt gas volume fraction

{\bar{c}}_{i g n}

defined by:

图 19. EQUATION_DISPLAY

{\bar{c}}_{i g n} = c_{0} \exp [- {(\frac{| x - x_{s p k} |}{{S p}_{C} d_{i e}})}^{2}]

(3985)

where

d_{i e}

is the spark gap,

x - x_{s p k}

the distance away from the spark plug location

x_{s p k}

, and the constants

c_{0}

and

{S p}_{C}

satisfy the conditions:

图 20. EQUATION_DISPLAY

\int \bar{ρ_{b}} {\bar{c}}_{i g n} d V = m_{b}^{i g n}

(3986)

and

${S p}_{C} = 0.7$ or $c_{0} = 1$

In order to assign a value of

{\bar{c}}_{i g n} (x, t)

in the 3D computational domain, the reaction rate of the progress variable

c

is written as:

图 21. EQUATION_DISPLAY

\bar{{\dot{ω}}_{c}} = \max (ρ_{u} S_{L} Σ, \frac{ρ_{b} ({\bar{c}}_{i g n} - \bar{c})}{d t})

(3987)

where

S_{L}

is the laminar flame speed,

Σ

the flame surface density,

ρ_{b}

and

ρ_{u}

are the density of the burnt and unburnt gas, respectively. In the above expression, the flame surface density (FSD) of a possible pre-existing flame is also taken into account. This allows the consideration of the interaction between flames coming from different spark ignitions.

The main purpose of ISSIM is to model the reaction rate of a growing flame kernel via the FSD equation starting from the very beginning of spark ignition. As a result, the equation needs to be initialised as soon as the initial burnt gas kernel of radius

r_{b}^{i g n}

is created. The following expression is proposed for this purpose:

图 22. EQUATION_DISPLAY

\bar{Σ} = \max (\bar{Σ}, \frac{3}{r_{b}^{i g n}} {\bar{c}}_{i g n})

(3988)

The term

\frac{3}{r_{b}^{i g n}} {\bar{c}}_{i g n}

is chosen in order to recover the initial total flame surface of a sphere of radius

r_{b}^{i g n}

so that:

图 23. EQUATION_DISPLAY

\int Σ d V = \int \frac{3}{r_{b}^{i g n}} {\bar{c}}_{i g n} d V = \frac{3}{r_{b}^{i g n}} \frac{4}{3} π {(r_{b}^{i g n})}^{2}

(3989)

It is also possible to provide ignition source terms after the initial burnt mass deposit from the resolution of the OD equation of the flame kernel radius, instead of using the source terms from the flame surface density equation. In this case, flame surface density (if available) is only initialized after transition from ignition to combustion if the flame surface density equation is solved. This option is mostly used with the complex chemistry model which does not evaluate flame surface density.

OD Flame Kernel Radius Equation

It is possible to solve the flame kernel radius from a 0D equation:

图 24. EQUATION_DISPLAY

\frac{d R_{k}}{d t} = \frac{ρ u}{ρ b} S_{T}

(3990)

where

ρ u

and

ρ b

are the unburnt density and burnt density respectively, and

S_{T}

is a dynamic turbulent flame speed given by:

图 25. EQUATION_DISPLAY

S_{T} = S_{L}^{M} \cdot X i

(3991)

where

S_{L}^{M}

is a modified effective laminar flame speed that is due to the influence of energy that is released by a spark and the curvature of the flame kernel.

X i

is a dynamic turbulent wrinkling factor that is given by:

图 26. EQUATION_DISPLAY

{X i}^{n + 1} = \min ({X i}_{\max}; 1 + [{X i}^{n} (1 + \Pr o d . d t) - 1] \frac{S_{k}^{n}}{S_{k}^{n + 1}})

(3992)

where

S_{k}^{n}

and

S_{k}^{n + 1}

denote the kernel radius surface at time step

n

and

n + 1

respectively.

\Pr o d

is the turbulent wrinkling production source, given by:

图 27. EQUATION_DISPLAY

\Pr o d = K_{x p} G (u^{'}, L_{int}, L_{k}, R_{k}, δ_{L})

(3993)

where

K_{x p}

is an adjustment constant, and

G

is a function of

u^{'}

L_{int}

L_{k}

R_{k}

, and

δ_{L}

, which denote the turbulent fluctuation velocity, turbulent integral length scale, Kolmogorov turbulent length scale, kernel radius, and laminar flame thickness, respectively.

The fully turbulent wrinkling factor

{X i}_{\max}

is given by:

图 28. EQUATION_DISPLAY

{X i}_{\max} = 1 + K_{x m} F (\frac{u^{'}}{S_{L}^{e f f}})

(3994)

where

K_{x m}

is an adjustment constant and

F

is a correlation function of

\frac{u^{'}}{S_{L}^{e f f}}

$S_{L}^{e f f}$ is the effective laminar flame speed.

Combustion heat release and associated reactions are fully driven by this flame kernel radius equation in the influence zone of the spark up to the transition to fully turbulent combustion that is driven by the combustion model itself.

Transition happens when the flame kernel radius

R_{k}

is bigger that the kernel radius transition

R_{k}^{t r}

, which is given by:

图 29. EQUATION_DISPLAY

R_{k}^{t r} = \min (R_{k}^{L}, \max (5 R_{k}^{i n i}, T_{r} L_{t}))

(3995)

where

R_{k}^{L}

is a maximum transition radius constant,

R_{k}^{i n i}

is the initial kernel radius given by the initial burnt mass deposit,

T_{r}

is a transition ratio constant, and

L_{t}

is the turbulent length scale.

The evolution of this kernel radius allows the combustion progress variable to be updated on the 3D mesh from the same distribution function that is used in Eqn. (3985), and burnt mass conservation, providing the species source terms with:

图 30. EQUATION_DISPLAY

\frac{d \tilde{c}}{d t} = \max (0, \frac{{\tilde{c}}_{i g n} - \tilde{c}}{d t})

(3996)

where

\tilde{c}

is the existing combustion progress variable, and

{\tilde{c}}_{i g n i}

is the updated combustion progress variable from ignition.

Modification of the FSD equation: Source terms at the spark plug in the presence of a spark

As long as the spark exists (

E_{s} > 0

), the ignition criterion defined by

E_{i g n} > E_{c r i t}

is calculated at the spark plug. Once this criterion is satisfied, we assume that a burnt gas kernel corresponding to the ignition burnt gas mass

m_{b}^{i g n}

exists at the spark plug. Without convection, the burnt gas mass is only produced at the spark timing interval and this condition is always satisfied. If the flame is convected away from the spark, it is replaced by incoming fresh gases which ignites in the vicinity of the spark. If the ignition criterion is satisfied, a target profile

{\bar{c}}_{i g n} (x, t)

is imposed at any instant. In the RANS approach, the profile of

{\bar{c}}_{i g n} (x, t)

diffuses rapidly so that

\bar{c} (x, t) \leq {\bar{c}}_{i g n} (x, t)

. As a result, the burnt gas mass that is produced in the vicinity of the spark is overestimated. The following diagram shows a one-dimensional illustration of the progress variable evolution at spark time:

A solution to this problem is to determine the right amount of burnt gas mass

d m_{b}^{i g n}

that is produced around the spark. This amount, which should be non-zero after spark ignition only if the flame kernel is convected, can be estimated using a phenomenological approach as illustrated in the diagram below which shows the burnt gas mass that is produced at the spark upon flame kernel convection:

The mass

d m_{b}^{i g n}

corresponds approximately to the mass of a parallelepiped of height equal to the spark gap (

d_{i e}

), depth

2 δ_{L}

, and width

δ = d + r_{b}^{i g n} - r_{b} (t)

and is given by:

图 31. EQUATION_DISPLAY

d m_{b}^{i g n} = {\bar{ρ}}_{u} 2 δ_{L} d_{i e} δ

(3997)

where

r_{b} (t)

is the result of the evolution of the flame kernel radius as a function of time. As the spark life is very short, we neglect the flame wrinkling and the pressure variation effect. Several tests have shown that these assumptions are reasonable. As for the initial burnt gas mass, the target burnt gas volume fraction

{\bar{c}}_{i g n}

is defined by requiring that the constant

c_{0}

satisfies the condition:

图 32. EQUATION_DISPLAY

d m_{b}^{i g n} = \int {\bar{ρ}}_{b} \max (({\bar{c}}_{i g n} - \bar{c}), 0) d V

(3998)

The spark source term

S_{c o n v}

in Eqn. (3909) for the ECFM model is expressed as:

图 33. EQUATION_DISPLAY

Φ (\bar{Σ}) = \max (\frac{3}{r_{b}^{i g n}} {\bar{c}}_{i g n} - \bar{Σ}, 0.0)

(3999)

A consequence of this modelling choice is that if ignition is successful at breakdown time and thereafter, the flame holder effect can be retrieved.

Modification of the FSD equation: Modification of the laminar flame speed near a spark plug

The energy

E_{i g n}

transferred to the gas column leads to an increase of the burnt gas temperature in the flame kernel. This increase in turn leads to an increase in the laminar flame speed, which in the vicinity of the spark can be expressed as:

图 34. EQUATION_DISPLAY

S_{L}^{e f f} = S_{L}^{0} [1 + K_{s l} (S_{L}^{c o r r} - 1) \exp (- {(\frac{2 | x - x_{s p k} |}{{S p}_{S L} l_{s p k}})}^{2}]

(4000)

where

S_{L}^{0}

is the standard laminar flame speed,

K_{s l}

and

{S p}_{S L}

are model constants.

图 35. EQUATION_DISPLAY

S_{L}^{c o r r} = S_{L}^{0} 0.5 (1 - \frac{4 δ_{L}^{0}}{l_{s p k}} + \sqrt{Δ})

(4001)

and

图 36. EQUATION_DISPLAY

Δ = {(\frac{4 δ_{L}^{0}}{l_{s p k}} - 1)}^{2} + 4 (\frac{4 δ_{L}^{0}}{l_{s p k}}) (1 + \frac{Δ T_{b}^{i g n}}{400})

(4002)

where

δ_{L}^{0}

is the flame thickness.

图 37. EQUATION_DISPLAY

Δ T_{b}^{i g n} = \min (\frac{E_{i g n}}{C_{p}^{g b} m_{b}}, T_{u}^{\max})

(4003)

where

T_{u}^{\max}

is a model constant,

c_{p}^{g b}

is the specific heat capacity of the burnt gases, and

m_{b}

represents the total burnt gas mass, given by:

图 38. EQUATION_DISPLAY

m_{b} = {\bar{ρ}}_{b} \frac{1}{6} π d_{i e}^{3}

(4004)