Engineering Applications of Computational Fluid Mechanics Vol. 1, No.3, pp.

207215 (2007)

NUMERICAL ANALYSIS OF THE FLOWFIELD IN A SUPERSONIC COIL

WITH AN INTERLEAVED JET CONFIGURATION
AND ITS EFFECT ON THE GAIN DISTRIBUTION
Zong-min Hu*#, Rho-shin Myong*+, Anh-thi Nguyen*, Zong-lin Jiang# and Tae-hwan Cho*
* Research Center for Aircraft Parts Technology and School of Mechanical and Aerospace Engineering,
Gyeongsang National University, Jinju 660-701, South Korea
#
LHD Laboratory, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100080, China
+
E-Mail: [email protected] (Corresponding Author)
ABSTRACT: In a supersonic chemical oxygen-iodine laser (COIL) operating without primary buffer gas, the features
of flowfield have significant effects on the Laser efficiency and beam quality. In this paper three-dimensional, multispecies, chemically reactive CFD technology was used to study the flowfield in mixing nozzle implemented with a
supersonic interleaving jet configuration. The features of the flowfield as well as its effect on the spatial distribution of
small signal gain were analyzed.
Keywords:

COIL without buffer, interleaved jet, flowfield, small signal gain

In this paper, a new jet configuration called

transverse interleaving jet was proposed and applied
to a primary-buffer-free COIL. Then, threedimensional (3D) CFD technology was used to
evaluate the flowfield and its effect on the gain field
in the mixing nozzle.

1. INTRODUCTION
Chemical oxygen-iodine laser (COIL) is a highpower chemical laser device that operates on a near
infrared radiation ( = 1.315 m) of atomic iodine.
The laser medium, I*, is produced through a
pumping reaction between singlet oxygen, O2(1),
and ground state iodine, I. Research on COILs has
accelerated owing to the potential military and
industrial applications. Recently, a new concept of
COIL operating without primary buffer gas was
proved to be of higher efficiency and stability, and
more applicability for mobile integration. This kind
of COIL uses nitrogen rather than helium as the
buffer gas in the secondary flow (Furman,
Barmashenko and Rosenwaks, 1998; Furman et al.,
2001; Fang, Sang and Chen, 2002 & 2003), which
results in characteristic velocity reduction. A series
of experimental studies were completed with a
supersonic COIL at Ben-Gurion University
(Furman, Barmashenko and Rosenwaks, 1998;
Furman et al., 2001; Bruins et al., 2002; Rybalkin
et al., 2002). They concluded that the gain
distribution is often inhomogeneous because of
inefficient mixing (Rybalkin et al., 2002).
Transverse jets in the transonic or supersonic
section chock the nozzle flow and induce
complicated flow discontinuities, which may in turn
significantly affect the mixing and pumping
processes.

2. NUMERICAL MODELS AND

ALGORITHMS
2.1

Governing equations and numerical models

The COIL flow with a low density and Reynolds

number from 102 (based on the jet flow conditions
and the orifice) to 103 (based on the primary flow
conditions and the inlet chamber size) is supposed
to be governed by 3D Navier-Stokes equations. In
Cartesian coordinates, the N-S equation for multicomponent, chemically reactive system is written
as:

U Fi Di
+
=
+S
t xi xi

(i = 1,2,3)

(1)

Here, U is the unknown variable vector

Fi ( i = 1,2,3 ) are convective flux vectors in xi (x, y
and z) directions, Di are the dissipation flux
vectors, and S is the chemical reaction source term,
where

Received: 9 Jan. 2007; Revised: 13 Mar. 2007; Accepted: 28 Apr. 2007

207

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

U = [1 2 L ns u1 u2 u3 e]

(1.1)

Fi = [1 ui 2ui ns ui uiu1 + i ,1 p

uiu2 + i , 2 p uiu3 + i ,3 p (e + p )ui ]

Di = 1v D1,i

h sp

S = [&1 & 2 L & ns 0 0 0 0]

(1.4)

S sp

sp , and the partial density of the

species sp ( sp = 1,2, L ns ) are denoted by sp .

a 2 sp
T

+ a 3 sp + a 4 sp T

a1sp
2T

T3 +

+ a 2 sp

a 7 sp

T3 +

a 6 sp

a1sp

a 7 sp
4

(2.1)

T3

a 4 sp
a 5 sp 2
ln T
+ a 3 sp +
T+
T
T
2
3

T4

a 2 sp

a8 sp

a 8 sp
2T 3

a 9 sp

+ a3 sp ln T + a 4 sp T +

T4 +

a 8sp
3T 3

(2.2)

T
a5 sp
2

+ a 10 sp

T2
(2.3)

In equations (2.1) to (2.3), parameters a1 to a10 are

coefficients for calculating thermodynamic
properties of individual species and can be found in
the work of Paschkewitz et al.(2000). Therefore, the
enthalpy of the mixture, h , is

u i ( i = 1,2,3 ) are velocity components in x, y and z

directions. e and p denote the total energy and
the pressure of the mixture respectively. v Dsp ,i
represents the diffusion velocity of the species sp
in xi direction. and q define the viscous stress
and heat conduction terms, and & sp is the variation
rate of the species sp . The transport properties of

ns

h = csp hsp

(2.4)

sp =1

where c sp = sp is the mass fraction of species

each pure component, such as viscosity and thermal

conductivity, and thermodynamic properties, such
as special heat capacity and enthalpies, are obtained
from relations given in the power series form of
temperature. The effective binary diffusion model
rather than Ficks model is applied to approximate
the diffusion velocities of all species. For the
chemical kinetic model, 21 elementary reactions
and 10 species, I2, I2*, I, I*, O2(1 ), O2(1),
O2(3), H2O, cl2 and N2 are taken into
consideration. The transport property model,
thermodynamics definitions, diffusion model and
chemical reactions were synthetically analyzed and
presented (Paschkewitz et al., 2000; Eppard et al.,
2000).
The total energy per specific volume is defined as:

e = h p + ui ui / 2, i = 1,2,3

a 6 sp

Rsp

In the above equations, the total density of the

ns

R sp T

(1.3)

(1.2)

2 v D 2,i ns v Dns ,i

+ a 5 sp T 2 + a 6 sp T 3 + a 7 sp T 4 +

i ,1 i , 2 i ,3 u j i , j + qi ], j = 1,2,3

mixture =

a1sp

Cp sp R sp =

sp . According to Dalton's law, pressure p is

ns

p = sp RspT

(2.5)

sp =1

where, R is the gas constant and T is the mixture

temperature.
The thermal conductivity coefficient and the
viscosity of each pure species can be approximated
by the following polynomials:

= a + bT + cT 2 + dT 3

(3.1)

= A + BT + CT 2

(3.2)

The coefficients, a , b , c , d , A , B , and C , for

calculating transport properties in equations (3.1)
and (3.2) can be found in a relevant literature
(Paschkewitz et al., 2000). And the transport
coefficient, = ( , ) of the mixture is given by
Wilkes equations:

(2)

The specific heat, enthalpy and entropy for each

species can be calculated by the following
polynomials:

208

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

i i

ns

mix =
i =1

ns

0.5

0.25 2

ssp
hsp

Rsp RspT

wi
81 +

w j

0.5

''
'
( sp
,r sp , r )
Patm sp=1

RT

(4.1)

(6.3)

where, the molecular weight and the molar fraction

of each species are denoted by w and .
An accurate description of the molecular diffusion
of chemical species is important in the low pressure
flowfields of chemical lasers. The RamshawDukowicz approximation for multi-component
diffusive transport is used, including concentrationgradient-driven
and
pressure-gradient-driven
contributions:

sp vsp = t Gspx

ns

w
ij = 1 + i j w

j
i

ns

sp =1

j ij

j =1

e, r = exp ( sp'' , r sp' , r )

(4)

'
In these equations, sp ,r

Gspp = wsp Dsp ( sp sp )P P

are the

stoichiometric coefficient of species sp in the r th

elementary reaction as the reactant and the product,
respectively. The pre-exponential factor, Cr , the
temperature index, nr , and the activation energy,
Ear are given by the detailed finite-rate chemical
reaction model. The main reaction kinetic processes
can be illustrated through Fig. 1. Reaction equations
(1) to (4) model the self-catalyzed chain reaction of
the dissociation of I2, where equation (1) is the
primary onset reaction and equation (3) represents
the resonant energy transfer or pumping reaction.
Equations (5), (6) and (7) denote different sorts of
quench processes by third-bodies. Equation (8) is
the stimulated laser radiation. For simplicity, other
reaction equations are not included in Fig.1 but can
be found in the work of Paschkewitz et al. (2000).

p sp ns p
sp ns x

G G
j =1 j t sp j =1 j

Gspx = wsp Dsp sp

"
and sp ,r

(5)
(5.1)
(5.2)

Here, t is the total molar concentration. The

average diffusion coefficient of species Dsp is
calculated by a binary effective diffusion model.
The production rate of species sp is determined by
NR

"
& sp = wsp sp' , r sp
,r
r =1

k f ,r

ns

( )
j

j =1

'
j ,r

kb , r

ns

( )
j

"
j ,r

j =1

(6)

The forward reaction rate constant, k f , r , of each

reaction is calculated by the Arrhenius equation and
the corresponding backward reaction rate constant,
k b , r , can be derived from the equilibrium
coefficient, ke , r , and the forwards rate constant:

k f ,r = C r T nr exp( Ear / RT )

b,r =

f ,r
e,r

Fig. 1

(6.1)
2.2

(6.2)

Sketch of the kinetic processes in the COIL

mixing nozzle.

Numerical algorithms

The semi-discretized
equation (1) is

difference

U
n
n
n

= CONVi , j, k + VISC i , j, k + Si , j, k
t
209

equation

for
(7)

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

The dissipative term, VISC, is discretized by a

second-order centered difference scheme, while the
convective term, CONV, is discretized by using a
second-order scheme, DCD (Dispersion-controlled
dissipative) scheme (Jiang, Takayama and Chen,
1995; Jiang, 2004; Hu et al., 2004) proposed by
Jiang et al. for shock-wave-capturing. The
computational domains used in present simulations
were all discretized on body-fitted, structured grids,
with local mesh refinement near the jet regions and
the solid boundaries.
DCD scheme uses the Steger-Warming (Steger and
Warming, 1981) flux splitting algorithm. In the
general coordinate system ( , , ) , the split
convective flux vectors are simple and symmetrical.
Here is the split convective flux vector of
~
flux F in direction:

a month on a self-parallelized 64-bit workstation

with 4-processors.
3. NUMERICAL RESULTS AND
DISCUSSIONS
3.1

To verify the numerical algorithms and validate the

numerical solutions, a cold flow with a transonic
injection in a RADICL (research assessment device
improvement chemical laser) slit nozzle has been
studied. The geometry used in the simulation was
based on the experimental device by Miller et al.
(2000). The complete RADICL nozzle has 115
large injectors and 230 small injectors on either side
of the nozzle blade, and the nominal nozzle width is
25.4 cm. The computational domain used the unit
cell approximation based on the symmetry of the
RADICL nozzle geometry. The symmetry planes
include the nozzle centerline. The centerline passes
through the large injector and the symmetry line
bisects the distance between two adjacent small
injectors. In other words, the computational domain
is a 1/460 cut of the entire nozzle flowfield.
Fig. 2(a) shows the computational domain and the
enlargement of the injectors of two symmetric unit
cells. A coarse grid (601 71 31) and a fine
grid (601 91 41) were used to discretize the
nozzle. The flow conditions were set according to
the experimental setup in the study of Miller et al.
(2000). Firstly, a secondary flow was
perpendicularly injected into the primary flow near
the nozzle throat. Secondly, helium was used for
both the primary and secondary flows. Then we
used a total temperature and pressure of Tp = 293 K
and pp = 7955 Pa for the primary flow (p) and
Ts = 403 K and ps = 32925 Pa for the secondary
flow (s).
Fig. 2(b) shows computed wall pressure together
with the measured experimental data. There is a
good agreement between the numerical results for
both the coarse and fine grids and the experiment.
The discrepancy in the pressure measurements is
highest near the downstream nozzle exit, which is
considered to be a consequence of injecting a purge
flow in the z direction downstream of the nozzle
throat in the experiment (Miller et al., 2000). The
application of the fine grid causes a minor change
in the wall pressure. From a direct comparison of
these results, it can be concluded that the numerical
solutions are well validated for the jet flow in the
COIL.

~
) ~ ~

C1 2( 1)1 + ns +3 + ns +4

~
) ~ ~

Cns 2( 1)1 + ns+3 + ns+4

~
~
~
) ~
F = ) u 2( 1)1 + (u ckx )ns+3 + (u + ckx )ns+4
2J
~
~
) ~
v 2( 1)1 + v cky ns +3 + v + cky ns +4

~
~
) ~
w 2( 1)1 + (w ckz )ns +3 + (w + ckz )ns +4
~
~
)
2 ~
2 ( 1)H c 1 + (H c )ns+3 + (H + c )ns+4
(8)

[
[
[

]
](
]
]

1
where i = (i i2 + 2 )
2

i = 1,2,L ns + 4 ,

and is very small number to keep all the

eigenvalues nonzero. = uk x + vk y + wk z

k x = x k y = y kz = z
= x2 + y2 + z2

Here,

x , y , z

are

transformation matrices between the physical and

computational space, and J is the Jacobian of the
transformation.
An explicit and second-order Runge-Kutta
integration method is used to advance the
computational solution in time toward steady-state.
The convergence criterion is defined as

Re s =

Validation

m& in m& out

10 4 , where m& in and m& out are
m& in

the mass flow rates of some important species at the

inflow and outflow surface of the nozzle
respectively, for example, helium for the validation
case and oxygen for other cases. One run took about

210

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

3.2 Computational setup for the transverse

interleaved jet

Fig. 2

Fig. 3

Downstream from the throat, two rows of jet

orifices are located on either side of the nozzle
blade. The jet orifices, which have a diameter of
1 mm, are interleaved and set 2 mm apart from each
other in the z direction. The height of the nozzle
throat is 6 mm. The angle between the secondary
jets and the nozzle axis is 45o. Fig. 3(a) shows the
computational domainthe unit cell approximation
enveloped by the thick lines, which include a half
jet orifice on either side of the nozzle blade.
Fig. 3(b) and (c) demonstrate the jet arrangement.
Total grid point of 501 121 31 is used to
discretize the domain, while mesh refinement is
handled near the jets and the nozzle walls. The
secondary flow is diluted by the buffer gas N2. Two
cases were considered, each with a different molar
ratio of I 2 /N 2 : namely case (1) 1/60 and case (2)
2/59. The total flow rate of the secondary jets for
the computational unit, n& s , is 1.86 10-5 mole/s,
and the stagnation temperature, Ts0, is 383 K. For
both cases considered, the molar ratio of the
primary flow, O2(3)/O2(1) /H2O/Cl2, is given as
1/1.5/0.278/0.278. The primary flow rate for the
computational unit, n& p , is 0.99 10-5 mole/s,

RADICL slit nozzle flow simulation for

validation and grid-independence analysis:
(a) computational domain and (b) comparison
between the numerical results and experimental
data in wall pressure.

while Tp0 = 271.3 K.

(a) Computational domain for the interleaved jet configuration in the supersonic section of the COIL nozzle;
(b) Enlargement of the jet region; (c) Arrangement of the jet orifices (length unit: cm).
211

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

3.3

shows the sharp velocity gradient in the near filed

of the interleaved jets, which facilities the mixing
between the primary and jet flows through shearing
entrainment. Near the upper wall, a recirculation is
observed. Between the lower wall surface and the
outmost trajectory of the upper jet, there is also a
region of reverse flow. The spanwise and
streamwise distribution of temperature is given in
Fig. 4(c). It appears relative homogeneous in
streamwise but heterogeneous in spanwise
direction. This distribution of temperature is also
consistent with the experimental measurement
documented by Rybalkin et al. (2002). As the
generation of laser medium is sensitive to the
temperature, such temperature profile will affect the
gain distribution. The flow Mach number given in
Fig. 4(d) indicates that a supersonic flowfield has
been set up within the core flow of the nozzle.

Features of the flowfield

For the jet mixing in a COIL systems operating

without primary buffer gas, the under-expanded
secondary flows are always injected from the
orifices located in the transonic region downstream
from the nozzle throat. Generally, the injection of
gas into the supersonic flow induces bow shock
waves in front of the jet trajectories. If the
structures block the stream to a certain extent, the
geometric throat doesnt work and a gasdynamic
throat takes over the responsibility of the subsonicsupersonic transition. Fig. 4 show the flowfield in
the z = 0 symmetric plane of the upper jet orifice. In
Fig. 4(a), an irregular gasdynamic throat can be
seen, as well as a barrel shock and a Mach disc
associated with the under-expanded jet condition.
The interleaved jets result in an asymmetric
flowfield with the sharp gradients of velocity,
temperature and species concentration. Fig. 4(b)

Fig. 4

Flow structures in the near field of the jets: (a) Mach number, (b) two-dimensional flow vector in the symmetric
plane of the upper jet orifice (the flow direction is from left to right), (c) temperature distribution and (d) Mach
number in the z-symmetric plane (length unit: cm).
212

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

Fig. 5

Fig. 6(c) and (d), the rich-iodine condition results in

a higher concentration of I* in the core flow than
that in case (1). Within both the jet trajectories,
there are two strip-like regions which are abundant
in ground state iodine atom as seen in Fig. 6(e) and
(f). Such regions indicate that the chain reactions of
dissociation illustrated by Fig. 1 go on successfully
benefiting from the sufficient molecular iodine
within the jet trajectories and efficient mixing
between the singlet oxygen flow and the jet flow.
Just in the same region, however, I* is infrequent as
shown in Fig. 6(d), which is ultimately different
from the more homogeneous distribution of I* in
Fig. 6(c). The nonlinear effects associated with the
concentration of the molecular iodine fed by jet
flows can hardly be predicted by using any
analytical methods. The three-dimensional CFD is
helpful to gain understanding of the physical and
chemical processes which take place in the COIL
mixing nozzles.

Three-dimensional distribution of I2 density

(normalized by the primary flow density).

The distribution contour of molecular iodine is

given in Fig. 5. In this figure, one can find the
three-dimensional structures of the under-expanded
jet and the interaction among the interleaved jet
trajectories. The blockage induced by the jet
structures should be considered carefully. The
subsonic jet mixing configurations which were used
in the traditional COIL design may worth further
studying for COIL operating without primary buffer
gas even though they were not recommended in
preceding research (Fang, Sang and Chen, 2002 &
2003).
The distribution of several species is shown in
Fig. 6 for case (1) and case (2). The molecular
iodine maintains mostly in the near field of jets as
indicated in Fig. 6(a) and (b), while in the richiodine case the iodine molecules extends more
deeply into the downstream field of the nozzle flow
than that in case (1). As shown in Fig. 6(c) and (f),
the interleaved jets cause the heterogeneous
distributions of laser media, the excited and ground
state iodine. Comparing the distributions of I* in

3.4

Spatial distribution of the small signal gain

The distribution of the small signal gain for the two

cases is plotted in Figs. 7. In Fig. 7(a), the small
signal gain degenerates to zero in the boundary
layers, which letdowns the laser energy density.
The spanwise distribution of the small signal gain in
case (1) appears more reasonable than that of case
(2) as shown in Fig. 7(b). Negative gain strips along
the jet trajectories can be seen in Fig. 7(b). The
distributional disparity between I* and I shown in
Fig. 6(d) and (f) accounts for this negative gain
since the small signal gain is proportional to the
molar concentration difference, I * 0.5[I ] . The
details of the spatial distribution and the evolution
of the small signal gain will be explained in a
follow-up research article in the future.

([ ]

213

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

Fig. 6

Distribution of species in the longitudinal planes: molecular iodine for (a) case 1 and (b) case 2; pumped iodine
atom for (c) case 1 and (d) case 2; ground-state iodine atom for (e) case 1 and (f) case 2 (length unit: cm, density
normalized by the primary flow density).

Fig. 7

Distribution of the small signal gain in y-z planes (x=8cm): (a) case 1I2/ N2=1/60 and (b) case 2I2/ N2=2/59
(gain unit: cm-1, length unit: cm).
214

Engineering Applications of Computational Fluid Mechanics Vol. 1, No. 3 (2007)

6. Furman D, Bruins E, Rybalkin V, Barmashenko

BD and Rosenwaks S (2001). Parameter study
of small-signal gain in s slit nozzle, supersonic
chemical oxygen-iodine laser operating without
primary buffer gas. IEEE Journal of Quantum
Electronics 37(2):174182.
7. Hu ZM, Gao YL, Zhang DL, Yang GW and
Jiang ZL (2004). Numerical simulation of
gaseous detonation reflection on wedges with a
detailed chemical reaction model. Acta
Mechanica Sinica 36(4):385392 (Chinese
series).
8. Jiang ZL, Takayama K and Chen YS (1995).
Dispersion conditions for non-oscillatory shock
capturing schemes and its applications.
Computational Fluid Dynamics Journal
4:137150.
9. Jiang ZL (2004). On dispersion-controlled
principles for non-oscillatory shock-capturing
schemes. Acta Mechanica Sinica 20(1):115.
10. Miller JH, Shang JS, Tomaro RF and Strang
WZ (2000). Computation of COIL nozzle
flowfields with transonic injection. AIAA Paper
20002575. 31st AIAA Plasmadynamics and
Lasers Conference, 19-22 June 2000, Denver,
CO, USA.
11. Paschkewitz J, Shang J, Miller J and Madden T
(2000). An assessment of COIL physical
property and chemical kinetic modelling
methodologies. AIAA Paper 20002574. 31st
AIAA Plasmadynamics and Lasers Conference,
19-22 June 2000, Denver, CO, USA.
12. Rybalkin V, Katz A, Furman D, Barmashenko
BD and Rosenwaks S (2002). Spatial
distribution of the gain and temperature across
the flow in a slitnozzle supersonic chemical
oxygen-iodine laser with transonic and
supersonic schemes of iodine injection.
IEEE Journal of Quantum Electronics
38(10):13981405.
13. Steger JL and Warming RF (1981). Flux vector
splitting of the inviscid gas dynamic equations
with applications to finite difference method.
Journal of Computational Physics 40:263293.

4. CONCLUSIONS

Through the numerical studies of the mixing and

reactive flow in the mixing nozzle of COIL
operating without primary buffer gas, some
important observations are summarized as follows:
A reasonable distribution and spatial evolution of
the small signal gain can be obtained through the
interleaved injection scheme if a proper molecular
iodine flow rate is used in the mixing nozzle. The
laser medium distribution suffers from the iodine
flow rate because of the nonlinear interaction
between the fluid dynamics of jets and the chain
reaction kinetics of the self-catalyzed I2
dissociation.
ACKNOWLEDGEMENTS

The authors would like to thank Prof. B.G. Wu and

Prof. G.W. Yang at the Institute of Mechanics. This
work was partially supported by Korea Research
Foundation Grant No. KRF-2005-005-J09901.
REFERENCES

1. Bruins E, Furman D, Rybalkin V, Barmashenko

BD and Rosenwaks S (2002). One-dimensional
modeling of the gain and temperature in a
supersonic chemical oxygen-iodine laser with
transonic injection of iodine. IEEE Journal of
Quantum Electronics 38(4):345352.
2. Eppard WM, McGrory WD, Godfrey AG, Cliff
EM and Borggaard JT (2000). Recent advances
in numerical techniques for the design and
analysis of COIL system. AIAA Paper
20002576. 31st AIAA Plasmadynamics and
Lasers Conference, 19-22 June 2000, Denver,
CO, USA.
3. Fang BJ, Sang FT and Chen F (2002).
Experimental study on kW COIL with nitrogen
buffer gases. High Power Laser and Particle
Beam 14(3):447450.
4. Fang BJ, Sang FT and Chen F (2003). Design
and experimental results of 2kW COIL with
nitrogen buffer gases. High Power Laser and
Particle Beam 15(12):11481150.
5. Furman D, Barmashenko BD and Rosenwaks S
(1998). Parameter study of an efficient
supersonic chemical oxygen-iodine laser/jet
generate system operating without buffer gas.
IEEE Journal of Quantum Electronics
34:10681074.
215