\titlea{Introduction}
Shock waves propagating in a particle-laden gas medium
have many applications.
\begexample{1.1}
 Rocket plumes from the metallized
 propellant and detonation in air laden with liquid droplets
 or solid particles such as coal and grain dust are two examples.
\endexample
Heating and acceleration of particles by shock waves is important
in these phenomena to understand the evaporation and / or
ignition process of the liquid droplets and small solid particles...

Numerous research papers have appeared on the subject of
gas-particle two phase flows. Carrier [4], Kriebel [13], and
Rudinger [19], studied the relaxation zone behind a normal shock
wave propagating in a gas-particle medium...

This phenomenon, known as
shock-wave diffraction, is influenced by many parameters such
as incident shock Mach number $M_s$, corner wedge angle
$\theta _w$, and
initial gas pressure and temperature, which are depicted in
Fig. 1. The particles are suspended uniformly in a limited region,
which would be an idealized case of the experiment...

In this paper we consider a gas-particle two phase flow and numerically
solve the two-dimensional Mach-reflection from a compression corner
using high resolution TVD scheme by Harten [10]. Figure 2
schematically shows a typical single Mach reflection of a moving shock
wave, consisting of incident shock, reflected shock, Mach stem and
slipline. We employed here the particle parameters satisfying the dilute
gas-particle flow proposed by Growe [9], and subsequently the two
fluid model.
%
\titlea{Governing equations}
The following assumptions are used in the present
investigation for the two-phase flow analysis. A similar
set of assumptions can be found in [3].
%
\medskip{\parindent=0pt
{\it a.} The gas is perfect.

{\it b.} The particles do not undergo a phase change.

{\it c.} The volume occupied by the particles is negligible.

{\it d.} The particles do not interact with each other.

{\it e.} The thermal and Brownian motions of the particles are
negligible.

{\it f.} The particles are solid spheres of uniform diameter and
have a constant material density...
}\medskip
%
The scales used  in  normalizing  the  variables  in   the
governing equations are $\bar {L}, \bar {U_r}, \bar {\rho _0},
\bar {p_0}, \bar {T_0}$ : $\bar L$ is the reference
length equal to the height of the planar shock tube\fonote{A tube is a
longish hollow object.}
and $\bar {U_r}$ is the reference velocity  equal  to
$(\bar {p_0}/\bar {\rho _0})^{1\over 2}$. Nondimensional
quantities can then be expressed as
%
$$ t = \bar t \bar {U_r}/\bar L, \ \ x = \bar x/\bar L
, \ \ y = \bar y/\bar L, $$
$$\rho = \bar \rho/\bar \rho_{0}, \ \ p = \bar p/\bar p_{0},
 \ \ T = \bar T/\bar T_{0},$$
$$u = \bar u/\bar U_{r}, \ \ etc. $$
%
Here $\rho$ is  the  gas density of the mixture, $u$ and $v$ the  gas
velocity components,  $p$  the  gas  pressure,  and $T$ the   gas
temperature...

\titleb{Gas phase equations}
The two-dimensional  unsteady Euler equations for the gas  phase
are in the Cartesian coordinate system are in vector form
%
$${\partial \vec Q\over \partial t}+{\partial \vec F\over \partial x}
+{\partial \vec G\over \partial y}+ \vec I = 0 \eqno(1)$$
%
\begfig 5.2cm
\figure{1}{Geometry of computational domain}
\endfig
%
\begfig 5.4cm
\figure{2}{Schematics of the reflected shock waves}
\endfig
%
Here $Q$ is the vector of conservation variables, and $F$ and $G$ are
flux vectors, and $I$ the source terms for the momentum and energy
balances given by
%
$$
\vec Q=
\left[\matrix{
\rho \cr
\rho u \cr
\rho v \cr
e \cr
}\right]
\ \ ,\ \ \ \ \ \
\vec F=
\left[\matrix{
\rho u \cr
\rho u^2+p \cr
\rho uv \cr
e(u+p) \cr
}\right]
$$
$$\eqno(2)$$
$$
\vec G=
\left[\matrix{
\rho v \cr
\rho uv \cr
\rho v^2+p \cr
e(v+p) \cr
}\right]
\ \ ,\ \ \ \ \ \
\vec I=
\left[\matrix{
0 \cr
f_{px} \cr
f_{py} \cr
q_{p}+u_{p}f_{px}+v_{p}f_{py} \cr
}\right]
$$
%
where
%
$$ e = {p\over (\gamma -1)}+{1\over 2}\rho (u^{2}+v^{2}) \eqno(3) $$
%
Here $\gamma$ is the gas specific heat ratio and $e$ is the total energy
of the gas.
%
\titleb{Particle phase equations}
The drag coefficient $C_d$ is usually expressed as a function of the
particle Reynolds number $\Rey_p$, for example, as done by Clift et
al. [6].
For $\Rey_p < 800$
%
$$ C_{d} = {24\over \Rey_{p}}(1+0.15 \Rey_{p}^{0.687})\eqno(11) $$
%
and for $800 < \Rey_p < 3\times 10^5$
%
$$ C_{d} ={24\over \Rey_{p}}(1+0.15 \Rey_{p}^{0.687})
+{0.42\over (1+42500 \Rey_{p}^{-1.16})} $$
%
Here the particle Reynolds number is defined by
%
$$ \Rey_{p}=(\rho \Delta U {\bar D}_{p}/{\bar \mu}){\bar U}_{r}\bar
{\rho _0}
\eqno(12) $$
%
Temperature dependence of the dynamic viscosity of the gas phase $\bar
\mu$, used here is similar to the one adopted by Martsiano et al. [16]
%
$$ {\bar \mu} = {\bar \mu}_{r}({\bar T}/{\bar Tr})^{0.65} \eqno(13)$$
%
where $\bar \mu _r$ is the dynamic viscosity of the gaseous phase at the
reference temperature $\bar T_r$.

The  heat exchange between the two phases is, on  the  other
hand,
%
$$ q_{p} = (6 Nu/Pr)({\bar \mu}\rho_{p}/{\bar D}_p{\bar \sigma}_p)
(\kappa {\bar L}/{\bar U}_r)(T-T_{p}) \eqno(14) $$
%
\begtheorem{2.1}
The Nusselt number can be expressed in terms  of  Reynolds
number and Prandtl number.
%
$$ Nu = 2 + 0.459 \Rey_{p}^{0.55}Pr^{0.33} \eqno(15)$$
%
where the Reynolds number is defined as above.
\endtheorem
This result was reported by Drake [8].
%
\titlea{Numerical approach}
\titleb{Numerical schemes}
\titlec{Gas phase}
The gas phase equations  (1), are a hyperbolic system. We
have used a high-resolution shock-capturing method, a class  of
total variation diminishing (TVD) scheme [10], and
Roe's average [18]. This explicit scheme is second-order
accurate both in time and in space coordinates.

It has the form in the $\xi$ direction:
%
$$ Q_{j,k}^{n+1} = Q_{j,k}^{n}-\lambda _{\xi}[\hat F_{j+{1 \over 2},k}^n
-\hat F_{j-{1 \over 2},k}^n] \eqno(16) $$
%
where $\lambda _{\xi}=\Delta t/\Delta \xi$, and $Q_{j,k}^{n+1}$ is the
numerical solution at space $\xi = j\Delta \xi$, $\eta = k\Delta \eta$
and time $t=(n+1)\Delta t$...

Also,
%
$$\eqalignno{&\gamma_{j+{1 \over 2},k}^l = \hfill\cr
             &\quad\cases{
                  \sigma({c_\xi} _{j+{1 \over 2},k}^l)
                   (g_{j+1}^l-g_{j}^l)/\alpha_{j+{1 \over 2},k}^l
                   & if $\alpha_{j+{1 \over 2},k}^l \ne 0$ \cr
                  0& if $\alpha_{j+{1 \over 2},k}^l = 0$ \cr}
                                                              &(20)}$$
%
$$ \psi(z) = \cases{|z|&if\ \ $|z| \geq \delta_{1}$ \cr\cr
(z^{2}+ \delta_{1}^2)/2\delta_{1}&if\ \ $|z| < \delta_{1}$ \cr}
\eqno(21) $$
%
The function $g_{j,k}^l$ is the limiter function and $\psi (z)$
is an entropy correction to $|z|$ where $\delta _1$ is a small
positive parameter. More
details can be found in [22]. The numerical approach in the $\eta$
direction is similar to the above procedures.
%
\titlec{Particle phase}
The particle properties do not change discontinuously across a shock
contrary to the gas phase. The MacCormack Scheme [14] is
therefore suitable to sovle the particle phase equations. It also has
second-order accuracy both in time and space and has already been
successfully applied to steady two-phase transonic nozzle flow problems
by Chang [5]. The particular scheme has the following predictor and
corrector expressions:
%
$$ \eqalign{Q_{p \ k,j}^* &=Q_{p \  k,j}^n \cr
&-\lambda_{\xi}\{(\xi_{x})_{k,j}(F_{p \ k,j}^n - F_{p \ k-1,j}^n)\cr
&\quad\quad
+(\xi_{y})_{k,j}(G_{p \ k,j}^n - G_{p \ k-1,j}^n)\} \cr
&-\lambda_{\eta}\{(\eta_{x})_{k,j}(F_{p \ k,j}^n - F_{p \ k,j-1}^n)\cr
&\quad\quad
+(\eta_{y})_{k,j}(G_{p \ k,j}^n - G_{p \ k,j-1}^n)\} \cr
&-\Delta t I_{p \ k,j}^n \cr} \eqno(22) $$
%
$$ \eqalign{Q_{p \ k,j}^{n+1} &={1\over 2}[Q_{p \  k,j}^n
+Q_{p \ k,j}^* \cr
&-\lambda_{\xi}\{(\xi_{x})_{k,j}(F_{p \ k+1,j}^* - F_{p \ k,j}^*)\cr
&\quad\quad
+(\xi_{y})_{k,j}(G_{p \ k+1,j}^* - G_{p \ k,j}^*)\} \cr
&-\lambda_{\eta}\{(\eta_{x})_{k,j}(F_{p \ k,j+1}^* - F_{p \ k,j}^*)\cr
&\quad\quad
+(\eta_{y})_{k,j}(G_{p \ k,j+1}^* - G_{p \ k,j}^*)\} \cr  &-\Delta t
I_{p \ k,j}^*] \cr} \eqno(23) $$
%
\begfig 11.3cm
\figure{3a,b}{Density distribution on the ramp for gas-only flow.
{\bf a} present result, {\bf b} experimental result.
Dashed line; Ben-Dor and Glass [2], Dotted line; Deschambault and
Glass [7]}
\endfig
%
\begfig 13.1cm
\figure{4a,b}{Mach number contours. {\bf a} Gas-only flow,
{\bf b} $\bar D_p$=1 $\mu$m ($\varphi$=0.23)}
\endfig
%
\titlec{Boundary conditions}
The flow data on the boundary CC',
once swept by the shock wave, are obtained from the moving shock
Rankine-Hugoniot relations:
%
$$\eqalign{
p_{1}/p_{0} &= [2 \gamma M_{s}-(\gamma -1)]/(\gamma +1)\cr
\rho _{1}/\rho _{0} &= [\Gamma (p_{1}/p_{0})+1]/
[\Gamma +(p_{1}/p_{0})]\cr
u_{1} &= M_{s}\{1-[(\gamma -1)M_{s}^2 +2]/[(\gamma
+1)M_{s}^2]\}c_{0}\cr} \eqno(24) $$
%
where $\Gamma =(\gamma +1)/(\gamma-1)$ and
$c_0$ is the speed of sound  in  the stagnation state...
%
\titleb{Coupling of the gas and particle equations}
Coupling of the two phase equations is possible using the method of
fractional steps. The gas phase equation in the conservation from is
split into two equations:
%
$$\eqalign{{\partial_{t}Q}+{\partial_{\xi}F} &= 0 \cr
%
{\partial_{t}Q}+{\partial_{\eta}G}+I &= 0 \cr} \eqno(25) $$
%

...At new time steps, (26) is solved using  the  particle  phase
properties of the previous calculation. (27) is then solved using
the calculated gas phase information.
%
\titlea{Results and discussion}
\titleb{Gas-only phase flow}
To test the validity of the present computational procedure,
a gas-only flow was first calculated for the  same  compression
corner, with the perfect-gas assumption and the  specific  heat
ratio $\gamma$=1.4. For the flow data $M_s$=2.03 and $\theta _w$ =
27$^\circ$, the shock wave...

The  shock  wave  configuration and  the  isopycnic  curves
plotted from the computed results, Fig. 3a, agreed well with the
experimental interferogram presented by Deschambault et
al. [7], shown in Fig. 3b. The slipline which is a
weak discontinuity which differentiates the thermodynamic regions
of different density, can be clearly seen...
%
\titleb{Particle-laden gas flow}
The  data  base of particles, uniformly distributed  in  the
region bounded by the two contact discontinuities, AA' and BB', is...

\begtabfull \tabcap{1}{Here we transfered a small portion of the article
into a table to demonstrate the use of the appropriate macros}
\halign {#\hfill\quad&#\hfill\quad&#\hfill\cr
\noalign{\hrule\smallskip}
&Symbol&value\cr
\noalign{\smallskip\hrule\smallskip}
Particle diameter,   &$\bar D_p$       &
   1 $\mu m$, 3 $\mu m$, 5 $\mu m$, and 10 $\mu m$.\cr
Mass fraction ratio, &$\varphi$        & 0.33, 0.23, 0.09.\cr
Mass density,        &$\bar \sigma _p$ & 4000 kg/m$^3$.\cr
Specific heat,       &$\bar C_s$       & 1380 J/kg-K.\cr
                     &$Pr$             &0.75,\cr
                     &$M_s$            &2.03\cr
                     &$\theta _w$      &27$^\circ$.\cr
\noalign{\smallskip\hrule}
}
\endtab

The overall configuration
appears to be somewhat similar to the  double  Mach  reflection
discussed by Ben-Dor and Glass [2], but the present case has the
peculiar cross-section of  the contact  discontinuity  and  the
reflected shock wave, a kind of weak triple point, without any
second Mach stem emanated from it....
%

It is observed that the particles are
accelerated more rapidly for the smaller particles, and that the
flow is therefore  divided  by  the  particle  free  zone,  the
relaxation zone, and the inactive particle zone...
%
\titlea{Conclusions}
The numerical  analyses  of  the  progressive   shock-wave
reflection in the  particle-laden  gas  medium,  performed  for
various particle sizes  and mass fraction ratios,  led  to  the
following conclusions:
\medskip
{\parindent=0pt
{\it 1.} The numerical results agreed well with the experimental results
for the case of gas-only phase flow. It exhibited supersonic flow in the
delta region enclosed by the Mach stem and the slipline.

{\it 2.} For the smaller particle sizes and higher mass fraction ratios,
the relaxation zone can include an extended supersonic flow region
generated by the reflected nearly straight local shock.

{\it 3.} The contact discontinuity becomes very sharp in the relaxation
zone for the reduced particle sizes and increased mass fraction ratios.

{\it 4.} For the smaller particle sizes and higher mass fraction ratios,
the shock wave is propagated with lower speed and the peak density in
the relaxation zone grows very rapidly with time.
} %end of list

\begref{References}{19.}
%
\noindent For examples of the references please refer to the journal.
The coding is discussed in detail in the documentation.
\endref
\bye