A Discrete-Ordinates Solution for the Strong Evaporation Problem in Rarefied Gas Dynamics

ABSTRACT

In this work we solve the nonlinear strong evaporation problem in rarefied gas dynamics. The analysis is based on the BGK model, with three-dimensional velocity, derived from the Boltzmann equation. We present the complete development of a closed form solution for evaluating density, velocity and temperature perturbations and the heat flux of a gas. Numerical results are presented and discussed.

Keywords:
rarefied gas dynamics; kinetic model; ADO method; strong evaporation

1 INTRODUCTION

In the previous work ¹³, the Analytical Discrete Ordinates (ADO) method ⁴ was used to solve the strong evaporation problem based on the BGK model ⁵, considering a monoatomic gas with one-dimensional velocity. This problem was formulated by Arthur and Cercignani ² and solved by Siewert and Thomas Jr. ¹⁸ where was used the Elementary Solutions method ⁶. With the purpose of improve the physical model, Siewert and Thomas Jr. ¹⁹ used the BGK equation, considering a gas with three-dimensional velocity, and the Elementary Solutions method for obtaining numerical results for the problem. In this work, the ADO method is used to solve the strong evaporation problem with the BGK equation, considering also a monoatomic gas, but with three-dimensional velocity.

In regard to the ADO solution, this method is an analytical version of the Discrete Ordinates method proposed by Chandrasekhar ⁸ for radiative transfer problems. The ADO method has been used to solve, in a unified manner, a class of flow ¹⁵, heat-transfer ¹⁶ and weak evaporation/condensation ¹⁴ problems in rarefied gas dynamics, where was considered a rectangular and one-dimensional geometry. Still, the ADO solution has been shown to be concise, accurate, easy to implement and efficient from a computational point of view. These good aspects of the solution can be attributed to the fact that the method is analytical in the spatial variable, generates results with accuracy of several digits and has a computational implementation constituted by simple steps, such as solve an eigenvalue problem or a linear system. Furthermore, the computational time of execution of the method is very small in a personal computer.

Research activity on the kinetics of vapors near to interphase boundaries has increased significantly over the last few years ^{1), (3), (10), (11), (14), (17), (21), (22), (24), (25), (26}, reflecting this way in a large number of problems in science and technology, to which evaporation and condensation phenomena are of practical importance. As noted in Ytrehus ²⁵ work, problems of this kind are found in diversified areas, such as upper atmosphere meteorology, the sodium cooling of nuclear reactors, projects of spacecraft experiments, petrochemical engineering, vacuum technology and the interaction of high-power laser radiation with metal surfaces.

Following Yasuda, Takata and Aoki ²⁴, it is widely accepted that gas flows around a body are usually studied by the Navier-Stokes set of equations with the nonslip condition for the flow velocity and the nonjump condition for the temperature. However, such conditions are invalid if a phase transition (evaporation or condensation) takes place at the surface of the body, i.e., if vapor flows around condensed phases are considered. In fact, it is known that there is a difference between the temperature of the surface of the condensed phase and the temperature of the gas in contact with this surface. The same is true for the pressure: the pressure of the vapor at the surface is different from its saturation pressure at the temperature of the surface. These differences are called jumps of temperature and pressure, respectively. In order to study gas flows with a phase transition by the fluid-dynamic set of equations, at physical problems of the previously mentioned areas, one needs boundary conditions that describe these jumps correctly. Such conditions are called jump conditions for evaporation and condensation.

Still following Yasuda, Takata and Aoki ²⁴, in kinetic theory the determination of the jumps for evaporation and condensation is done solving a half-space boundary-value problem of the Boltzmann equation ^{7), (23}. When the evaporation or condensation is strong, either the nonlinear Boltzmann equation or nonlinear models should be used to obtain the jump conditions for the compressible Euler set of equations.

In the strong evaporation problem the distribution function, in the original nonlinear equation, is linearized around a downstream Maxwell distribution containing a drift velocity u _∞. This procedure is described in Ref. ². In this work, firstly, the linearized version of the problem is solved to evaluate density, velocity and temperature perturbations and the heat flux of the gas. Then, the analytical ADO solution is used in the original nonlinear version of the model, in order to get a second set of results for the quantities of interest. Numerical results are presented for both approaches.

2 THE KINETIC MODEL

Following Scherer ¹³, we consider here the steady-state limit of the one-dimensional and time dependent problem proposed by Ytrehus ²⁵: a liquid (or solid) is initially in equilibrium with its pure vapor which occupies the half-space $x\ge 0$ at uniform temperature and pressure T ₀ and p ₀, respectively. At time $t=0$ the pressure level in the vapor changes discontinuously to the value p _∞ and it is kept at this (constant) value. Then, evaporation or condensation begins, through the phase boundary, according to whether the pressure level p _∞ is below or above the saturation pressure p ₀, respectively.

Continuing to follow Ytrehus ²⁵ analysis, it is reasonable to assume that far downstream of the boundary, after a sufficiently long time, a steady-state will be accomplished. The flow far from the phase boundary is then a uniform equilibrium flow with constant parameters n _∞, u _∞ and T _∞. In this way, a kinetic boundary layer will form between the phase boundary and the downstream equilibrium region, in which, nonequilibrium effects may influence significantly the motion of the vapor ²⁵.

In this context, we describe here the state of gas by the nonlinear BGK model ⁵, with three degrees of freedom, which can be written as

$\mathbf{v}·{\nabla }_{r}f(\mathbf{r},\mathbf{v})=\eta \left[\varphi \right(\mathbf{r},\mathbf{v})-f(\mathbf{r},\mathbf{v}\left)\right],$ (2.1)

where $f(\mathbf{r},\mathbf{v})$ is the distribution function, $\mathbf{r}=(x,y,z)$ is the spatial vector, $\mathbf{v}=(v_x,v_y,v_z)$ is the molecular velocity vector, η is an appropriate collision frequency, $\varphi (\mathbf{r},\mathbf{v})$ is a local Maxwell distribution given by

$\varphi (\mathbf{r},\mathbf{v})=n\left(\mathbf{r}\right){\left[\frac{m}{2\pi kT\left(\mathbf{r}\right)}\right]}^{3/2}\text{exp}\left\{-\frac{m|\mathbf{v}-\mathbf{u}(\mathbf{r}){|}^2}{2kT\left(\mathbf{r}\right)}\right\},$ (2.2)

m is the mass of a gas particle and k is the Boltzmann constant ⁷. Following Williams ²³, the density n(r), mass velocity u(r) and temperature T (r) in Eq. (2.2) are given by

$n\left(\mathbf{r}\right)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(\mathbf{r},\mathbf{v})\text{d}\mathbf{v},$ (2.3)

$\mathbf{u}\left(\mathbf{r}\right)=\frac{1}{n\left(\mathbf{r}\right)}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathbf{v}f(\mathbf{r},\mathbf{v})\text{d}\mathbf{v}$ (2.4)

and

$T\left(\mathbf{r}\right)=\frac{m}{3kn\left(\mathbf{r}\right)}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }|\mathbf{v}-\mathbf{u}(\mathbf{r}\left){|}^{2}f\right(\mathbf{r},\mathbf{v})\text{d}\mathbf{v}.$ (2.5)

Still, the heat flux is given by

$\mathbf{q}\left(\mathbf{r}\right)=\frac{m}{2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }|\mathbf{v}-\mathbf{u}(\mathbf{r}\left){|}^2\right[\mathbf{v}-\mathbf{u}\left(\mathbf{r}\right)\left]f\right(\mathbf{r},\mathbf{v})\text{d}\mathbf{v}.$ (2.6)

We assume that far downstream the gas relaxes to an equilibrium distribution characterized by steady drift velocity u _∞, density n _∞ and temperature T _∞,

$f_{\infty }\left(\mathbf{v}\right)=n_{\infty }{\left[\frac{m}{2\pi k{T}_{\infty }}\right]}^{3/2}\text{exp}\left\{-\frac{m|\mathbf{v}-{\mathbf{u}}_{\infty }{|}^2}{2k{T}_{\infty }}\right\},$ (2.7)

where ${\mathbf{u}}_{\infty }=(u_{\infty },0,0)$. At this point, we follow Arthur and Cercignani ² and we linearize $f(\mathbf{r},\mathbf{v})$ and $\varphi (\mathbf{r},\mathbf{v})$ around $f_{\infty }\left(\mathbf{v}\right)$. In this manner, we write $f(\mathbf{r},\mathbf{v})$ as

$f(\mathbf{r},\mathbf{v})=f_{\infty }\left(\mathbf{v}\right)[1+h(\mathbf{r},\mathbf{v}\left)\right],$ (2.8)

where $h(\mathbf{r},\mathbf{v})$ is a perturbation to the absolute Maxwellian $f_{\infty }\left(\mathbf{v}\right)$. Thus, we substitute Eq. (2.8) into Eq. (2.1) and linearize $\varphi (\mathbf{r},\mathbf{v})$ around $f_{\infty }\left(\mathbf{v}\right)$ to obtain, for the x direction (normal in relation to the interface), the linearized equation written in terms of the perturbation function h,

$c_x\frac{\partial }{\partial \tau }h(\tau ,\mathbf{c})+h(\tau ,\mathbf{c})={\pi }^{-3/2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-|\mathbf{c}\mathbf{\text{'}}-\mathbf{u}{|}^2}F(\mathbf{c}\mathbf{\text{'}},\mathbf{c}:\mathbf{u})h(\tau ,\mathbf{c}\mathbf{\text{'}})\text{d}\mathbf{c}\mathbf{\text{'}},$ (2.9)

where

$F(\mathbf{c}\mathbf{\text{'}},\mathbf{c}:\mathbf{u})=1+2(\mathbf{c}\mathbf{\text{'}}-\mathbf{u})·(\mathbf{c}-\mathbf{u})+(2/3)\left(\right|\mathbf{c}\mathbf{\text{'}}-\mathbf{u}{|}^2-3/2\left)\right(|\mathbf{c}-\mathbf{u}{|}^2-3/2),$ (2.10)

is the scattering kernel and

$\tau =x\eta [m/(2k{T}_{\infty }){]}^{1/2},$ (2.11 a)

$\mathbf{c}=\mathbf{v}[m/(2k{T}_{\infty }){]}^{1/2}$ (2.11 b)

and

$\mathbf{u}={\mathbf{u}}_{\infty }[m/(2k{T}_{\infty }){]}^{1/2},$ (2.11 c)

are dimensionless variables. We note that $\mathbf{u}=(u,0,0)$ is the normalized downstream drift velocity.

2.1 Boundary Conditions

To obtain the boundary condition for the interface, in terms of the perturbation h, we follow Siewert and Thomas Jr. ¹⁸ and set $x=0$ in Eq. (2.8), to find (for $v_x>0$)

$h(0,\mathbf{v})=\frac{f_0\left(\mathbf{v}\right)-f_{\infty }\left(\mathbf{v}\right)}{f_{\infty }\left(\mathbf{v}\right)},$ (2.12)

where $f_0\left(\mathbf{v}\right)$ is the Maxwellian distribution, given by Eq. (2.2), evaluated at $x=0$

$f_0\left(\mathbf{v}\right)=n_0{\left[\frac{m}{2\pi k{T}_0}\right]}^{3/2}\text{exp}\left\{-\frac{m|\mathbf{v}-{\mathbf{u}}_0{|}^2}{2k{T}_0}\right\},$ (2.13)

where ${\mathbf{u}}_0=(u_0,0,0)$ is the mass velocity at the interface. We then linearize $f_0\left(\mathbf{v}\right)$ around $f_{\infty }\left(\mathbf{v}\right)$ to obtain the dimensionless boundary condition (for $c_x>0$)

$h(0,\mathbf{c})=\Delta N_0+2(c_x-u)(u_w-u)+\left(\right|\mathbf{c}-\mathbf{u}{|}^2-3/2)\Delta T_0,$ (2.14)

with dimensionless variables defined by

$u_w=u_0[m/(2k{T}_{\infty }){]}^{1/2},$ (2.15 a)

$\Delta N_0=\frac{n_0-n_{\infty }}{n_{\infty }}$ (2.15 b)

and

$\Delta T_0=\frac{T_0-T_{\infty }}{T_{\infty }}.$ (2.15 c)

On the other hand, when $x\to \infty ,f(\mathbf{r},\mathbf{v})$ approaches $f_{\infty }\left(\mathbf{v}\right)$ and, looking back at Eq. (2.8), we then find the condition

$\underset{\tau \to \infty }{\lim }h(\tau ,\mathbf{c})=0.$ (2.16)

2.2 Physical Quantities of Interest

We substitute Eq. (2.8) into Eqs. (2.3) to (2.6) to find, in terms of h, the density perturbation

$\Delta N\left(\tau \right)={\pi }^{-3/2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-|\mathbf{c}-\mathbf{u}{|}^2}h(\tau ,\mathbf{c})\text{d}\mathbf{c},$ (2.17)

the velocity perturbation

$\Delta U\left(\tau \right)=\frac{{\pi }^{-3/2}}{u}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-|\mathbf{c}-\mathbf{u}{|}^2}(c_x-u)h(\tau ,\mathbf{c})\text{d}\mathbf{c},$ (2.18)

the temperature perturbation

$\Delta T\left(\tau \right)=\frac{2}{3}{\pi }^{-3/2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-|\mathbf{c}-\mathbf{u}{|}^2}\left(\right|\mathbf{c}-\mathbf{u}{|}^2-3/2\left)h\right(\tau ,\mathbf{c})\text{d}\mathbf{c}$ (2.19)

and the heat flux

$Q_x\left(\tau \right)={\pi }^{-3/2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-|\mathbf{c}-\mathbf{u}{|}^2}\left(\right|\mathbf{c}-\mathbf{u}{|}^2-5/2\left)\right(c_x-u\left)h\right(\tau ,\mathbf{c})\text{d}\mathbf{c}.$ (2.20)

3 A REFORMULATION IN TERMS OF MOMENT EQUATIONS

Since the quantities we want to evaluate, Eqs. (2.17) to (2.20), are defined in terms of moments of the h function, we look for simpler problems written in terms of those moments. In this way, we first define

$h_1(\tau ,c_x)={\pi }^{-1/2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-(c_y^2+c_z^2)}h(\tau ,c_x,c_y,c_z)\text{d}{c}_y\text{d}{c}_z$ (3.1)

and

$h_2(\tau ,c_x)={\pi }^{-1/2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\text{e}}^{-(c_y^2+c_z^2)}(c_y^2+c_z^2-1)h(\tau ,c_x,c_y,c_z)\text{d}{c}_y\text{d}{c}_z.$ (3.2)

Then, we multiply Eq. (2.9) by

${\varphi }_1(c_y,c_z)={\text{e}}^{-(c_y^2+c_z^2)}$ (3.3)

and we integrate over all c _y and c _z , such that, (introducing the new notation $c_x=\xi$) we obtain a new equation for h ₁

$\xi \frac{\partial }{\partial \tau }{h}_1(\tau ,\xi )+h_1(\tau ,\xi )={\pi }^{-1/2}{\int }_{-\infty }^{\infty }{\text{e}}^{-(\xi \text{'}-u{)}^2}\times \left[k_{11}(\xi \text{'},\xi :u)h_1(\tau ,\xi \text{'})+k_{12}(\xi \text{'},\xi :u)h_2(\tau ,\xi \text{'})\right]\text{d}\xi \text{'},$ (3.4)

where

$k_{11}(\xi \text{'},\xi :u)=1+2(\xi \text{'}-u)(\xi -u)+(2/3)\left[(\xi \text{'}-u{)}^2-1/2\right]\left[(\xi -u{)}^2-1/2\right]$ (3.5)

and

$k_{12}(\xi \text{'},\xi :u)=(2/3)\left[(\xi -u{)}^2-1/2\right].$ (3.6)

Following analogous procedure, we multiply Eq. (2.9) now by

${\varphi }_2(c_y,c_z)=(c_y^2+c_z^2-1){\text{e}}^{-(c_y^2+c_z^2)}$ (3.7)

and we integrate over all c _y and c _z to obtain (for $c_x=\xi$)

$\xi \frac{\partial }{\partial \tau }{h}_2(\tau ,\xi )+h_2(\tau ,\xi )={\pi }^{-1/2}{\int }_{-\infty }^{\infty }{\text{e}}^{-(\xi \text{'}-u{)}^2}\times \left[k_{21}(\xi \text{'},\xi :u)h_1(\tau ,\xi \text{'})+k_{22}(\xi \text{'},\xi :u)h_2(\tau ,\xi \text{'})\right]\text{d}\xi \text{'},$ (3.8)

where

$k_{21}(\xi \text{'},\xi :u)=(2/3)\left[(\xi \text{'}-u{)}^2-1/2\right]$ (3.9)

and

$k_{22}(\xi \text{'},\xi :u)=2/3.$ (3.10)

Equations (3.4) and (3.8) can be rewritten in vector form as

$\xi \frac{\partial }{\partial \tau }\mathbf{H}(\tau ,\xi )+\mathbf{H}(\tau ,\xi )={\int }_{-\infty }^{\infty }\psi (\xi \text{'}:u)\mathbf{K}(\xi \text{'},\xi :u)\mathbf{H}(\tau ,\xi \text{'})\text{d}\xi \text{'},$ (3.11)

where $\mathbf{H}(\tau ,\xi )$ has components $h_1(\tau ,\xi )$ and $h_2(\tau ,\xi )$,

$\psi (\xi :u)={\pi }^{-1/2}{\text{e}}^{-(\xi -u{)}^2}$ (3.12)

and the components $k_{i,j}(\xi \text{'},\xi :u)$ of the $2\times 2$ matrix $\mathbf{K}(\xi \text{'},\xi :u)$ are defined in Eqs. (3.5), (3.6), (3.9) and (3.10).

Analogously, the boundary conditions given by Eqs. (2.14) and (2.16) are rewritten, for $\xi >0$, as

$\mathbf{H}(0,\xi )=\Delta N_0{\pi }^{1/2}\left[\begin{array}{c}1\\ 0\end{array}\right]+2{\pi }^{1/2}(\xi -u)(u_w-u)\left[\begin{array}{c}1\\ 0\end{array}\right]+\Delta T_0{\pi }^{1/2}\left[\begin{array}{c}(\xi -u{)}^2-1/2\\ 1\end{array}\right]$ (3.13)

and

$\underset{\tau \to \infty }{\lim }\mathbf{H}(\tau ,\xi )=\mathbf{0}.$ (3.14)

In terms of the definitions given in Eqs. (3.1) and (3.2) the quantities of interest are also expressed in a vector form. In this way, the solution for the H problem is used to evaluate the density perturbation

$\Delta N\left(\tau \right)={\pi }^{-1}{\int }_{-\infty }^{\infty }{\text{e}}^{-(\xi -u{)}^2}{\left[\begin{array}{c}1\\ 0\end{array}\right]}^T\mathbf{H}(\tau ,\xi )\text{d}\xi ,$ (3.15)

the velocity perturbation

$\Delta U\left(\tau \right)=\frac{{\pi }^{-1}}{u}{\int }_{-\infty }^{\infty }{\text{e}}^{-(\xi -u{)}^2}{\left[\begin{array}{c}\xi -u\\ 0\end{array}\right]}^T\mathbf{H}(\tau ,\xi )\text{d}\xi ,$ (3.16)

the temperature perturbation

$\Delta T\left(\tau \right)=\frac{2}{3}{\pi }^{-1}{\int }_{-\infty }^{\infty }{\text{e}}^{-(\xi -u{)}^2}{\left[\begin{array}{c}(\xi -u{)}^2-1/2\\ 1\end{array}\right]}^T\mathbf{H}(\tau ,\xi )\text{d}\xi$ (3.17)

and the heat flux

$Q_x\left(\tau \right)={\pi }^{-1}{\int }_{-\infty }^{\infty }{\text{e}}^{-(\xi -u{)}^2}(\xi -u){\left[\begin{array}{c}(\xi -u{)}^2-3/2\\ 1\end{array}\right]}^T\mathbf{H}(\tau ,\xi )\text{d}\xi ,$ (3.18)

where T denote the transpose operation.

To develop an analytical solution to the problem defined by Eq. (3.11), it is convenient to introduce a new function

$\mathbf{G}(\tau ,\xi )={\text{e}}^{-(\xi -u{)}^2}\mathbf{H}(\tau ,\xi )$ (3.19)

such that, we rewrite Eq. (3.11) in the form

$\xi \frac{\partial }{\partial \tau }\mathbf{G}(\tau ,\xi )+\mathbf{G}(\tau ,\xi )=\psi (\xi :u){\int }_{-\infty }^{\infty }\mathbf{K}(\xi \text{'},\xi :u)\mathbf{G}(\tau ,\xi \text{'})\text{d}\xi \text{'},$ (3.20)

where $\mathbf{G}(\tau ,\xi )$ has components $g_1(\tau ,\xi )$ and $g_2(\tau ,\xi )$, with boundary conditions given, for $\xi >0$, by

$\mathbf{G}(0,\xi )={\text{e}}^{-(\xi -u{)}^2}\mathbf{H}(0,\xi )$ (3.21)

and

$\underset{\tau \to \infty }{\lim }\mathbf{G}(\tau ,\xi )=\mathbf{0}.$ (3.22)

In the same way, based on the definition given in Eq. (3.19), we express the density perturbation as

$\Delta N\left(\tau \right)={\pi }^{-1}{\int }_{-\infty }^{\infty }{\left[\begin{array}{c}1\\ 0\end{array}\right]}^T\mathbf{G}(\tau ,\xi )\text{d}\xi ,$ (3.23)

the velocity perturbation

$\Delta U\left(\tau \right)=\frac{{\pi }^{-1}}{u}{\int }_{-\infty }^{\infty }{\left[\begin{array}{c}\xi -u\\ 0\end{array}\right]}^T\mathbf{G}(\tau ,\xi )\text{d}\xi ,$ (3.24)

the temperature perturbation

$\Delta T\left(\tau \right)=\frac{2}{3}{\pi }^{-1}{\int }_{-\infty }^{\infty }{\left[\begin{array}{c}(\xi -u{)}^2-1/2\\ 1\end{array}\right]}^T\mathbf{G}(\tau ,\xi )\text{d}\xi$ (3.25)

and the heat flux

$Q_x\left(\tau \right)={\pi }^{-1}{\int }_{-\infty }^{\infty }(\xi -u){\left[\begin{array}{c}(\xi -u{)}^2-3/2\\ 1\end{array}\right]}^T\mathbf{G}(\tau ,\xi )\text{d}\xi .$ (3.26)

We develop the solution for the G problem: Eq. (3.20) supplemented by Eqs. (3.21) and (3.22), in the next section.

4 A DISCRETE ORDINATES SOLUTION

Equation (3.20) clearly has separable exponential solutions, so we use ν as a separation constant and we seek solutions of the form

$\mathbf{G}(\tau ,\xi )=\mathbf{\Phi }(\nu ,\xi ){\text{e}}^{-\tau /\nu },$ (4.1)

where $\mathbf{\Phi }(\nu ,\xi )$ has components ${\Phi }_1(\nu ,\xi )$ and ${\Phi }_2(\nu ,\xi )$. Thus, we substitute Eq. (4.1) into Eq. (3.20) to obtain

$(1-\xi /\nu )\mathbf{\Phi }(\nu ,\xi )=\psi (\xi :u){\int }_{-\infty }^{\infty }\mathbf{K}(\xi \text{'},\xi :u)\mathbf{\Phi }(\nu ,\xi \text{'})\text{d}\xi \text{'}.$ (4.2)

Following Scherer ¹³ and Siewert and Thomas Jr. ¹⁸, we can use some moments of the function $\mathbf{\Phi }(\nu ,\xi )$ to simplify Eq. (4.2). In this way, firstly, we integrate Eq. (4.2) over all ξ to find

${\int }_{-\infty }^{\infty }{\left[\begin{array}{c}\xi \\ 0\end{array}\right]}^T\mathbf{\Phi }(\nu ,\xi )\text{d}\xi =0.$ (4.3)

Continuing, we can multiply Eq. (4.2) by $(\xi -u)$ and integrate the resulting equation over all ξ to find

${\int }_{-\infty }^{\infty }{\left[\begin{array}{c}{\xi }^2\\ 0\end{array}\right]}^T\mathbf{\Phi }(\nu ,\xi )\text{d}\xi =0.$ (4.4)

Thus, looking to Eqs. (4.3) and (4.4), we rewrite Eq. (4.2), as

$(1-\xi /\nu )\mathbf{\Phi }(\nu ,\xi )=\psi (\xi :u)\mathbf{Q}(\xi :u){\int }_{-\infty }^{\infty }\mathbf{\Phi }(\nu ,\xi \text{'})\text{d}\xi \text{'},$ (4.5)

where the $2\times 2$ matrix

$\mathbf{Q}(\xi :u)=\left[\begin{array}{cc}{q}_{11}(\xi :u)& q_{12}(\xi :u)\\ q_{21}(\xi :u)& q_{22}(\xi :u)\end{array}\right]$ (4.6)

have the components

$q_{11}(\xi :u)=1-2u\left(\xi -u\right)+(2/3)\left(u^2-1/2\right)\left[(\xi -u{)}^2-1/2\right],$ (4.7 a)

$q_{12}(\xi :u)=(2/3)\left[(\xi -u{)}^2-1/2\right],$ (4.7 b)

$q_{21}(\xi :u)=(2/3)\left(u^2-1/2\right)$ (4.7 c)

and

$q_{22}(\xi :u)=2/3.$ (4.7 d)

Still, at this point we note that the exponential term, in Eq. (4.5), can be expressed as

$\psi (\xi :u)={\pi }^{-1/2}\text{cosh}\left(2\xi u\right){\text{e}}^{-({\xi }^2+u^2)}+{\pi }^{-1/2}\text{sinh}\left(2\xi u\right){\text{e}}^{-({\xi }^2+u^2)}$ (4.8)

and, in this manner, we write the final convenient form of Eq. (4.5)

$(1-\xi /\nu )\mathbf{\Phi }(\nu ,\xi )=P(\xi :u)\left[\mathbf{A}\right(\xi :u)+\mathbf{B}(\xi :u\left)\right]{\int }_{-\infty }^{\infty }\mathbf{\Phi }(\nu ,\xi \text{'})\text{d}\xi \text{'},$ (4.9)

where

$P(\xi :u)={\pi }^{-1/2}{\text{e}}^{-({\xi }^2+u^2)}$ (4.10)

and the $2\times 2$ matrixes

$\mathbf{A}(\xi :u)=\left[\begin{array}{cc}{a}_{11}(\xi :u)& a_{12}(\xi :u)\\ a_{21}(\xi :u)& a_{22}(\xi :u)\end{array}\right]$ (4.11)

and

$\mathbf{B}(\xi :u)=\left[\begin{array}{cc}{b}_{11}(\xi :u)& b_{12}(\xi :u)\\ b_{21}(\xi :u)& b_{22}(\xi :u)\end{array}\right]$ (4.12)

have the components

$a_{11}(\xi :u)=[1+2u^2+(2/3\left)\left(u^2-1/2\right)\left({\xi }^2+u^2-1/2\right)\right]\text{cosh}\left(2\xi u\right)-\left[\right(4/3)u^3\xi +(4/3\left)u\xi \right]\text{sinh}\left(2\xi u\right),$ (4.13 a)

$a_{12}(\xi :u)=(2/3)\left({\xi }^2+u^2-1/2\right)\text{cosh}\left(2\xi u\right)-(4/3)u\xi \text{sinh}\left(2\xi u\right),$ (4.13 b)

$a_{21}(\xi :u)=(2/3)\left(u^2-1/2\right)\text{cosh}\left(2\xi u\right),$ (4.13 c)

$a_{22}(\xi :u)=(2/3)\text{cosh}\left(2\xi u\right),$ (4.13 d)

$b_{11}(\xi :u)=[1+2u^2+(2/3\left)\left(u^2-1/2\right)\left({\xi }^2+u^2-1/2\right)\right]\text{sinh}\left(2\xi u\right)-\left[\right(4/3)u^3\xi +(4/3\left)u\xi \right]\text{cosh}\left(2\xi u\right),$ (4.14 a)

$b_{12}(\xi :u)=(2/3)\left({\xi }^2+u^2-1/2\right)\text{sinh}\left(2\xi u\right)-(4/3)u\xi \text{cosh}\left(2\xi u\right),$ (4.14 b)

$b_{21}(\xi :u)=(2/3)\left(u^2-1/2\right)\text{sinh}\left(2\xi u\right)$ (4.14 c)

and

$b_{22}(\xi :u)=(2/3)\text{sinh}\left(2\xi u\right).$ (4.14 d)

Now we rewrite the integral term in Eq. (4.9) as

$(1-\xi /\nu )\mathbf{\Phi }(\nu ,\xi )=P(\xi :u)\left[\mathbf{A}\right(\xi :u)+\mathbf{B}(\xi :u\left)\right]{\int }_0^{\infty }\left[\mathbf{\Phi }\right(\nu ,\xi \text{'})+\mathbf{\Phi }(\nu ,-\xi \text{'}\left)\right]\text{d}\xi \text{'}.$ (4.15)

Then we introduce a (half-range) quadrature scheme $\left[0,\infty \right)$, to approximate the integral term of the above equation, such that

$(1-\xi /\nu )\mathbf{\Phi }(\nu ,\xi )=P(\xi :u)\left[\mathbf{A}\right(\xi :u)+\mathbf{B}(\xi :u\left)\right]\sum _{k=1}^N{w}_k\left[\mathbf{\Phi }(\nu ,{\xi }_k)+\mathbf{\Phi }(\nu ,-{\xi }_k)\right].$ (4.16)

Here ξ _k and w _k are, respectively, N nodes and weights of the (arbitrary) quadrature scheme. If we now evaluate Eq. (4.16) in $\xi =\pm {\xi }_i$, for $i=1,...,N$, and note that $P(\xi :u)$ and $\mathbf{A}(\xi :u)$ are even functions,

$P(\xi :u)=P(-\xi :u),$ (4.17)

$\mathbf{A}(\xi :u)=\mathbf{A}(-\xi :u)$ (4.18)

and $\mathbf{B}(\xi :u)$ is an odd function,

$\mathbf{B}(\xi :u)=-\mathbf{B}(-\xi :u),$ (4.19)

we obtain the discrete-ordinates version of Eq. (4.16) as

$(1\mp {\xi }_i/\nu )\mathbf{\Phi }(\nu ,\pm {\xi }_i)=P({\xi }_i:u)\left[\mathbf{A}\right({\xi }_i:u)\pm \mathbf{B}({\xi }_i:u\left)\right]\sum _{k=1}^N{w}_k\left[\mathbf{\Phi }(\nu ,{\xi }_k)+\mathbf{\Phi }(\nu ,-{\xi }_k)\right].$ (4.20)

We express now Eq. (4.20) in a matrix form, as

$(\mathbf{I}-\mathbf{M}/\nu ){\mathit{\Phi }}_{+}\left(\nu \right)=\mathbf{P}[\mathbf{A}+\mathbf{B}]\mathbf{W}\left[{\mathit{\Phi }}_{+}\right(\nu )+{\mathit{\Phi }}_{-}(\nu \left)\right]$ (4.21)

and

$(\mathbf{I}+\mathbf{M}/\nu ){\mathit{\Phi }}_{-}\left(\nu \right)=\mathbf{P}[\mathbf{A}-\mathbf{B}]\mathbf{W}\left[{\mathit{\Phi }}_{+}\right(\nu )+{\mathit{\Phi }}_{-}(\nu \left)\right],$ (4.22)

where I is the $2N\times 2N$ identity matrix, M, P, A, B and W are $2N\times 2N$ matrices defined by

$\mathbf{M}=\text{diag}\left\{{\xi }_1,\dots ,{\xi }_N,{\xi }_1,\dots ,{\xi }_N\right\},$ (4.23)

$\mathbf{P}=\text{diag}\left\{P({\xi }_1:u),\dots ,P({\xi }_N:u),P({\xi }_1:u),\dots ,P({\xi }_N:u)\right\},$ (4.24)

$\mathbf{A}=\left[\begin{array}{cc}{\mathbf{a}}_{11}& {\mathbf{a}}_{12}\\ {\mathbf{a}}_{21}& {\mathbf{a}}_{22}\end{array}\right],$ (4.25)

$\mathbf{B}=\left[\begin{array}{cc}{\mathbf{b}}_{11}& {\mathbf{b}}_{12}\\ {\mathbf{b}}_{21}& {\mathbf{b}}_{22}\end{array}\right],$ (4.26)

where in Eqs. (4.25) and (4.26) the $N\times N$ submatrices are defined by

${\mathbf{a}}_{ij}=\text{diag}\left\{a_{ij}({\xi }_1:u),\dots ,a_{ij}({\xi }_N:u)\right\}$ (4.27)

and

${\mathbf{b}}_{ij}=\text{diag}\left\{b_{ij}({\xi }_1:u),\dots ,b_{ij}({\xi }_N:u)\right\},$ (4.28)

for $i,j=1,2$. Furthermore,

$\mathbf{W}=\left[\begin{array}{cc}\mathbf{w}& \mathbf{0}\\ \mathbf{0}& \mathbf{w}\end{array}\right]$ (4.29)

is the $2N\times 2N$ matrix with components

$[\mathbf{w}{]}_{i,j}=w_j,$ (4.30)

for $i,j=1,...,N$. Still, Φ _±(ν) are 2N vectors, such that

${\mathbf{\Phi }}_{\pm }\left(\nu \right)={\left[\begin{array}{cccccc}{\mathrm{\Phi }}_1(\nu ,\pm {\xi }_1)& \cdots & {\mathrm{\Phi }}_1(\nu ,\pm {\xi }_N)& {\mathrm{\Phi }}_2(\nu ,\pm {\xi }_1)& \cdots & {\mathrm{\Phi }}_2(\nu ,\pm {\xi }_N)\end{array}\right]}^T.$ (4.31)

Now we rewrite Eqs. (4.21) and (4.22) as the $4N\times 4N$ eigenvalue problem

$\left[\begin{array}{cc}\mathbf{C}& \mathbf{D}\\ \mathbf{E}& \mathbf{F}\end{array}\right]\left[\begin{array}{c}{\mathbf{\Phi }}_{+}\left(\mathrm{\nu }\right)\\ {\mathbf{\Phi }}_{-}\left(\mathrm{\nu }\right)\end{array}\right]=\lambda \left[\begin{array}{c}{\mathbf{\Phi }}_{+}\left(\nu \right)\\ {\mathbf{\Phi }}_{-}\left(\mathrm{\nu }\right)\end{array}\right],$ (4.32)

where the eigenvalues are given by

$\lambda ={\nu }^{-1},$ (4.33)

$\mathbf{C}={\mathbf{M}}^{\mathbf{-}\mathbf{1}}[\mathbf{I}-\mathbf{P}(\mathbf{A}+\mathbf{B}\left)\mathbf{W}\right],$ (4.34)

$\mathbf{D}=-{\mathbf{M}}^{\mathbf{-}\mathbf{1}}\mathbf{P}(\mathbf{A}+\mathbf{B})\mathbf{W},$ (4.35)

$\mathbf{E}={\mathbf{M}}^{\mathbf{-}\mathbf{1}}\mathbf{P}(\mathbf{A}-\mathbf{B})\mathbf{W}$ (4.36)

and

$\mathbf{F}=-{\mathbf{M}}^{\mathbf{-}\mathbf{1}}[\mathbf{I}-\mathbf{P}(\mathbf{A}-\mathbf{B}\left)\mathbf{W}\right].$ (4.37)

From Eq. (4.32) we obtain a set of 4N eigenvalues associated with their respective eigenvectors. How the eigenvalues are related with the separation constant in the form given by Eq. (4.33), we obtain, not only one, but 4N separation constants ν _j and 4N elementary solutions Φ _±(ν _j ) (see Eq. (4.31)). In this way, we are ready to write the general solution of the discrete-ordinates version of the G problem, given by Eq. (3.20), as

$\mathbf{G}(\tau ,\pm {\xi }_i)=\sum _{j=1}^{4N}{A}_j\mathbf{\Phi }({\nu }_j,\pm {\xi }_i){\text{e}}^{-\tau /{\nu }_j}.$ (4.38)

Although there are no mathematical evidences in the formulation presented until here, we observe in our numerical results of the computational procedure (presented in Section 6) that, for any N, all eigenvalues of the problem given by Eq. (4.32) are real.

As was mentioned by Arthur and Cercignani ², this problem have solution only for values of the downstream drift velocity $u^2<5/6$. This critical value of the parameter u (see Eq. (2.11 c)) is such that the corresponding Mach number is unity ². Still, this result was confirmed by Siewert and Thomas Jr. ¹⁹ and Sone and Sugimoto ²¹ calculations.

Continuing, since this is a conservative problem, we have to deal with the issue of having degenerate eigenvalues, that approach zero (separation constants going to infinity) as N tends to infinity ⁶. For the specific problem we are solving in this work, the number of eigenvalues with this behavior depends on the value of the parameter u. Thus, for $u=0$ and $u^2=5/6$ we find four degenerate eigenvalues and, for other values of u, we find three, when N increase. Because of that, we neglect these eigenvalues and look for exact solutions of the problem given by Eq. (3.20) to add to the general discrete ordinates solution. Therefore, Eq. (4.38) can be rewrite as

$\mathbf{G}(\tau ,\pm {\xi }_i)=A_1^{*}{\mathbf{G}}_1(\pm {\xi }_i)+A_2^{*}{\mathbf{G}}_2(\pm {\xi }_i)+A_3^{*}{\mathbf{G}}_3(\pm {\xi }_i)+A_4^{*}{\mathbf{G}}_4(\tau ,\pm {\xi }_i)+\sum _{j=1}^{4N-4}{A}_j\mathbf{\Phi }({\nu }_j,\pm {\xi }_i){\text{e}}^{-\tau /{\nu }_j},$ (4.39)

for $u=0$ and $u^2=5/6$, and

$\mathbf{G}(\tau ,\pm {\xi }_i)=A_1^{*}{\mathbf{G}}_1(\pm {\xi }_i)+A_2^{*}{\mathbf{G}}_2(\pm {\xi }_i)+A_3^{*}{\mathbf{G}}_3(\pm {\xi }_i)+\sum _{j=1}^{4N-3}{A}_j\mathbf{\Phi }({\nu }_j,\pm {\xi }_i){\text{e}}^{-\tau /{\nu }_j},$ (4.40)

for $0<u^2<5/6$ and $u^2>5/6$, where the introduced exact solutions are given by ¹⁹

${\mathbf{G}}_1\left(\xi \right)=\left[\begin{array}{c}1\\ 0\end{array}\right]{\text{e}}^{-(\xi -u{)}^2},$ (4.41 a)

${\mathbf{G}}_2\left(\xi \right)=\left[\begin{array}{c}\xi -u\\ 0\end{array}\right]{\text{e}}^{-(\xi -u{)}^2}$ (4.41 b)

and

${\mathbf{G}}_3\left(\xi \right)=\left[\begin{array}{c}(\xi -u{)}^2-1/2\\ 1\end{array}\right]{\text{e}}^{-(\xi -u{)}^2}.$ (4.41 c)

In Eq. (4.39), the fourth exact solution is ¹² (for $u=0$)

${\mathbf{G}}_4(\tau ,\xi )=\left[\begin{array}{c}\tau {\xi }^2-3\tau /2-{\xi }^3\\ \tau -\xi \end{array}\right]{\text{e}}^{-{\xi }^2},$ (4.41 d)

and (for $u^2=5/6$)

${\mathbf{G}}_4(\tau ,\xi )=\left[\begin{array}{c}\tau (\xi -u{)}^2-3u\tau (\xi -u)+\tau -(\xi -u{)}^3+3u(\xi -u{)}^2\\ \tau -(\xi -u)\end{array}\right]{\text{e}}^{-(\xi -u{)}^2}.$ (4.41 e)

The next step is to determine the arbitrary constants present in the solution (Eq. (4.39) or (4.40)). We use the boundary conditions for doing that.

Still in relation to the eigenvalue problem given by Eq. (4.32), the sign of eigenvalues is relevant. The number of eigenvalues (positives or negatives) also depend on the value of parameter u. For $u^2<5/6$ we find $2N-2$ positive eigenvalues and for $u^2\ge 5/6$ we obtain $2N-3$ eigenvalues with this sign. How the separation constants have the same sign of eigenvalues (see Eq. (4.33)), we then substitute the general solution, Eqs. (4.39) and (4.40), into Eq. (3.22) to obtain that the coefficients, $A_1^{*},A_2^{*},A_3^{*},A_4^{*}$ and A _j , these last associated with the negative separation constants, must be equal to zero. Then, for $u^2<5/6$ we write

$\mathbf{G}(\tau ,\pm {\xi }_i)=\sum _{j=1}^{2N-2}{A}_j\mathbf{\Phi }({\nu }_j,\pm {\xi }_i){\text{e}}^{-\tau /{\nu }_j}$ (4.42)

and for $u^2\ge 5/6$

$\mathbf{G}(\tau ,\pm {\xi }_i)=\sum _{j=1}^{2N-3}{A}_j\mathbf{\Phi }({\nu }_j,\pm {\xi }_i){\text{e}}^{-\tau /{\nu }_j},$ (4.43)

where ν _j are only the positive separation constants. In addition, the discrete-ordinates version of the interface boundary condition, Eq. (3.21), is

$\mathbf{G}(0,{\xi }_i)={\pi }^{1/2}\left[\Delta N_0{\mathbf{G}}_1\right({\xi }_i)+2(u_w-u\left){\mathbf{G}}_2\right({\xi }_i)+\Delta T_0{\mathbf{G}}_3({\xi }_i\left)\right],$ (4.44)

for $i=1,...,N$. In this way, if we substitute Eq. (4.42) into Eq. (4.44), we obtain for $u^2<5/6$ the square linear system $2N\times 2N$

$\sum _{j=1}^{2N-2}{A}_j\mathbf{\Phi }({\nu }_j,{\xi }_i)-{\pi }^{1/2}\left[\Delta N_0{\mathbf{G}}_1\right({\xi }_i)+\Delta T_0{\mathbf{G}}_3({\xi }_i\left)\right]=2{\pi }^{1/2}(u_w-u){\mathbf{G}}_2\left({\xi }_i\right),$ (4.45)

for $i=1,...,N$. If we substitute Eq. (4.43) into Eq. (4.44), we obtain for $u^2\ge 5/6$ the undetermined linear system $2N\times (2N-1)$

$\sum _{j=1}^{2N-3}{A}_j\mathbf{\Phi }({\nu }_j,{\xi }_i)-{\pi }^{1/2}\left[\Delta N_0{\mathbf{G}}_1\right({\xi }_i)+\Delta T_0{\mathbf{G}}_3({\xi }_i\left)\right]=2{\pi }^{1/2}(u_w-u){\mathbf{G}}_2\left({\xi }_i\right),$ (4.46)

for $i=1,...,N$. Once we solve Eqs. (4.45) and (4.46) we find the coefficients A _j and the quantities ∆N ₀ and ∆T ₀ defined in Eqs. (2.15b) and (2.15c).

Here, we note that, although, in principle, we can deal with the system given by Eq. (4.46), with the least squares method, we have not found convergence, when we increase N, to suggest a possible numerical reliable result. So, we consider this fact as an indication of the known theoretical result ² and assume Eq. (4.42) as our general solution for the strong evaporation problem, where the arbitrary constants are given by the solution of the square linear system defined by Eq. (4.45).

Thus, we substitute Eq. (4.42) into Eqs. (3.23) to (3.25) and we use Eqs. (4.3) and (4.4) to express the final form of the density perturbation

$\Delta N\left(\tau \right)={\pi }^{-1}\sum _{j=1}^{2N-2}{A}_j{\text{e}}^{-\tau /{\nu }_j}\mathbf{X}\left({\nu }_j\right),$ (4.47)

the velocity perturbation

$\Delta U\left(\tau \right)=-\Delta N\left(\tau \right)$ (4.48)

and the temperature perturbation

$\Delta T\left(\tau \right)=\frac{2}{3}{\pi }^{-1}\sum _{j=1}^{2N-2}{A}_j{\text{e}}^{-\tau /{\nu }_j}\mathbf{Y}\left({\nu }_j\right).$ (4.49)

For we obtain the expression for the heat flux, we multiply the first component of Eq. (4.2) by $(\xi -u{)}^2$ and we integrate the resulting equation over all ξ. In the same way, we integrate the second component of Eq. (4.2) over all ξ and, using Eqs. (4.3) and (4.4), we find

${\int }_{-\infty }^{\infty }{\left[\begin{array}{c}{\xi }^3\\ \xi \end{array}\right]}^T\mathbf{\Phi }(\nu ,\xi )\text{d}\xi =0.$ (4.50)

Thus, we substitute Eq. (4.42) into Eq. (3.26) and, using Eqs. (4.3), (4.4) and (4.50), we find for the heat flux

$Q_x\left(\tau \right)=-u{\pi }^{-1}\sum _{j=1}^{2N-2}{A}_j{\text{e}}^{-\tau /{\nu }_j}\left[\mathbf{Y}\right({\nu }_j)-\mathbf{X}({\nu }_j\left)\right].$ (4.51)

Still, in Eqs. (4.47) and (4.51)

$\mathbf{X}\left({\nu }_j\right)=\sum _{k=1}^N{w}_k{\left[\begin{array}{c}1\\ 0\end{array}\right]}^T\left[\mathbf{\Phi }\right({\nu }_j,{\xi }_k)+\mathbf{\Phi }({\nu }_j,-{\xi }_k\left)\right]$ (4.52)

and in Eqs. (4.49) and (4.51)

$\mathbf{Y}\left({\nu }_j\right)=\sum _{k=1}^N{w}_k{\left[\begin{array}{c}{u}^2-1/2\\ 1\end{array}\right]}^T\left[\mathbf{\Phi }\right({\nu }_j,{\xi }_k)+\mathbf{\Phi }({\nu }_j,-{\xi }_k\left)\right].$ (4.53)

5 A NONLINEAR APPROACH

Once the discrete-ordinates solution for the linearized version of the strong evaporation problem is completely established, we follow Scherer ^{12), (13} to define what we call a “post-processing (PP)” procedure. In this sense, we consider the proposed nonlinear model, given by Eq. (2.1) to (2.5), with boundary conditions defined in Eqs. (2.7) and (2.13), rewritten in terms of the dimensionless variables given in Eqs. (2.11). We then use the quantities evaluated by the ADO method, Eqs. (4.47) to (4.49), into Eq. (2.2), which defines the Maxwellian distribution. Then, we substitute this distribution into the right-hand side of Eq. (2.1), which is then solved for a known distribution $\varphi (\mathbf{r},\mathbf{v})$. The solution defines the original f distribution, which is then used to evaluate again Eqs. (2.3) to (2.6) - the macroscopic quantities for the gas. We present in the next section some numerical results and comparisons between this procedure and the linear version.

6 COMPUTATIONAL ASPECTS AND NUMERICAL RESULTS

To start the computational procedure, the first step is to define the quadrature scheme. Then, once we have the N quadrature points ξ _k and the weights w _k defined, the solution is concise and easy to implement. We proceed with:

• the solution of an eigenvalue problem, Eq. (4.32), to obtain the separation constants ν _j and the elementary solutions Φ _±(ν _j );

• the solution of a linear system, given by Eq. (4.45);

• the evaluation of the density, velocity and temperature perturbations (Eqs. (4.47) to (4.49)) and the heat flux (Eq. (4.51)). Still, from the solution of Eq. (4.45) we are able to get the quantities ∆N ₀ and ∆T ₀ in Eqs. (2.15b) and (2.15c).

• The quantities listed above are then used, in what we call “post-processing” procedure, in Eqs. (2.1) to (2.6).

We use the subroutine RG from the EISPACK package ²⁰ to solve the eigenvalue problem (Eq. (4.32)). To find the solution of the linear systems (Eqs. (4.45) and (4.46)) we then use the subroutines DGEFA and DGESL, for the determined case, and the subroutines DQRDC and DQRSL for the undetermined case, all from the LINPACK package ⁹.

The numerical results considered here were obtained by a FORTRAN program, using, in general, $N=80$ quadrature points. The computational time required for generating all quantities of interest for one value of u is a few seconds on a Pentium (R) Dual-Core (2.10 GHz, 3.0 GB RAM). If we increase N up to $N=200$, all digits listed in the tables are preserved (plus or minus one in the last digit). Still, in regard to our quadrature scheme, we define a half-range scheme in $[0,\infty )$. We use the transformation

$\gamma \left(\xi \right)={\text{e}}^{-\xi }$ (6.1)

to map the interval $[0,\infty )$ into $[0,1]$, where we are then able to make use of the usual Gauss-Legendre quadrature scheme, after using a new change of variable

$v\left(\gamma \right)=2\gamma -1.$ (6.2)

In regard to the numerical results, we first checked some previous results available in the literature ¹⁹, for the linearized problem, for the ratios $n_{\infty }/n_0$ and $T_{\infty }/T_0$. We obtained agreement with all digits listed in that reference. Results we generate from the linearized version are listed in Tables 1 and 2, referred as “Linear” case, where we also show the results we generate with the “Post-processing (PP)” approach. Both are compared with the results of the previous work ¹³ and with the Ytrehus ²⁵ results for a numerical treatment of a set of moment equations which take into account nonlinear terms. These results ²⁵ seem to be in agreement with experimental works. Looking at these tables it is apparent that the BGK model considered here is a significant improvement over the one-dimensional gas considered in Ref. ¹³.

Table 1:
Ratio $n_{\infty }/n_0$.

	Previous work 13		Present work		Ytrehus 25
u	Linear	PP	Linear	PP
0.0	1.00000	1.00000	1.00000	1.00000	1.000
0.1	8.81170(-1)	8.73324(-1)	8.61415(-1)	8.48948(-1)	8.494(-1)
0.2	7.96123(-1)	7.72378(-1)	7.65173(-1)	7.27601(-1)	7.283(-1)
0.3	7.32184(-1)	6.90327(-1)	6.94942(-1)	6.29386(-1)	6.303(-1)
0.4	6.82430(-1)	6.22488(-1)	6.41788(-1)	5.49278(-1)	5.501(-1)
0.5	6.42727(-1)	5.65542(-1)	6.00448(-1)	4.83468(-1)	4.841(-1)
0.6	6.10439(-1)	5.17090(-1)	5.67628(-1)	4.29068(-1)	4.292(-1)
0.7	5.83801(-1)	4.75381(-1)	5.41176(-1)	3.83870(-1)	3.834(-1)
0.8	5.61584(-1)	4.39135(-1)	5.19638(-1)	3.46183(-1)	3.447(-1)
0.9	5.42908(-1)	4.07414(-1)	5.02020(-1)	3.14713(-1)	3.120(-1)

Table 2:
Ratio $T_{\infty }/T_0$.

	Previous work 13		Present work		Ytrehus 25
u	Linear	PP	Linear	PP
0.0	1.00000	1.00000	1.00000	1.00000	1.000
0.1	9.19458(-1)	9.20433(-1)	9.55222(-1)	9.56753(-1)	9.567(-1)
0.2	8.42106(-1)	8.45652(-1)	9.10103(-1)	9.15865(-1)	9.152(-1)
0.3	7.69530(-1)	7.76876(-1)	8.64777(-1)	8.76913(-1)	8.756(-1)
0.4	7.02368(-1)	7.14499(-1)	6.41788(-1)	8.39659(-1)	8.378(-1)
0.5	6.40748(-1)	6.58453(-1)	7.74781(-1)	8.03915(-1)	8.016(-1)
0.6	5.84522(-1)	6.08403(-1)	7.30756(-1)	7.69501(-1)	7.671(-1)
0.7	5.33391(-1)	5.63854(-1)	6.87722(-1)	7.36245(-1)	7.342(-1)
0.8	4.86987(-1)	5.24234(-1)	6.45835(-1)	7.03974(-1)	7.028(-1)
0.9	4.44912(-1)	4.88933(-1)	6.05152(-1)	6.72524(-1)	6.729(-1)

We found more significant difference, when comparing results from the linearized model with the PP approach, for the density ratio than the temperature ratio. As expected, major variation is noted when u increases and the nonlinear modeling should be more effective. However, we also see that the PP procedure seems generate results in better agreement with Ytrehus ²⁵ results.

Still, in Tables 3 to 6, we list results for the distribution profiles, which were not provided in previous references ^{19), (25}, where only the ratios showed in Tables 1 and 2 were evaluated.

Table 3:
Density perturbation ∆N(τ).

	$u=0.5$		$u=0.8$		$u=0.9$
τ	Linear	PP	Linear	PP	Linear	PP
0.0	3.74170(-1)	4.078(-1)	8.08298(-1)	1.070	9.77886(-1)	1.521
0.1	3.10046(-1)	3.216(-1)	7.69429(-1)	9.258(-1)	9.72630(-1)	1.361
0.2	2.72970(-1)	2.763(-1)	7.42591(-1)	8.459(-1)	9.68785(-1)	1.277
0.3	2.44505(-1)	2.433(-1)	7.19507(-1)	7.869(-1)	9.65348(-1)	1.216
0.4	2.21154(-1)	2.174(-1)	6.98666(-1)	7.399(-1)	9.62139(-1)	1.170
0.5	2.01347(-1)	1.961(-1)	6.79404(-1)	7.008(-1)	9.59081(-1)	1.133
0.6	1.84195(-1)	1.782(-1)	6.61356(-1)	6.675(-1)	9.56131(-1)	1.102
0.7	1.69133(-1)	1.627(-1)	6.44292(-1)	6.384(-1)	9.53265(-1)	1.076
0.8	1.55770(-1)	1.492(-1)	6.28057(-1)	6.128(-1)	9.50465(-1)	1.054
0.9	1.43819(-1)	1.373(-1)	6.12538(-1)	5.898(-1)	9.47720(-1)	1.035
1.0	1.33065(-1)	1.268(-1)	5.97650(-1)	5.690(-1)	9.45020(-1)	1.018
2.0	6.57847(-2)	6.254(-2)	4.73550(-1)	4.274(-1)	9.19560(-1)	9.226(-1)
5.0	1.10003(-2)	1.073(-2)	2.47078(-1)	2.256(-1)	8.50492(-1)	8.224(-1)
7.0	3.74891(-3)	3.691(-3)	1.61927(-1)	1.513(-1)	8.07905(-1)	7.802(-1)

Table 4:
Velocity perturbation ∆U(τ).

	$u=0.5$		$u=0.8$		$u=0.9$
τ	Linear	PP	Linear	PP	Linear	PP
0.0	-3.74170(-1)	-2.897(-1)	-8.08298(-1)	-5.170(-1)	-9.77886(-1)	-6.034(-1)
0.1	-3.10046(-1)	-2.433(-1)	-7.69429(-1)	-4.807(-1)	-9.72630(-1)	-5.766(-1)
0.2	-2.72970(-1)	-2.164(-1)	-7.42591(-1)	-4.582(-1)	-9.68785(-1)	-5.608(-1)
0.3	-2.44505(-1)	-1.957(-1)	-7.19507(-1)	-4.403(-1)	-9.65348(-1)	-5.489(-1)
0.4	-2.21154(-1)	-1.786(-1)	-6.98666(-1)	-4.252(-1)	-9.62139(-1)	-5.392(-1)
0.5	-2.01347(-1)	-1.640(-1)	-6.79404(-1)	-4.120(-1)	-9.59081(-1)	-5.311(-1)
0.6	-1.84195(-1)	-1.512(-1)	-6.61356(-1)	-4.003(-1)	-9.56131(-1)	-5.243(-1)
0.7	-1.69133(-1)	-1.399(-1)	-6.44292(-1)	-3.896(-1)	-9.53265(-1)	-5.183(-1)
0.8	-1.55770(-1)	-1.298(-1)	-6.28057(-1)	-3.799(-1)	-9.50465(-1)	-5.132(-1)
0.9	-1.43819(-1)	-1.207(-1)	-6.12538(-1)	-3.709(-1)	-9.47720(-1)	-5.086(-1)
1.0	-1.33065(-1)	-1.125(-1)	-5.97650(-1)	-3.626(-1)	-9.45020(-1)	-5.046(-1)
2.0	-6.57847(-2)	-5.885(-2)	-4.73550(-1)	-2.994(-1)	-9.19560(-1)	-4.798(-1)
5.0	-1.10003(-2)	-1.062(-2)	-2.47078(-1)	-1.840(-1)	-8.50492(-1)	-4.512(-1)
7.0	-3.74891(-3)	-3.678(-3)	-1.61927(-1)	-1.314(-1)	-8.07905(-1)	-4.382(-1)

Table 5:
Temperature perturbation ∆T(τ).

	$u=0.5$		$u=0.8$		$u=0.9$
τ	Linear	PP	Linear	PP	Linear	PP
0.0	9.63884(-2)	7.013(-2)	4.27184(-1)	2.066(-1)	6.35434(-1)	2.861(-1)
0.1	8.48578(-2)	6.418(-2)	4.10018(-1)	2.156(-1)	6.32504(-1)	3.092(-1)
0.2	7.66260(-2)	5.927(-2)	3.97059(-1)	2.188(-1)	6.30202(-1)	3.221(-1)
0.3	6.98330(-2)	5.491(-2)	3.85566(-1)	2.200(-1)	6.28092(-1)	3.315(-1)
0.4	6.40061(-2)	5.098(-2)	3.75001(-1)	2.200(-1)	6.26094(-1)	3.387(-1)
0.5	5.89040(-2)	4.740(-2)	3.65116(-1)	2.191(-1)	6.24172(-1)	3.445(-1)
0.6	5.43773(-2)	4.413(-2)	3.55771(-1)	2.177(-1)	6.22306(-1)	3.492(-1)
0.7	5.03237(-2)	4.113(-2)	3.46876(-1)	2.159(-1)	6.20484(-1)	3.530(-1)
0.8	4.66683(-2)	3.837(-2)	3.38367(-1)	2.137(-1)	6.18697(-1)	3.561(-1)
0.9	4.33540(-2)	3.584(-2)	3.30199(-1)	2.113(-1)	6.16939(-1)	3.587(-1)
1.0	4.03359(-2)	3.349(-2)	3.22335(-1)	2.088(-1)	6.15206(-1)	3.609(-1)
2.0	2.06982(-2)	1.763(-2)	2.56089(-1)	1.797(-1)	5.98749(-1)	3.685(-1)
5.0	3.62361(-3)	3.154(-3)	1.33810(-1)	1.069(-1)	5.53819(-1)	3.534(-1)
7.0	1.24777(-3)	1.093(-3)	8.76890(-2)	7.444(-2)	5.26085(-1)	3.392(-1)

Table 6:
Heat flux Q _x (τ).

	$u=0.5$		$u=0.8$		$u=0.9$
τ	Linear	PP	Linear	PP	Linear	PP
0.0	1.14793(-1)	1.027(-1)	1.34017(-1)	3.025(-1)	2.22613(-2)	4.221(-1)
0.1	9.13798(-2)	8.094(-2)	1.23520(-1)	2.441(-1)	2.14856(-2)	3.438(-1)
0.2	7.90159(-2)	6.964(-2)	1.17601(-1)	2.114(-1)	2.11341(-2)	2.992(-1)
0.3	6.98781(-2)	6.147(-2)	1.12925(-1)	1.875(-1)	2.08887(-2)	2.661(-1)
0.4	6.25727(-2)	5.506(-2)	1.08931(-1)	1.687(-1)	2.06976(-2)	2.399(-1)
0.5	5.64955(-2)	4.981(-2)	1.05384(-1)	1.535(-1)	2.05399(-2)	2.186(-1)
0.6	5.13148(-2)	4.538(-2)	1.02159(-1)	1.409(-1)	2.04046(-2)	2.008(-1)
0.7	4.68241(-2)	4.158(-2)	9.91830(-2)	1.303(-1)	2.02854(-2)	1.858(-1)
0.8	4.28838(-2)	3.826(-2)	9.64052(-2)	1.213(-1)	2.01782(-2)	1.730(-1)
0.9	3.93941(-2)	3.533(-2)	9.37919(-2)	1.136(-1)	2.00804(-2)	1.620(-1)
1.0	3.62807(-2)	3.272(-2)	9.13181(-2)	1.068(-1)	1.99899(-2)	1.525(-1)
2.0	1.73687(-2)	1.663(-2)	7.15326(-2)	6.974(-2)	1.92932(-2)	1.022(-1)
5.0	2.78246(-3)	2.958(-3)	3.70904(-2)	3.632(-2)	1.77870(-2)	7.754(-2)
7.0	9.38625(-4)	1.020(-3)	2.43155(-2)	2.466(-2)	1.68996(-2)	7.650(-2)

In regard to parameter u _w (see Eq. (2.15 a)) associated to mass velocity at interface, how, in the begin of problem, the gas is in equilibrium with the phase boundary, in the same conditions of temperature and saturation pressure, and no there evaporation or condensation, we consider u _w equal to zero in our formulation. Thus, we set $u_w=0$ in our computational code to generate all numerical results presented. If a non zero value of u _w was used, the results would certainly change. Still, in Refs. ^{13), (19), (25}, which were used to make comparisons with numerical results, the mass velocity at interface also is considered as zero.

7 CONCLUDING REMARKS

The ADO method was used to develop a closed form solution for the nonlinear BGK version of the strong evaporation problem in rarefied gas dynamics, considering a monoatomic gas with three-dimensional velocity. The analytical discrete-ordinates solution obtained for the linearized version of the problem was associated with a re-evaluation of the quantities of interest, in order to take into account the nonlinear effects inherent to the problem. The new approach seemed to improve the results of the linearized version, mainly when the values of the drift velocity increase, when compared with results available in the literature. The BGK model considered here is also a significant improvement over the one-dimensional case considered in work ¹³.

REFERENCES

• ¹
  K. Aoki & N. Masukawa. Gas Flows Caused by Evaporation and Condensation on Two Parallel Condensed Phases and the Negative Temperature Gradient: Numerical Analysis by Using a Nonlinear Kinetic Equation. Phys. Fluids, 6 (1994), 1379-1395.
• ²
  M.D. Arthur & C. Cercignani. Non-Existence of a Steady Rarefied Supersonic Flow in a Half-Space. Z. Angew. Math. Phys., 31 (1980), 634-645.
• ³
  E. Aursand & T. Ytrehus. Comparison of Kinetic Theory Evaporation Models for Liquid Thin-Films. Int. J. Multiphase Flow, 116 (2019), 67-79.
• ⁴
  L.B. Barichello & C.E. Siewert. A Discrete-Ordinates Solution for a Non-Grey Model with Complete Frequency Redistribution. JQSRT, 62 (1999), 665-675.
• ⁵
  P.L. Bhatnagar, E.P. Gross & M. Krook. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev., 94 (1954), 511-525.
• ⁶
  K.M. Case & P.F. Zweifel. “Linear Transport Theory”. Addison-Wesley, Massachusetts (1967).
• ⁷
  C. Cercignani . “The Boltzmann Equation and its Applications”. Springer-Verlag, New York (1988).
• ⁸
  S. Chandrasekhar. “Radiative Transfer”. Dover, New York (1960).
• ⁹
  J.J. Dongarra, J.R. Bunch, C.B. Moler & G.W. Stewart. “LINPACK User’s Guide”. SIAM, Philadelphia (1979).
• ¹⁰
  A. Frezzotti. A Numerical Investigation of the Strong Evaporation of a Polyatomic Gas. Euro. J. Mechanics B/Fluids, 26 (2007), 93-104.
• ¹¹
  A. Frezzotti & P. Barbante. Kinetic Theory Aspects of Non-Equilibrium Liquid-Vapor Flows. Mech. Eng. Reviews, 4(2) (2017), 16-00540.
• ¹²
  C.S. Scherer. “Efeitos de Evaporação em Gases Rarefeitos”. Ph.D. thesis, Programa de Pós-Graduação em Engenharia Mecânica. Universidade Federal do Rio Grande do Sul, Porto Alegre, RS (2009).
• ¹³
  C.S. Scherer . An Analytical Approach to the Strong Evaporation Problem in Rarefied Gas Dynamics. Z. Angew. Math. Phys., 66 (2015), 1821-1833.
• ¹⁴
  C.S. Scherer & L.B. Barichello . An Analytical Approach to the Unified Solution of Kinetic Equations in Rarefied Gas Dynamics. III. Evaporation and Condensation Problems. Z. Angew. Math. Phys., 61 (2010), 95-117.
• ¹⁵
  C.S. Scherer , J.F. Prolo Filho & L.B. Barichello . An Analytical Approach to the Unified Solution of Kinetic Equations in Rarefied Gas Dynamics. I. Flow Problems. Z. Angew. Math. Phys., 60 (2009), 70-115.
• ¹⁶
  C.S. Scherer , J.F. Prolo Filho & L.B. Barichello . An Analytical Approach to the Unified Solution of Kinetic Equations in the Rarefied Gas Dynamics. II. Heat Transfer Problems. Z. Angew. Math. Phys. , 60 (2009), 651-687.
• ¹⁷
  C.E. Siewert . Heat Transfer and Evaporation/Condensation Problems Based on the Linearized Boltzmann Equation. Euro. J. Mechanics B/Fluids, 22 (2003), 391-408.
• ¹⁸
  C.E. Siewert & J.R. Thomas Jr. Strong Evaporation into a Half Space. Z. Angew. Math. Phys. , 32 (1981), 421-433.
• ¹⁹
  C.E. Siewert & J.R. Thomas Jr. Strong Evaporation into a Half Space. II. The Three-Dimensional BGK Model. Z. Angew. Math. Phys. , 33 (1982), 202-218.
• ²⁰
  B.T. Smith, J.M. Boyle, J.J. Dongarra , B.S. Garbow, Y. Ikebe, V.C. Klema & C.B. Moler . “Matrix Eigensystem Routines - EISPACK Guide”. Springer-Verlag, Berlin (1976).
• ²¹
  Y. Sone & H. Sugimoto. Strong Evaporation from a Plane Condensed Phase. In “Adiabatic Waves in Liquid-Vapor Systems”. Springer, Berlin (1990).
• ²²
  Y. Sone , S. Takata & F. Golse. Notes on the Boundary Conditions for Fluid-Dynamic Equations on the Interface of a Gas and its Condensed Phase. Phys. Fluids, 13 (2001), 324-334.
• ²³
  M.M.R. Williams. “Mathematical Methods in Particle Transport Theory”. Butterworth, London (1971).
• ²⁴
  S. Yasuda, S. Takata & K. Aoki . Evaporation and Condensation of a Binary Mixture of Vapors on a Plane Condensed Phase: Numerical Analysis of the Linearized Boltzmann Equation. Phys. Fluids , 17 (2005), 047105(1)-047105(19).
• ²⁵
  T. Ytrehus . Theory and Experiments on Gas Kinetics in Evaporation. In “10th International Symposium on Rarefied Gas Dynamics”. Aspen (1976), pp. 1197-1212.
• ²⁶
  Y.B. Zudin. “Non-Equilibrium Evaporation and Condensation Processes. Analytical Solutions”. Springer, Cham (2019).

Publication Dates

• Publication in this collection
  02 July 2021
• Date of issue
  May-Aug 2021

History

• Received
  15 Jan 2020
• Accepted
  06 Jan 2021