2.1. Enskog operator

The short-range interactions have two main influences on evolution of the velocity distribution: relaxation towards equilibrium and transport of conserved quantities. Note that the latter is only the case for dense flows. The Enskog collision operator derived by Enskog (Reference Enskog1922) and then modified by Van Beijeren & Ernst (Reference Van Beijeren and Ernst1973), Resibois (Reference Resibois1978) is an extension of the Boltzmann collision operator by including the volume of particles. The standard Enskog collision operator for hard-sphere collisions reads as

(2.4)\begin{align} S^{{Ensk}}(\mathcal{F}) &= \frac{1}{m} \int_{\mathbb{R}^{3}} \int_{0}^{2{\rm \pi} } \int_0^\sigma \left[ Y \left( {x}+\frac{1}{2} \sigma \hat{k}\right) \mathcal{F}( {v}^*, {x})\mathcal{F}( {v}^*_1,\ {x}+\sigma \hat{k})\right. \nonumber\\ &\quad -\left.Y \left( {x}-\frac{1}{2} \sigma \hat{k}\right)\mathcal{F}( {v}, {x})\mathcal{F}( {v}_1, {x}-\sigma \hat{k}) \right] g \mathcal{H}( g_i \hat{k}_i) \hat{b} \, \textrm{d}\hat{b} \, \textrm{d} \hat{\epsilon}\, \textrm{d} {v}_1, \end{align}

where $\hat {k}$ indicates the unit vector connecting the centre of a particle to its colliding pair with subscript $(\cdot )_1$, $\mathcal {H}(\cdot )$ is the Heaviside step function, $g$ is the relative velocity of colliding pair, and $\hat {b}$ and $\hat {\epsilon }$ are the impact parameter and scattering angle, respectively. Note that in the case of the hard-sphere potential the impact parameter relates to deflection angle $\chi$ via $\hat {b}=\sigma \cos (\chi /2)$, where $\chi \in [0,{\rm \pi} )$. For further details of the Enskog and Boltzmann collision operators, see, e.g. Chapman & Cowling (Reference Chapman and Cowling1970), Hirschfelder, Bird & Curtiss (Reference Hirschfelder, Bird and Curtiss1964) and Bird (Reference Bird1994).

In comparison to the dilute limit the Enskog collision operator admits a diameter distance between the colliding pair as well as an increase in collision rate by $Y$. The standard Enskog equation estimates the pair correlation factor $Y$ to be a function of density at the contact of colliding spheres, which can be chosen such that the correct equilibrium pressure is guaranteed. By means of the Carnahan–Starling equation of state, we get

(2.5)$$\begin{gather} Y = \frac{1}{2}\frac{2-\eta}{(1-\eta)^3}, \end{gather}$$

(2.6)$$\begin{gather}\textrm{with}\quad \eta := nb/4, \end{gather}$$

where $b:=2 {\rm \pi}\sigma ^3/3$ is the second virial coefficient (Carnahan & Starling Reference Carnahan and Starling1969). Furthermore, the number density is denoted by $n$ which gives the density $\rho =mn$.

In addition to the curse of dimensionality associated with $\mathcal {F}$, the increase of collision rate along with the non-locality of collision operator lead to further computational complications. Although direct Monte-Carlo particle methods can cope with the high dimensionality at the cost of introducing statistical noise, they inherit the issues with non-locality of the collision step along with resolving collisional scales.

2.2. Vlasov integral

By applying the mean-field limit, the attractive force of (2.3) is represented as the derivative of the Vlasov integral

(2.7)\begin{equation} \varPhi(x,t) = \int_{r>\sigma} \phi(r) n(y,t)\, {\textrm{d} y},\end{equation}

where $r:=|y-x|$ and $| \cdot |$ denotes Euclidean distance. Although the Vlasov integral can be computed using quadrature rules as discussed in Appendix A.2, long tails and mesh refinement near $r=\sigma$ make the numerical integration expensive in practice.

3.1. Short-range interactions

The computational cost associated with resolving the Enskog collision operator, i.e. right-hand side of (2.3), can be decoupled from collisional scales once the underlying jump process is described by an equivalent continuous process. Furthermore, if the stochastic paths generated by the latter process depends only on local moments of the distribution, efficient cell-based parallelization can be exploited.

Let us consider an Itô drift-diffusion process as a general prototype of such a continuous process. Then, the task becomes to find a drift-diffusion closure such that the resulting single-particle velocity distribution evolves consistent with the underlying jump process, up to a desired degree of moments. From a physical point of view, such a drift-diffusion process expresses all the forces acted on a sample particle into a portion which is correlated to the particle's state and a remainder being purely random. The former is governed by the drift whereas the latter is controlled by the diffusion coefficient.

Consider the position and velocity of particles given by random variables $(X(t),V(t))$. Suppose that the evolution of particles can be modelled through stochastic differential equations (SDEs) with drift $A$ and diffusion $D$ in velocity space, as well as a spatial drift $\hat {A}$ in the physical space

(3.1)

where $i\in \{1,2,3\}$. Here $\textrm {d} W_{t,i}:=W_i(t+\textrm {d} t)-W_i(t)$ indicates the increment of the Wiener process for the $i$th component of particle velocity which follows a Gaussian distribution with zero mean and variance of $\textrm {d} t$, i.e. $\textrm {d} W_{t,i}\sim \mathcal {N}(0,\textrm {d} t)$. Using Itô's lemma, evolution of the corresponding $\mathcal {F}$ follows a FP equation of the form

(3.2)\begin{equation} \dfrac{\partial \mathcal{F}}{\partial t} + \dfrac{\partial (\mathcal{F} v_i)}{\partial x_i} = S^{{DFP}}(\mathcal{F}),\end{equation}

where

(3.3)\begin{equation} S^{{DFP}} (\mathcal{F}) = \underbrace{-\frac{ \partial (\mathcal{F} A_i )}{\partial v_i} +\frac{1}{2}\frac{\partial^2 ( D^2 \mathcal{F}) }{\partial v_i \partial v_i}}_{\text{relaxation rates}} \underbrace{-\frac{ \partial (\mathcal{F} \hat{A}_i )}{\partial x_i}}_{\text{col. transfer}}.\end{equation}

For a relationship between the FP equation and SDEs, see, e.g. Oksendal (Reference Oksendal2013). The drift $A$ and diffusion $D$ allow us to fix the rate at which non-conservative velocity moments including the kinetic stress tensor

(3.4)\begin{equation} {\rm \pi}_{ij}^{kin} =\int_{\mathbb{R}^3} v_{\langle i}^\prime v_{j \rangle}^\prime \mathcal{F} \, \textrm{d} v \end{equation}

and the kinetic heat flux vector

(3.5)\begin{equation} q_{i}^{kin} = \tfrac{1}{2} \int_{\mathbb{R}^3} v_i^\prime v_j^\prime v_j^\prime \mathcal{F} \, \textrm{d} v \end{equation}

relax due to collisions. Note that subscript $\langle . \rangle$ denotes the deviatoric part of the tensor and $v^\prime :=v-U$ is the fluctuating velocity with $U=(1/\rho )\int v \mathcal {F} \, \textrm {d} v$ being the bulk velocity.

The spatial drift $\hat {A}$ accounts for the collisional transfer resulting from dense effects. From a physical point of view, the introduced drift in the position update corrects the particle position due to the volume exclusion. In other words, in the absence of such a correction the described SDEs only provide the velocity and position of particles with a vanishing diameter.

By considering the spatially homogeneous setting, we fix $A$ and $D$ such that the resulting moment system remains consistent with the Enskog equation up to the heat fluxes. Accordingly, the constraints

(3.6)\begin{gather} \int_{\mathbb{R}^{3}} v^\prime_{\langle i} v^\prime_{ j \rangle} S^{{DFP}}(\mathcal{F})\, \textrm{d} {v} = \int_{\mathbb{R}^{3}} v^\prime_{\langle i} v^\prime_{ j \rangle} S^{{Ensk}}(\mathcal{F})\, \textrm{d} {v} \quad \text{and} \end{gather}

(3.7)\begin{gather}\int_{\mathbb{R}^{3}} v^\prime_{ i} v^\prime_{ j} v^\prime_{ j} S^{{DFP}} (\mathcal{F})\, \textrm{d} {v} = \int_{\mathbb{R}^{3}} v^\prime_{ i} v^\prime_{ j} v^\prime_{ j} S^{{Ensk}} (\mathcal{F})\, \textrm{d} {v} \end{gather}

are devised. In order to obtain a tractable expression for the right-hand side of (3.6) and (3.7), here we approximate moments of the Enskog collision operator by applying Maxwell molecules. Therefore, we get

(3.8)\begin{gather} \int_{\mathbb{R}^3} \left( A_i v^\prime_j + A_j v^\prime_i - \frac{2}{3} A_k v^\prime_k \delta_{ij} + \frac{2}{3} \delta_{ij} D^2 \right) \mathcal{F} \, \textrm{d} v ={-} Y \frac{p^{{kin}}}{\mu^{{kin}} } {\rm \pi}_{ij}^{{kin}} \quad \textrm{and} \end{gather}

(3.9)\begin{gather}\int_{\mathbb{R}^3} ( A_i v^\prime_j v^\prime_j + 2 A_j v^\prime_j v^\prime_i ) \mathcal{F} \, \textrm{d} v ={-} Y \frac{2}{3} \frac{p^{{kin}}}{\mu^{{kin}} } q_i^{{kin}}. \end{gather}

Note that here the kinetic ideal gas viscosity coefficient is denoted by $\mu ^{{kin}}$, the ideal gas equilibrium pressure is $p^{{kin}}=nk_b T$ and the superscript $(\cdot )^{kin}$ indicates the kinetic contribution.

Next, the extra transport resulting from spatial dependencies of the Enskog operator is accounted by the spatial drift. The degrees of freedom in extra streaming induced by $\hat {A}$ are set such that the resulting FP model remains consistent with the Enskog collision operator up to the heat fluxes. Let $\psi =[1,v_i, v_jv_j/2]^\textrm {T}$, by the Taylor expansion of $\mathcal {F}(x\pm \sigma \hat {k})$ and $Y(x\pm \sigma \hat {k}/2)$ around $x$, the conservative moments of the Enskog collision operator can be formulated as fluxes

(3.10)\begin{gather} \int_{\mathbb{R}^3} \psi_l S^{{Ensk}}(\mathcal{F})\, \textrm{d} {v} ={-} \dfrac{\partial \varPsi^\phi_{lk}}{\partial x_k}, \quad \textrm{where} \end{gather}

(3.11) \begin{gather} \hspace{-4pc}\varPsi_{lk}^\phi = \frac{Y \sigma}{2m} \iiiint (\psi^*_l-\psi_l) \mathcal{F} \mathcal{F}_1 \hat{k}_k g \hat{b} \, \textrm{d} \hat{b} \, \textrm{d} \hat{\epsilon}\, \textrm{d} {v}_1\, \textrm{d} {v} \nonumber\\ \hspace{4.5pc} +\frac{Y \sigma^2}{4m} \iiint (\psi^*_l-\psi_l) \hat{k}_k \hat{k}_j \mathcal{F} \mathcal{F}_1 \frac{\partial}{\partial x_j} \ln \left(\frac{\mathcal{F}}{\mathcal{F}_1}\right) g \hat{b} \, \textrm{d} \hat{b} \, \textrm{d} \hat{\epsilon}\, \textrm{d} {v}_1\, \textrm{d} {v} \nonumber\\ \hspace{-13pc} + \cdots , \end{gather}

as shown by Chapman & Cowling (Reference Chapman and Cowling1970). In particular, the corrections accounted by the spatial drift $\hat {A}$ become

(3.12)\begin{equation} \frac{\partial }{\partial x_k} \int_{\mathbb{R}^{3}} \hat{A}_k \psi_l \mathcal{F} \, \textrm{d} v \approx \dfrac{\partial \varPsi^\phi_{lk}}{\partial x_k} \quad \textrm{for }l=1,\ldots,\dim( \psi). \end{equation}

Assuming Maxwellian for $\mathcal {F}$ in $\partial _{x_j} \ln ({\mathcal {F}}/{\mathcal {F}_1})$ inside the second integral of (3.11) and ignoring higher-order terms, a closure for collisional transfer can be derived,

(3.13)\begin{gather} \int_{\mathbb{R}^3} \hat{A}_i \mathcal{F} \, \textrm{d} v = 0, \end{gather}

(3.14)\begin{gather}\int_{\mathbb{R}^3} \hat{A}_i v_j \mathcal{F} \, \textrm{d} v = nbY(p^{kin}\delta_{ij}+2/5{\rm \pi}^{kin}_{ij}) -w\left(\frac{\partial U_k}{\partial x_k}\delta_{ij}+\frac{5}{6}\frac{\partial U_{\langle i}}{ \partial x_{j \rangle}}\right), \end{gather}

(3.15)\begin{gather}\text{and} \quad \frac12 \int_{\mathbb{R}^3} \hat{A}_i v_j v_j \mathcal{F} \, \textrm{d} v = \frac{3}{5} nbY q_i^{{kin}} - c_v w \frac{\partial T}{\partial x_i}. \end{gather}

Here $c_v=3 k_b /(2m)$ is the heat capacity at constant volume,

(3.16)\begin{equation} w = (nb)^2 Y \sqrt{mk_bT}/({\rm \pi}^{3/2} \sigma^2)\end{equation}

is the second viscosity coefficient, $T= \int _{\mathbb {R}^3} v^\prime _j v^\prime _j \mathcal {F} \, \textrm {d} v/(3nk_b)$ is the temperature and $k_b$ denotes the Boltzmann constant.

A simple set of polynomial ansatz to fulfil (3.8) and (3.9) in combination with (3.13)–(3.15) is adopted. Following Sadr & Gorji (Reference Sadr and Gorji2017), the cubic FP model for dense collisions

(3.17)\begin{gather} A_i=c_{ij}v^\prime_j+\gamma_i\left(v^\prime_jv^\prime_j-\frac{3k_bT}{m}\right)+\varLambda \left(v^\prime_iv^\prime_j v^\prime_j-\frac{2q^{kin}_i}{\rho}\right), \end{gather}

(3.18)\begin{gather}D =\sqrt{\frac{k_bT}{\tau m}}, \end{gather}

(3.19)\begin{gather}\textrm{and} \quad \hat{A}_i = \hat{c}_{ij} v^\prime_j + \hat{\gamma}_i \left(v^\prime_j v^\prime_j - \frac{3k_bT}{m}\right) + \hat{\varLambda} \left( v^\prime_i v^\prime_j v^\prime_j - \frac{2q^{kin}}{\rho} \right) \end{gather}

is devised. The time scale $\tau$ is set using kinetic scales $\tau =2\mu ^{{kin}}/(p^{{kin}}Y)$. The coefficients of the tensors $c$ and $\hat {c}$ along with the vectors $\gamma$ and $\hat {\gamma }$ need to be computed based on velocity moments. The cubic term is responsible to make sure that the SDEs remain stable (see Risken Reference Risken1989; Gorji, Torrilhon & Jenny Reference Gorji, Torrilhon and Jenny2011). The coefficients of the cubic terms are set to

(3.20)\begin{gather} \varLambda={-}\frac{|\det (u_{ij})|}{ (u^{(2)})^4 \tau }\quad \textrm{and} \end{gather}

(3.21)\begin{gather}\hat{\varLambda} ={-} \epsilon \frac{nbY}{k_b T/m }, \end{gather}

where $\epsilon =10^{-3}$, $\det (\cdot )$ indicates the determinant and

(3.22)\begin{equation} u_{i_1\ldots i_n}^{(k)}:=\frac{1}{\rho}\int_{\mathbb{R}^{3}}| {v^\prime}|^kv^\prime_{i_1}v^\prime_{i_2}\ldots v^\prime_{i_n}\mathcal{F}\, \textrm{d} v; \quad i_k\in\{1,2,3\}. \end{equation}

By substituting the ansatz for $A$ and $D$ in equations for the relaxation rates (3.8) and (3.9), one can show that

(3.23)\begin{align} &c_{ik}u^{(0)}_{kj}+c_{jk}u^{(0)}_{ki}+\gamma_iu^{(2)}_{j}+\gamma_ju^{(2)}_{i}={-}2\varLambda u^{(2)}_{ij} \quad \textrm{and} \end{align}

(3.24)\begin{align} &c_{ij}u^{(2)}_j+2c_{jk}u^{(0)}_{ijk}+\gamma_i(u^{(4)}-(u^{(2)})^2) +2\gamma_j(u^{(2)}_{ij}-u^{(2)}u_{ij}) \nonumber\\ &\quad ={-}\varLambda (3u^{(4)}_i-u^{(2)}_iu^{(2)}-2u^{(2)}_ju^{(0)}_{ij}) +\frac{5}{6}\frac{Yp^{kin}}{\mu^{{kin}}}q_i^{kin}. \end{align}

Equivalently, for $\hat {A}$, we get

(3.25)\begin{gather} \hspace{-9.5pc}\hat{c}_{jk}{\rm \pi}_{ik}^{kin}+\hat{c}_{ji}p^{kin}+2\hat{\gamma}_jq_i^{kin} ={-}\rho \hat{\varLambda} u^{(2)}_{ij}\\ + nbY(p^{kin}\delta_{ij}+2/5{\rm \pi}^{kin}_{ij}) -w\left(\frac{\partial U_k}{\partial x_k}\delta_{ij}+ \frac{5}{6}\frac{\partial U_{\langle i}}{ \partial x_{j \rangle}}\right) \quad \textrm{and} \end{gather}

(3.26)\begin{gather} \hat{c}_{ij} q_j^{kin} +\frac{1}{2}\rho \hat{\gamma}_i ( u^{(4)} - (u^{(2)}) ^2) =\frac{3}{5}nbYq_i^{kin} - w c_v\frac{\partial T}{\partial x_i} - \frac{1}{2}\rho \hat{\varLambda} (u^{(4)}_i - u^{(2)}_i u^{(2)} ),\end{gather}

where $i,j,k\in \{1,2,3\}$. The unknowns $c, \gamma , \hat {c}$ and $\hat {\gamma }$ can be found by solving two linear systems,

(3.27a,b)\begin{equation} A_v \begin{pmatrix} c_{1,1}\\ \vdots\\ c_{3,3} \\ \gamma_{1}\\ \vdots\\ \gamma_{3} \end{pmatrix} = b_v \quad \mathrm{and} \quad A_x \begin{pmatrix} \hat{c}_{11}\\ \vdots\\ \hat{c}_{3,3} \\ \hat{\gamma}_{1}\\ \vdots\\ \hat{\gamma}_{3} \end{pmatrix}= b_x,\end{equation}

where the matrices $A_v$ and $A_x$ together with the right-hand side matrices $b_v$ and $b_x$ are read from (3.23) and (3.24) and (3.25) and (3.26), respectively.

As shown in Sadr & Gorji (Reference Sadr and Gorji2017), the DFP model, besides reproducing near-equilibrium results, accurately captures non-equilibrium effects within the transition regime for moderate densities. Specifically, we expect a reasonable accuracy for the DFP model in the range of density where $nb\leq 1$.

3.2. Long-range interactions

The contribution of the long-range molecular potential in the evolution of molecular velocity distribution function can be accurately accounted by the Vlasov mean-field limit. Although the direct numerical evaluation of long-ranged forces through the Vlasov integral is grid-independent however, it can lead to high-computational costs due to a relatively large support of the kernel. In practice, global solutions over a spatial discretization turn out to be less expensive (Sadr & Gorji Reference Sadr and Gorji2019). Consider an approximation $\tilde {\phi }(r)$ to the attractive part of molecular potential $\phi (r)$ where

(3.28)\begin{gather} \tilde{\phi}( r) = a G( r) \hspace{0.5cm} \text{and} \end{gather}

(3.29)\begin{gather}G( r) = \frac{\textrm{e}^{-\lambda r}}{4 {\rm \pi}r}. \end{gather}

An estimate of fitting parameters $a$ and $\lambda$ can be obtained by solving an optimization problem of the type

(3.30)\begin{equation} (a,\lambda) = \underset{r \in (\sigma, \infty)}{\operatorname{argmin}}( \| \partial_r \phi(r)-\partial_r \tilde{\phi}(r) \|_2^2) \end{equation}

for a given molecular potential. Here, $||.||_2$ indicates the $L^2$ norm. This ansatz is motivated by the fact that $G(r)$ is the Green function of the SP equation

(3.31)\begin{equation} ({\rm \Delta}-\lambda^2)\tilde{\varPhi}_{r>0}( x,t)= n( x,t); \quad (\forall x\in \mathbb{R}^3), \end{equation}

in three-dimensional space. As shown above, the parameters $a$ and $\lambda$ give us flexibility to fit $\tilde {\phi }(r)$ to our desired potential. Hence, by substituting $\tilde {\phi }(r)$ in the Vlasov integral, we obtain

(3.32)\begin{align} \varPhi( x, t) &= \int_{r>\sigma}\phi( r) n( y, t)\, {\textrm{d} y} \nonumber\\ &\approx a\underbrace{\int_{r>0}G( r) n( y, t)\, {\textrm{d} y}}_{\tilde{\varPhi}_{r>0}} - a\underbrace{\int_{r<\sigma}G( r) n( y, t)\, {\textrm{d} y}}_{\tilde{\varPhi}_{r<\sigma}} \nonumber\\ &= a( \tilde{\varPhi}_{r>0} - \tilde{\varPhi}_{r<\sigma} ). \end{align}

Furthermore, by assuming that the number density is smooth within $r\in (0,\sigma )$, an explicit expression for $\tilde {\varPhi }_{r<\sigma }$ can be derived. In particular, using change of variables

(3.33)\begin{align} \tilde{\varPhi}_{r<\sigma} (x,t) &= \int_{r<\sigma} G(| y - x|) n( y, t) \, \textrm{d}^3 y \nonumber\\ &= \int_{| y|<\sigma} G(| y| ) n( y- x, t) \, \textrm{d}^3 y, \end{align}

and deploying Taylor's expansion for $n(y-x,t)$ around $(x,t)$, it can be shown that

(3.34)\begin{align} \tilde{\varPhi}_{r<\sigma} (x,t)&= n( x,t) \underbrace{\int_{| y|<\sigma} G(| y| )\, \textrm{d}^3 y}_{\text{analytical expr.}} + \frac{1}{2}\frac{\partial^2 n ( x,t)}{\partial x_j \partial x_j} \underbrace{\int_{| y|<\sigma} G(| y|) y_j y_j \, \textrm{d}^3 y}_{\text{analytical expr.}} +\cdots \nonumber\\ &\approx n(x,t) [\textrm{e}^{- \lambda \sigma} {(- \lambda \sigma + \textrm{e}^{\lambda \sigma} - 1) }/{\lambda^{2}}] \nonumber\\ &\quad + \frac{\partial^2 n(x,t)}{\partial x_j \partial x_j} [ \textrm{e}^{- \lambda \sigma} {(- \lambda^{3} \sigma^{3}/6 - \lambda^{2} \sigma^{2}/2 - \lambda \sigma + \textrm{e}^{\lambda \sigma} - 1) }/{\lambda^{4}}], \end{align}

by second-order truncation.

Therefore, instead of solving the Vlasov integral (2.7) directly, we find solution of the SP equation globally where the evaluated number density is fed as the source term. Fast Poisson solvers can be adjusted to the SP problem in the bounded domain $\varOmega _x \subset \mathbb {R}^{3}$ with Dirichlet boundary conditions that provide the unique solution. The values of the potential on the boundaries of the domain $\partial \varOmega _x$ are computed using the kernel $G(r)$ along with the boundary condition of $\mathcal {F}$. The detailed numerics of the SP model is presented § 4.2.

As a simpler alternative to the SP model, one can obtain an explicit approximation of the Vlasov integral by expanding the number density around the point of interest and truncating the higher-order terms. Although simple expressions can be obtained, this approach requires strong regularity assumptions behind the number density and also introduces further error as long-range contributions of the potential are computed locally; see Sadr & Gorji (Reference Sadr and Gorji2019).

4.1. Stochastic particle method

Since the SDEs (3.1a) and (3.1b) do not admit closed-form time integration (due to nonlinearity of the drift), a numerical scheme has to be adopted. As explained by Jenny et al. (Reference Jenny, Torrilhon and Heinz2010), Gorji et al. (Reference Gorji, Torrilhon and Jenny2011), Gorji & Jenny (Reference Gorji and Jenny2014), the simple Euler–Maruyama scheme does not guarantee energy conservation on average. Given the fact that the Langevin equation with linear drift allows for explicit time integration, a decomposition of the SDEs into a linear drift portion and a remainder is pursued.

Following Sadr & Gorji (Reference Sadr and Gorji2017) we consider a spatial discretization where each cell $c$ contains $N_{p,{(\textit {c})}}$ particles. Then the evolution of each particle is described via an explicit time-stepping scheme, where the solution at time step $(n+1)$ is computed based on the observables estimated at the $(n)$th step.

Realizations of the SDEs (3.1a) and (3.1b), through particles with weights $w^{(k)}$, velocities $V^{(k)}$ and positions $X^{(k)}$ with $k\in \{1,\ldots ,N_p\}$, provide an estimate

(4.1)\begin{equation} \mathcal{F}(v,x,t)\approx \sum_{k=1}^{N_p}w^{(k)} \delta ( V^{(k)}(t)- v)\delta( X^{(k)}(t)- x) \end{equation}

of the distribution which can be used to obtain various observables. By generalizing particle discretization introduced by Gorji & Jenny (Reference Gorji and Jenny2014), the velocity and position of each particle gets updated as

(4.2) \begin{equation} V_i^{(*)}= \langle V_i \rangle^{(n)}_{(\textit{c})} + \alpha_{(\textit{c})}\tilde{V}^{(n)}_i, \\ X_i^{(n+1)}= X_i^{(n)} + V_i^{(n+1)}{\rm \Delta} t+\delta \tilde{X}^{(n)}_i, \end{equation}

where

(4.3)\begin{align} \tilde{V}^{(n)}_i&= V_i^{\prime, (n)} \,\textrm{e}^{-{\rm \Delta} tY^{(n)}/\tau^{(n)}} + c^{(n)}_{ij} V_j^{\prime, (n)}+ \gamma^{(n)}_i ( V_j^{\prime, (n)} V_j^{\prime, (n)} - 3k_bT^{(n)}/m ) \nonumber\\ &\quad + \varLambda ( V^{\prime, (n)}_i V^{\prime, (n)}_j V^{\prime, (n)}_j - 2q^{(n)}_i/\rho^{(n)}) + \sqrt{\frac{kT^{(n)}Y^{(n)}}{\tau^{(n)} m}({1-\textrm{e}^{{-2{\rm \Delta} tY^{(n)}}/{ \tau^{(n)}}}})} \xi_i, \end{align}

and

(4.4)\begin{align} \delta \tilde{X}^{(n)}_i&= ( \hat{c}^{(n)}_{ij} V^{\prime, (n)}_j + \hat{\gamma}^{(n)}_i (V^{\prime, (n)}_j V^{\prime, (n)}_j - {3k_bT^{(n)}}/{m}) \nonumber\\ &\quad + \hat{\varLambda}^{(n)} ( V^{\prime, (n)}_i V^{\prime, (n)}_j V^{\prime, (n)}_j - 2q^{(n)}_i/\rho^{(n)} ) {\rm \Delta} t. \end{align}

Notice that $V^{(*)}$ is modified accordingly once external or long-range forces are considered; otherwise, we have $V^{(n+1)} = V^{(*)}$. Here $\xi _i$ is a random number drawn from the standard normal distribution, i.e. $\xi _i \sim \mathcal {N}(0,1)$, and $\alpha ^{(\textit {c})}$ is the scaling factor which guarantees conservation of kinetic energy at each cell for a given time step,

(4.5)\begin{equation} \alpha_{(\textit{c})} =\frac{\langle V_{j}^{\prime, (n)} V_{j}^{\prime, (n)} \rangle_{(\textit{c})} } {\langle \tilde{V}_j^{(n)} \tilde{V}_{j}^{(n)} \rangle_{(\textit{c})}}.\end{equation}

Note that $\langle \cdot \rangle ^{(\textit {c})}$ indicates a moment estimator. Using (4.1), we observe that

(4.6)\begin{equation} \langle \hat{\phi} \rangle_{(\textit{c})} = \frac{1}{\hat{\rho}} \sum_{i\in \mathcal{S}_c} w^{(i)}{\hat{\phi}}^{(i)},\end{equation}

where ${\hat {\phi }}^{(i)}$ indicates $\phi$ evaluated by the $i$th particle, $\mathcal {S}_c$ is a set of particles in the cell $(\textit {c})$ and

(4.7)\begin{equation} \hat{\rho}= \frac{1}{\delta V_c}\sum_{i\in \mathcal{S}_c}w^{(i)} \end{equation}

is the estimated density based on the cell volume $\delta V_c$.

In comparison to conventional Langevin simulations, here the macroscopic coefficients $\{ c, \gamma , { \hat {c}}, { \hat {\gamma }}\}$ which close the drift term are computed from linear systems (3.27a,b) for each computational cell and per each time step. Hence, correct relaxation rates in the velocity update besides accurate collisional transfer in the position update are obtained.

Algorithm 1

DFP-SP solution algorithm

4.2. Discretization of the screened-Poisson model

In order to translate the long-range interactions into a suitable elliptic problem we first need to represent (3.31) in a bounded domain $\varOmega _x \subset \mathbb {R}^3$. Following uniqueness of the solution, the SP equation can be reset in a bounded domain $\varOmega _x$ as

(4.8)\begin{equation} \left\{\begin{array}{@{}ll} ({\rm \Delta} -\lambda^2) \tilde{\varPhi}_{r>0} ={-}n & \text{in }\varOmega_x, \\ \tilde{\varPhi}_{r>0} = \tilde{\varPhi}_D & \text{on }\partial \varOmega_x, \end{array}\right.\end{equation}

where the Dirichlet condition $\tilde {\varPhi }_D$ is evaluated via

(4.9)\begin{equation} \tilde{\varPhi}_D(x) = \int_{r>0} G(r) n(y)\, {\textrm{d} y}, \quad \forall\, x\in\partial \varOmega_x. \end{equation}

The boundary term is computed numerically using quadrature estimation of the Vlasov integral and introducing a cut-off to contain the computational cost, i.e. $r\in (0, r_{{cut}})$. Note that the singularity in $G(r)$ is avoided once the integral is formulated in the spherical coordinates. Here, we have set $r_{{cut}} = 3 \sigma$ and used $N_{{quad}}$ Legendre–Gauss quadrature points for the numerical integration. Once the values of the unique global solution on the boundary $\partial \varOmega _x$ is computed, this elliptic differential equation on the bounded domains can be solved for the interior of the domain using standard numerical methods. In Appendix B we present a short description for the adopted hybridized discontinuous Galerkin solver. Once our estimate of the attractive force at time step $(n)$ is obtained, the velocity update takes the form

(4.10)\begin{equation} V^{(n+1)}= V^{(*)}+\hat{H}^{(n)}(X^{(n)}){\rm \Delta} t,\end{equation}

where $\hat {H}^{(n)}(X^{(n)})$ is the approximated attractive force at the particle position $X^{(n)}$ (see Appendix B for details).

4.3. Solution algorithm

It is evident from the time-stepping scheme given by (4.2a) and (4.2b) that the evolution of the position and velocity of particles is decoupled during each time step. While more accurate schemes such as those discussed by Jenny et al. (Reference Jenny, Torrilhon and Heinz2010), Gorji et al. (Reference Gorji, Torrilhon and Jenny2011), Jenny & Gorji (Reference Jenny and Gorji2019) can be devised for a joint position-velocity update, here for simplicity, the decoupled version is utilized. Note that, as shown by Pfeiffer & Gorji (Reference Pfeiffer and Gorji2017), Jun et al. (Reference Jun, Pfeiffer, Mieussens and Gorji2019), this decoupling is accurate and justified as long as the gradients of macroscopic fields are resolved. In comparison to DSMC (Bird Reference Bird1994), where the collisional scales have to be resolved, here we obtain a significant computational advantage by relaxing the discretization constraints.

Algorithm 1 provides a summary of the DFP-SP solution method for simulations of multi-phase flows. Prior to the simulation, the domain is discretized into a set of computational cells. Furthermore, an appropriate time step size and particle weight are chosen. After initializing the particles according to their weights and initial density field, the particles evolve according to algorithm 1. Note that the particles boundary conditions are implemented similar to DSMC, where depending on the distribution of incoming or reflected particles, open or wall boundaries are handled, respectively.

5.1. Isothermal co-existing states

In this section we investigate the accuracy of the DFP-SP model in predicting the co-existing liquid–vapour states at constant temperatures against ESMC-Vlasov simulation results, besides the Maxwell area reconstruction solution (Li, Luo & Li Reference Li, Luo and Li2013). We impose a thermostat every $N_{therm}$ steps until the solution reaches stationary state, in order to keep the temperature of the system constant (Frezzotti et al. Reference Frezzotti, Gibelli and Lorenzani2005). The thermostat is then removed and afterwards the system reaches the equilibrium. The thermostat phase provides us with a decent initial distribution of particles in physical space for a given temperature. Furthermore, note that the solution after the thermostat phase is already very close to the steady-state equilibrium solution.

In order to achieve simulation results of co-existing states, we consider a liquid slab $[0,1] \times [0,40\sigma ] \times [0,1]\ \mathrm {m}^3$ with number density $n_l^{(0)}$ inside vapour with number density $n_v^{(0)}$. Here, the total volume $[0,1] \times [0,L] \times [0,1]\ \mathrm {m}^3$, with $L=100 \sigma$, is enclosed in a box with specular walls. Due to symmetry of the solution space, we only need to simulate the particle motion in the $x_2$ direction, yet velocity is treated in three-dimensional space. The spatial dimension is decomposed to $100$ cells and ${\rm \Delta} t$ was chosen to satisfy a CFL number of $0.01$. Initially, in each cell $N_{p/c}$ particles per cell with a Gaussian distribution in velocities, i.e. $V_i \sim \mathcal {N}(v_i|0,3k_bT^{(0)}/m)$, and a uniform distribution in position are sampled. The statistical weight of particles $w$ was taken constant and chosen such that $N_{p/c}\geq 300$ in the vapour phase. The velocity of particles are thermostated every $N_{therm}=10$ steps for the first $1000$ steps. Then, the simulation is continued without the thermostat until the stationary state is achieved ($\sim 2000$ steps). Then, the solution over next $3000$ steps is averaged for the post-processing.

Note that the Maxwell solution leads to a minimization problem, i.e. find the pressure $p_0 = {\operatorname {argmin}} ( \int _{n_v}^{n_l} p(n,T)-p_0\, \textrm {d} n)$ for a given temperature. This pressure leads to two equal separated areas enclosed between $p(n,T)$ and the constant $p(n_l,T)=p(n_v,T)=p_0$. In order to compute the integral, one needs to find $n_l$ and $n_v$ for any guessed pressure $p_0$ which leads to finding roots of a nonlinear equation $p(n,T)=p_0$ for $n$. In order to simplify computation of the Maxwell solution, we approximate the pair correlation function with a polynomial $Y(n)\approx 1 + 0.625n b + 0.2869 (n b)^2 + 0.115 (n b)^3+\cdots$, instead of (2.6); see Chapman & Cowling (Reference Chapman and Cowling1970). As shown in figure 1, the polynomial approximation has a reasonable accuracy for the range of $n \sigma ^3\leq 0.45$. Hence, we consider the equilibrium pressure

(5.3)\begin{gather} p_{eq}^{Max., Suth.} = nk_bT(1+nb\tilde{Y}(n)) - \frac{4 {\rm \pi}}{3} \phi_0 \sigma^2 n^2, \end{gather}

(5.4)\begin{gather}p_{eq}^{Max., SP} = nk_bT(1+nb\tilde{Y}(n)) + {a}\, {\textrm{e}^{-\lambda \sigma} } (\lambda^2 \sigma^2 + 3 \sigma \lambda +3 ) n^2/(6 \lambda^2), \quad \mathrm{where} \end{gather}

(5.5)\begin{gather}\tilde{Y}(n) = 1 + 0.625n b + 0.2869 (n b)^2 + 0.115 (n b)^3, \end{gather}

for the Maxwell area reconstruction solution of co-existing states with validity range of $0< n \sigma ^3\leq 0.45$.

Figure 1.

Pair correlation function $Y(n)$ and polynomial approximation $\tilde {Y}(n)$ used in the Maxwell solution against normalized density for $n\sigma ^3\in (0,0.6)$.

Here, we consider several equilibrium states within $T \in [130,185]\ \mathrm {K}$. The initial number density of vapour and liquid for each case is provided in table 1. The density and temperature profiles for each case are shown in figures 2–4, respectively, where a reasonable consistency between DFP-SP and ESMC-Vlasov solutions is observed. The obtained $(n,T)$ pair for liquid and vapour of each case is plotted against the Maxwell area reconstruction solution in figure 5, where a reasonable agreement between the kinetic solution and Maxwell area reconstruction can be seen for its range of validity.

Table 1.

Initial number density and temperature of liquid slab and vapour for the co-existing state simulations.

Figure 2.

Number density profiles for simulation of the DFP-SP model (in red) and ESMC-Vlasov solution algorithm (in blue) for the obtained isothermal equilibrium simulation in (a) $T = 180\ \mathrm {K}$, (b) $T = 170\ \mathrm {K}$, (c) $T = 160\ \mathrm {K}$ and (d) $T = 150\ \mathrm {K}$.

Figure 3.

Number density profiles for simulation of the DFP-SP model (in red) and ESMC-Vlasov solution algorithm (in blue) for the obtained isothermal equilibrium simulation with (a) $T=140\ \mathrm {K}$, (b) $T=140\ \mathrm {K}$, log-scale, (c) $T=130\ \mathrm {K}$ and (d) $T=130\ \mathrm {K}$, log-scale.

Figure 4.

Temperature profiles obtained from simulation of the DFP-SP model (a) and ESMC-Vlasov solution algorithm (b) for co-existing state test case at temperatures $T\in \{130, 140,150,160,165,170,175,180,185\}\ \mathrm {K}$ shown in rainbow spectrum, from dark blue to red.

5.2. Evaporation near equilibrium

In order to present the DFP-SP prediction of the surface tension, evaporation of a liquid argon slab in contact with the vapour phase near equilibrium is studied. In particular, we consider the phase transition of the liquid slab located within $[0,1]\times [0,30 \sigma ] \times [0,1]\ \textrm {m}^3$ which is surrounded by the vapour confined in $[0,1]\times [0,L] \times [0,1]\ \textrm {m}^3$, with $L=90 \sigma$. Since only the $x_2$ direction is relevant, the other dimensions in the physical space are ignored in the evolution of particles positions. The liquid slab and the vapour around it are initialized with number densities $n_{{liq}}$ and $n_{{vap}}$, respectively. In order to maintain the stationary scenario, the liquid slab is heated up every $N_{{therm}}$ steps, where the particles inside $x_2\in [L/2-L_{{therm}}/2,L/2+L_{{therm}}/2]$ are thermostated with the initial temperature $T^{(0)}$. This is achieved by rescaling the particles fluctuating velocities. Furthermore, the specular reflection boundary conditions are deployed here.

After conducting convergence studies, it is found that $N_{{cell}}=180$ computational cells and the CFL number $0.01$ provide converged results. Moreover, a fixed statistical weight of $w=10^{14}$ was picked for the test cases with $T^{(0)} \in \{150, 160, 170, 180 \}\ \mathrm {K}$, $n_{{liq}}^{(0)}\in \{0.6,0.526,0.45,0.36 \} \times \sigma ^{-3} \ \mathrm {m}^{-3}$ and $n_{{vap}}^{(0)} \in \{ 0.05,0.077,0.12,0.15 \} \times \sigma ^{-3} \ \mathrm {m}^{-3}$. The stationary distribution is obtained after $5000$ steps and the time averaging is performed over $15\,000$ steps onward. Here, $L_{{therm}}=10 \sigma$ was chosen as the length of the thermostated sub-domain which was carried out once every $N_{{therm}}=1000$ steps. The initial densities of vapour and liquids for each simulation were picked such that unstable regions of the phase diagram are not met, as depicted in figure 6. Here, the total equilibrium pressure is $p_{{eq}}^{{tot}}(n,T) :=nk_bT(1+nbY)+p^{{att}}$, where the attractive part

(5.6)\begin{equation} p^{att}_{{eq}}(n)={-}({2{\rm \pi} }/{3})n^2 \int_0^\infty r^3 \frac{\partial \phi(r)}{\partial r} \, \textrm{d} r \end{equation}

becomes $p^{att}_{{eq}}(n)=-({4 {\rm \pi}}/{3}) \phi _0 \sigma ^3 n^2$ and $p^{{att}}_{{eq}} (n) = {a } { \textrm {e}^{-\lambda \sigma } } (\lambda ^2 \sigma ^2 + 3 \sigma \lambda +3 ) n^2/(6 \lambda ^2)$ in the case of Sutherland and screened-Poisson potentials, respectively (Sadr & Gorji Reference Sadr and Gorji2019).

For the liquid slab at equilibrium, the surface tension coefficient $\gamma$ is defined as a norm of the difference between the diagonal element of the kinetic pressure tensor and the equilibrium kinetic pressure across the phase transition,

(5.7)\begin{equation} \gamma = \tfrac{1}{2} \int_0^L (p_{22}-p_0) \, \textrm{d} x_2. \end{equation}

For details, see Ghoufi & Malfreyt (Reference Ghoufi and Malfreyt2010) and Hirschfelder et al. (Reference Hirschfelder, Bird and Curtiss1964). The integration over the domain can be taken numerically, since the discrete values of the pressure tensor are estimated at every cell.

Comparison between simulation results obtained from the DFP-SP model with respect to the benchmark ESMC-Vlasov is shown in figure 7. Although the DFP-SP model underpredicts the peak in estimating $|p_{22}-p_0|$ at the interface, overall a very good agreement between the two is observed. Note that the discrepancy in $\gamma$ between both methods increases as the liquid gets denser at lower temperatures. Furthermore, the computational costs of the discussed simulation settings are presented in figure 8. Accurate quantification of the relative computational costs would require a well-designed prospective study with detailed numerical investigations. However, a general comparison is possible by comparing the complexity order of ESMC-Vlasov and DFP-SP methods. For the former, observe that the number of operations follows $\mathcal {O}(n\sqrt {T})$. Yet the cost of the latter is independent of the density, which makes our assembled approach attractive for efficient Monte-Carlo simulations of multi-phase flows.

Figure 5.

Co-existing states obtained from simulations of the DFP-SP model (in red), ESMC-Vlasov solution algorithm (in blue), and Maxwell area reconstruction solution of Sutherland (in black) and screened-Poisson potential (in green) for equilibrium at $T\in \{130, 140,150,160,165,170,175,180,185\} \ \mathrm {K}$ with relative error in density of vapour $|n_v^{DFP-SP}-n_v^{ESMC-Vlasov}|/n_v^{ESMC-Vlasov}={0.052, 0.037,0.002,0.004,0.06,0.039,0.001,0.05,0.01}$, respectively. Density for the Monte Carlo solution is obtained by averaging over cells within $x_2\in (0,5\sigma )\cup (L-5\sigma ,L)$ and $x_2\in (L/2-5\sigma , L/2+5\sigma )$ for vapour and liquid, respectively. (a) Co-existing states for $n_l \sigma ^3\leq 0.55$ in log-scale and (b) co-existing states extended.

Figure 6.

Contours of total equilibrium pressure $p_{{eq}}^{{tot}}(n,T){:=nk_bT(1+nbY)+p^{{att}}_{{eq}}}$ at constant temperatures as well as the initial state pairs $(n_{{liq}},T_{{liq}})$ and $(n_{{vap}},T_{{vap}})$ for the evaporation test case with initial temperatures of $T^{(0)}\in \{ 150,160,170,180 \} \ \mathrm {K}$ and corresponding vapour (in red) and liquid (in blue) number densities of $n_{{liq}}^{(0)}=\{0.6,0.526,0.45,0.36 \} \times \sigma ^{-3} \ \mathrm {m}^{-3}$ and $n_{{vap}}^{(0)}=\{ 0.05,0.077,0.12,0.15 \} \times \sigma ^{-3} \ \mathrm {m}^{-3}$, respectively. Note that here the total pressure with attractive part associated with screen-Poisson's potential is used, $p^{att}_{{eq}}(n)={{a } { \textrm {e}^{-\lambda \sigma } } (\lambda ^2 \sigma ^2 + 3 \sigma \lambda +3 ) n^2/(6 \lambda ^2)}$.

Figure 7.

Phase transition of argon at $T^{(0)}=160\ \mathrm {K}$. (a) Normalized number density. (b) Normalized temperature. (c) Normalized difference between diagonal and equilibrium kinetic pressure. (d) Normalized surface tension for liquid temperatures of $T^{(0)}\in \{ 150,160, 170, 180 \}\ \mathrm {K}$. Here, a maximum relative error of 7% in evaluating $\gamma$ is observed.

Figure 8.

Execution time per particle $\tau :=t^{{ex.}}/N_p$ for evaporation simulation at the equilibrium using ESMC-Vlasov and DFP-SP solution algorithms. The cost is normalized by $\tau _0$ denoting DFP-SP execution time for the case $T^{(0)}=180$ K.

Next, we compare the obtained numerical results in estimating the surface tension coefficient with respect to an analytical approximation based on the free surface energy $S(T)$. Following Elsner (Reference Elsner1991),

(5.8)\begin{equation} \gamma(T) = T \int_T^{T_{crit}} \frac{S(T)}{T^2} \, \textrm{d} T, \end{equation}

where

(5.9)\begin{equation} S(T) = \frac{\alpha}{3} \int_{\vartheta _v}^{\vartheta _l} \left(\frac{\mu}{\vartheta }-p_v\right)\frac{\textrm{d} \vartheta }{\vartheta^{2/3}}. \end{equation}

Here, $\vartheta =1/n$ denotes the reciprocal number density, $\mu$ is the chemical potential, $p_v$ is the vapour pressure and $\alpha$ a proportionality factor. Note that $\mu -p_v \vartheta$ gives the Helmholtz free energy per particle. For simplicity, similar to Hadjiconstantinou et al. (Reference Hadjiconstantinou, Garcia and Alder2000), we set $\alpha =0.6$ (for all temperatures). Furthermore, we calculate the chemical potential through

(5.10)\begin{equation} \frac{\mu}{T} = \frac{\mu_{crit.}}{T_{crit.}} - \int_{p_v/T}^{p_v/T} \vartheta \, \textrm{d} (p_v/T), \end{equation}

where

(5.11)\begin{gather} \frac{\mu_{crit.}}{T_{crit.}} ={-} k_b (2\xi_{crit.} + \ln(\vartheta _{crit.}/\vartheta _{q,\textit{crit.}})), \end{gather}

(5.12)\begin{gather}\xi_c = \frac{\vartheta _{crit.} p_{crit.}}{k_b T_{crit.}}, \end{gather}

(5.13)\begin{gather}\textrm{and} \quad \vartheta _{q} = \frac{\hbar^3}{(2{\rm \pi} m k_b T)^{3/2}}. \end{gather}

Note that $\hbar$ is the Planck constant (see Elsner (Reference Elsner1989) for details). Figure 7(d) shows reasonable agreement between all three considered methods, i.e. analytical Elsner's approach, the DFP-SP model and the ESMC-Vlasov simulation results.

5.3. Evaporation to vacuum

Since streaming of the stochastic particle in the DFP model is consistent with the Enskog equation only on the moment level, the measurement of evaporation rate based on particle trajectory as suggested in Yasuoka, Matsumoto & Kataoka (Reference Yasuoka, Matsumoto and Kataoka1994), Kon et al. (Reference Kon, Kobayashi and Watanabe2014), Kobayashi et al. (Reference Kobayashi, Hori, Kon, Sasaki and Watanabe2016) would not lead to an accurate estimation of the desired property. Alternatively, the evaporation rates can be approximated via moments as a liquid slab evaporates into the vacuum; see Frezzotti et al. (Reference Frezzotti, Gibelli and Lorenzani2005). The outflow mass flux $\dot {m}_{{out}}$ is indeed a function of the evaporation rate $\sigma _e$,

(5.14)\begin{equation} \dot{m}_{{out}} =A \sigma_e n_g(T_L) \sqrt{\frac{k_b T_L}{2 {\rm \pi}m}}, \end{equation}

where $A$ indicates the area of cross-section normal to the interphase, $n_g(T_L)$ is the number density of the gas computed at the temperature of the liquid phase $T_L$. Hence, validation for the evaporation rate translates into the consistency of temperature and density profiles, as well as the outgoing mass flux with respect to the benchmark solution. Except for the boundary conditions, we consider similar simulation settings as the test case with the evaporation of liquid at equilibrium, i.e. § 5.2. Particles are reinitialized in the liquid slab as they leave the domain with the new position sampled randomly from a uniform distribution on the interval $x_2\in [L/2-L_{{therm}}/2,L/2+L_{{therm}}/2]$. For the reinitialized particle, the velocity component in the relevant flow direction, i.e. $V_2$, is sampled again from an equilibrium distribution with mean and temperature of the local computational cell, as the other components of the velocity remain unchanged.

In comparison to ESMC-Vlasov, we obtain accurate predictions of temperature, density, velocity and the outflow mass flux, as shown in figure 9. Similar test cases were performed at other temperatures and again good agreements between the DFP-SP model and the benchmark were observed, indicating consistent evaporation rates.

Figure 9.

Simulation results for evaporation of liquid argon into the vacuum obtained from ESMC-Vlasov and DFP-SP models. (a) Normalized number density. (b) Temperature. (c) Mean velocity. (d) Outflow mass flux.