Quantifying Compressibility and Slip in Multiparticle Collision (MPC) Flow Through a Local Constriction

Abstract

The flow of a compressible fluid with slip through a cylinder with an asymmetric local constriction has been considered both numerically, as well as analytically. For the numerical work, a particle-based method whose dynamics is governed by the multiparticle collision (MPC) rule has been used together with a generalized boundary condition that allows for slip at the wall. Since it is well known that an MPC system corresponds to an ideal gas and behaves like a compressible, viscous flow on average, an approximate analytical solution has been derived from the compressible Navier–Stokes equations of motion coupled to an ideal gas equation of state using the Karman–Pohlhausen method. The constriction is assumed to have a polynomial form, and the location of maximum constriction is varied throughout the constricted portion of the cylinder. Results for centerline densities and centerline velocities have been compared for various Reynolds numbers, Mach numbers, wall slip values and flow geometries.

1. Introduction

Flows through microchannels and microtubes have become recent areas of interest due to new developments in the fabrication technology of microfluidic devices. Examples of applications include micro-gas turbine generators and bio-analytical devices. In order to implement flow scontrol measures or to optimize the design of bio-analytical devices, for example, a proper understanding of the flow through the device has to be developed. On the other hand, in gas microflows, compressibility effects can be important, and wall slip can be measurable, requiring incorporation of these in any numerical or analytical studies in this field. Particle-based methods, such as multiparticle collision dynamics (MPCD), are a means to simulate flows of a Newtonian, compressible, ideal gas, and slip effects can be incorporated very easily. Additionally, a constricted geometry is an ideal flow domain where compressibility effects can be important, for which an analytical solution is feasible. Our goal in this paper is to develop a better understanding, both theoretically and numerically, of the effects of compressibility and wall slip in a flow through a local constriction.

Flows through constrictions are popular in blood flow studies, and the analytical method used in this paper is an extension of the pioneering analysis carried out in [1–3]. The method used is called the Karman–Pohlhausen method, which essentially leads to the determination of the axial velocity profile. In [1–3], the fluid is considered to be Newtonian and incompressible, and the no-slip assumption is used, as would be common for blood flow applications. A more accurate pressure distribution was later developed for the same flow problem and presented in [4]. The same method was also used in [5], where the flow of an incompressible couple-stress fluid through a constriction was developed. In [6], a modified Karman–Pohlhausen method was proposed, and a general (2M)-degree polynomial was used for the flow field rather than a fourth degree polynomial, as per the original Karman–Pohlhausen method. In [7], slip was incorporated for incompressible, Newtonian flow through a local constriction. Weakly compressible flow with slip was later considered by [8,9], who also allowed for a flow geometry that is not necessarily symmetric about the location of maximum constriction. The results presented here are extensions of the results given in [8,9], giving more accurate expressions for the axial velocity profile.

Numerical works for flow through constrictions are two-fold. Discretization of the Navier–Stokes equations of motion for steady flow through stenoses was carried out by a number of authors for a Newtonian fluid [10–16]. Non-Newtonian models were considered numerically in [17–19] to name a few. All but [19] used the no-slip boundary condition. All of these works are for incompressible flows as they are applied to blood flow studiess. Particle-based numerical methods, such as the Lattice-Boltzmann method [20], dissipative particle dynamics (DPD) [21,22] and multiparticle collision (MPC) dynamics [8,9,16], have more recently led to numerical solutions for flow through a local constriction. The Lattice-Boltzmann method has also recently been used for blood flow studies in complex flow geometries for realistic cardiovascular flow domains [23–25]. The method has been reviewed recently in [26], and its use for complex flows has been reviewed in [27]. Except for [8,9], the results are numerical. Since compressible flows through constrictions can exhibit significant compressibility effects and since particle-based methods have compressibility built-in, such methods are ideal numerical means for simulating compressible flow through local constrictions.

Additional particle-based methods applied to blood flow studies in microvessels, for which deformable particles are modeled separately from the fluid in which they are suspended, include simulations with MPC [28,29] and DPD [30–33]. In [32], a Y-shaped bifurcation is considered, and [29] consider a complex flow domain. The simulations in this paper differ from these references in that the MPC fluid in this paper has point particles that neither deform nor aggregate, and there is only one type of particle in the system.

In this paper, the Karman–Pohlhausen method is used to develop the axial velocity distribution for steady, Newtonian flow through a stenosed vessel, allowing for slip at the wall, as well as compressibility. The analysis is a natural extension of [1] and an improvement to the results given in [8,9]. The flow geometry considered is axisymmetric, but asymmetric about the location of maximum constriction. Effects of compressibility, slip and flow geometry are assessed. Numerical results for flow through the same geometry using multiparticle collision (MPC) dynamics are also obtained and compared to the analytical solution.

2. Multiparticle Collision Dynamics

The particle system contains N identical point particles of unit mass that are distributed uniformly over cells on a regular three-dimensional lattice. Each cell, ξ, contains n particles on average. At discrete time intervals, Δt, the continuous positions, r_i, and velocities, v_i (i = 1,... N), are updated according to the multiparticle collision (MPC) dynamics originally developed in [34]. So as to ensure Galilean invariance, a random grid shift is implemented prior to each collision step as first introduced in [35]. The idealized collisions of the MPC algorithm then update the velocity of particle i according to:

$$\mathbf{v}\to {\mathbf{V}}_{\xi }+{\widehat{\omega }}_{\xi }\left({\mathbf{v}}_i-{\mathbf{V}}_{\xi }\right)$$

(1)

where ω̂_ξ is a stochastic rotation matrix that rotates the velocities by either +π/2 or −π/2 about a randomly chosen axis that varies from cell to cell and in time, and V_ξ is the average velocity of all particles in cell ξ in the pre-collision state [34].

Next, a constant external force accelerates the post-collision velocity of particle i in the z-direction according to:

$$v_z^i\to v_z^i+g\Delta t$$

(2)

where $v_z^i$ is the z-component of the velocity of particle i and g is the acceleration value.

To simulate isothermal flow conditions, a thermostat is applied to the system, so as to remove the energy that the external force pumps into the system. The velocity of each particle is rescaled according to a profile-unbiased Galilean invariant thermostat first introduced by [36], the details of which can be found in [8,9,16].

Finally, free-streaming of the particles updates the positions according to:

$${\mathbf{r}}_i\to {\mathbf{r}}_i+{\mathbf{v}}_i\Delta t$$

(3)

where the velocity here is the velocity after the collision, acceleration and thermostatting steps have taken place.

2.1. Boundary Conditions

Periodic boundary conditions are applied in the z-direction, and collisions with the cylinder walls follow the generalized boundary condition [8,9,16,37,38]:

$${\mathbf{v}}_n\to -{\mathbf{v}}_n$$

(4)

$${\mathbf{v}}_t\to \left(2\lambda -1\right){\mathbf{v}}_t$$

(5)

which is capable of incorporating macroscopic slip by means of changing the value of λ ∈ [0, 1]. No-slip flow is obtained with the λ = 0 bounce-back rule, while elastic collisions (λ = 1) would result in uniform flow through the pipe. For our simulations, we use λ ∈ [0, 0.5].

In order to compare the particle-based method with the analytical results, the particle-system is subjected to a cumulative averaging procedure as outlined in [16], where it was found that the averaging method is ideal for determining the macroscopic velocity profile for MPC flows.

Theoretical expressions for the viscosity coefficient of an MPC flow have been developed, and it has been shown that for our choice in ω̂:

$$\mu ={\mu }_{\mathit{kin}}+{\mu }_{\mathit{coll}}$$

(6)

where:

$${\mu }_{\mathit{kin}}=\left(\frac{n{k}_{B}T}{m{(\Delta x)}^3}\right)\Delta t\left[\frac{5n}{6(n-1+e^{-n})}-\frac{1}{2}\right]$$

(7)

$${\mu }_{\mathit{coll}}=\frac{m}{18\Delta x\Delta t}\left(n-1+e^{-n}\right)$$

(8)

and k_B is the Boltzmann constant, T the system temperature, Δx the length of a cubic cell in the lattice and n the average number of particles in a cell [34,35,39–43].

3. Theoretical Analysis

The governing equations of motion for a compressible, isothermal, viscous flow of an ideal gas are given by:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\text{conservation}\hspace{0.17em}\text{of}\hspace{0.17em}\text{mass})$$

(9)

$$\begin{array}{l}\begin{array}{l}\rho \frac{D}{Dt}\mathbf{u}=-\nabla P+\rho \mathbf{f}+\mu {\nabla }^2\mathbf{u}+\left(\kappa -\frac{2}{3}\mu \right)\nabla \left(\nabla \cdot \mathbf{u}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\text{conservation}\hspace{0.17em}\text{of}\hspace{0.17em}\text{momentum})\end{array}\hfill \end{array}$$

(10)

$$P=\frac{k_{B}T}{m}\rho ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\text{equation}\hspace{0.17em}\text{of}\hspace{0.17em}\text{state})$$

(11)

where ρ is the density, t is time, D/Dt = ∂/∂t + u · ∇ is the material derivative, u is the velocity vector, P is the pressure, f corresponds to an external force, μ is the viscosity, κ is the bulk viscosity, m is the mass of the fluid particle, k_B is the Boltzmann constant and T is the constant fluid temperature.

Assuming steady-state and axisymmetry, the velocity vector in cylindrical coordinates is assumed to have the form:

$$\mathbf{u}=\left(u_r,u_{\theta },u_z\right)=\left(u\left(r,z\right),0,w\left(r,z\right)\right)$$

(12)

together with ρ = ρ(r, z). Under the Stokes assumption (κ = 0), the governing equations, with an external force in the form f = (f_r, f_θ, f_z) = (0, 0, ρg) become:

$$\frac{\partial }{\partial r}(\rho u)+\frac{\partial }{\partial z}(\rho w)+\frac{\rho u}{r}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\text{mass})$$

(13)

$$\begin{array}{l}\rho \left(u\frac{\partial u}{\partial r}+w\frac{\partial u}{\partial z}\right)=-\frac{\partial P}{\partial r}+\mu \left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{{\partial }^{2}u}{\partial z^2}-\frac{u}{r^2}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{\mu }{3}\frac{\partial }{\partial r}(\nabla \cdot \mathbf{v}),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(r-\text{momentum})\end{array}$$

(14)

$$\begin{array}{l}\rho \left(u\frac{\partial w}{\partial r}+w\frac{\partial w}{\partial z}\right)=\rho g-\frac{\partial P}{\partial z}+\mu \left(\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial w}{\partial r}+\frac{{\partial }^{2}w}{\partial z^2}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{\mu }{3}\frac{\partial }{\partial z}(\nabla \cdot \mathbf{v}),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(z-\text{momentum})\end{array}$$

(15)

$$P(r,z)=\frac{k_{B}T}{m}\rho (r,z),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(\text{equation}\hspace{0.17em}\text{of}\hspace{0.17em}\text{state})$$

(16)

where:

$$\nabla \cdot \mathbf{v}=\frac{u}{r}+\frac{\partial u}{\partial r}+\frac{\partial w}{\partial z}$$

(17)

and the θ-momentum equation is identically satisfied.

As per [1], for a mild stenosis geometry, the r-momentum Equation (14) can be approximated as $\frac{\partial P}{\partial r}=0$, in which case, Equation (16) implies ρ = ρ(z), which can be used in Equation (13) to give:

$$\frac{u}{r}+\frac{\partial u}{\partial r}=-\frac{1}{\rho }\frac{\partial }{\partial z}(\rho w)$$

(18)

Using this in the last term of Equation (15), together with the assumption that $u\frac{\partial w}{\partial r}\ll w\frac{\partial w}{\partial z}$ allows us to write the system for determining w(r, z) and P (z) as:

$$\rho w\frac{\partial w}{\partial z}=\rho g-\frac{dP}{dz}+\mu \left(\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial w}{\partial r}+\frac{4}{3}\frac{{\partial }^{2}w}{\partial z^2}\right)-\frac{\mu }{3}\frac{\partial }{\partial z}\left(\frac{1}{\rho }\frac{\partial }{\partial z}(\rho w)\right)$$

(19)

$$P(z)=\frac{k_{B}T}{m}\rho (z)$$

(20)

Following [1], we now assume that the radial dependence of the axial velocity, w, is a fourth-order polynomial in the form:

$$\frac{w}{W}=A\eta +B{\eta }^2+C{\eta }^3+D{\eta }^4+E$$

(21)

where $\eta =\frac{R-r}{R}$, and W = W (z) is the as yet unknown centerline velocity. Constants A to E are determined by imposing:

(i)

$w=\frac{w_s}{\sqrt{1+{R^{\prime }}^2}}$ at r = R (slip boundary condition),

(ii)

$\frac{\partial w}{\partial r}=0$ at r = 0 (axisymmetric flow),

(iii)

w = W at r = 0 (by definition of centerline velocity W),

(iv)

$\frac{{\partial }^{2}w}{\partial r^2}=-\frac{2\left(W-w_s\right)}{R^2}$ at r = 0 (nearly parabolic flow with slip),

(v)

$\frac{dP}{dz}\approx \rho g+\mu \left(\frac{{\partial }^{2}w}{\partial r^2}+\frac{1}{r}\frac{\partial w}{\partial r}\right)$ at r = R (using (19)).

Condition (i) follows from solving u · n = 0 (the vanishing normal component of velocity) and u · t = w_s (the tangential component of velocity is w_s) for w, while (iv) comes from the velocity profile:

$$w^{\mathit{poi}}(r)=(W-w_s)\hspace{0.17em}\hspace{0.17em}\left[1-\frac{r^2}{R^2}\right]+w_s$$

(22)

which is Poiseuille flow in an unconstricted tube with slip, w_s, at the wall (r = R) and W is centerline velocity.

Imposing (i)–(v) and solving for the unknown constants gives:

$$A=\frac{1}{7}\left(-\lambda +10-12E+T+2\frac{w_s}{W}\right)$$

(23)

$$B=\frac{1}{7}\left(3\lambda +5-6E-3T+\frac{w_s}{W}\right)$$

(24)

$$C=\frac{1}{7}\left(-3\lambda -12+20E+3T-8\frac{w_s}{W}\right)$$

(25)

$$D=\frac{1}{7}\left(\lambda +4-9E-T+5\frac{w_s}{W}\right)$$

(26)

$$E=\frac{w_s}{W\sqrt{1+{R^{\prime }}^2}}$$

(27)

where:

$$\lambda =\frac{R^2}{\mu W}\frac{dP}{dz}$$

(28)

and:

$$T=\frac{\rho g{R}^2}{\mu W}$$

(29)

By definition, the flow rate is given by:

$$Q=\pi \rho R^2\overline{W}=\underset{0}{\overset{R}{\int }}2\pi \rho wr\hspace{0.17em}dr$$

(30)

Substituting Equation (21) for w, using Equations (23)–(27) and solving for W in terms of W̄, gives the relationship:

$$W=\frac{210}{97}\overline{W}+\frac{2}{97}\frac{R^2}{\mu }\frac{dP}{dz}-\frac{2}{97}\frac{R^2\rho g}{\mu }-\frac{11}{97}{w}_s-\frac{102}{97}\frac{w_s}{\sqrt{1+{R^{\prime }}^2}}$$

(31)

The details pertaining to the next step involving the derivation of the equation for $\frac{dP}{dz}$ are outlined in Appendix A. The result is:

$$\begin{array}{l}\frac{R^2}{\mu \overline{W}}\frac{dP}{dz}\left(1-\frac{388}{225}{\mathit{Ma}}^2+\frac{97}{225}\frac{w_s}{\overline{W}}{\mathit{Ma}}^2-\frac{194}{225}\frac{w_s}{\overline{W}}\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}\frac{{\mathit{Ma}}^2}{\mathit{Re}}\right)=\frac{388}{225}{R}^{\prime }\mathit{Re}+\frac{gR}{{\overline{W}}^2}\mathit{Re}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-8-\frac{8}{25}\frac{w_s}{\overline{W}}-\frac{97}{225}\frac{w_s^2}{{\overline{W}}^2}{R}^{\prime }\mathit{Re}+\frac{{\mathit{w}}_{\mathit{s}}}{\overline{\mathit{W}}}\frac{1}{\sqrt{1+{R^{\prime }}^2}}\left(\frac{208}{25}-\frac{194}{75}\frac{d}{dz}(R{R}^{\prime })\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{97}{75}\frac{R^{\prime }}{1+{R^{\prime }}^2}\frac{w_s^2}{{\overline{W}}^2}\mathit{Re}-\frac{388}{75}R{R}^{\prime }\frac{{\mathit{w}}_{\mathit{s}}}{\overline{\mathit{W}}}\frac{d}{dz}{(1+{R^{\prime }}^2)}^{-1/2}\end{array}$$

(32)

where we have defined the local Reynolds and Mach numbers as:

$$\mathit{Re}=\frac{\rho \overline{W}R}{\mu }$$

(33)

and:

$$\mathit{Ma}=\frac{\overline{W}}{\sqrt{\frac{k_{B}T}{m}}}$$

(34)

respectively.

Finally, substitution of Equations (31) and (32) in Equation (21) and subsequent simplification gives the axial velocity as:

$$\frac{w(\eta ,z)}{\overline{W}}=\frac{G\eta +H{\eta }^2+I{\eta }^3+J{\eta }^4}{\left(1-\frac{388}{225}{\mathit{Ma}}^2+\frac{97}{225}\frac{w_s}{\overline{W}}{\mathit{Ma}}^2-\frac{194}{225}\frac{w_s}{\overline{W}}\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}\frac{{\mathit{Ma}}^2}{\mathit{Re}}\right)}+K$$

(35)

where G, H, I, J and K are given in Appendix B.

Substituting η = 1 and simplifying gives the centerline velocity as:

$$\begin{array}{l}\left[\frac{w\left(\eta =1,z\right)}{\overline{W}}-\frac{w_s}{\overline{W}\sqrt{1+{R^{\prime }}^2}}\right]\left(1-\frac{388}{225}{\mathit{Ma}}^2+\frac{97}{225}\frac{w_s}{\overline{W}}{\mathit{Ma}}^2-\frac{194}{225}\frac{w_s}{\overline{W}}\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}\frac{{\mathit{Ma}}^2}{\mathit{Re}}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=2+\mathit{Re}\frac{\mathit{dR}}{\mathit{dz}}\left[\frac{8}{225}-\frac{2}{225}\frac{w_s^2}{{\overline{W}}^2}+\frac{6}{225}\frac{1}{1+{R^{\prime }}^2}\frac{w_s^2}{{\overline{W}}^2}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{75}\frac{w_s}{\overline{W}}\left[-9-\frac{141}{\sqrt{1+{R^{\prime }}^2}}-\frac{4}{\sqrt{1+{R^{\prime }}^2}}\left({R^{\prime }}^2+R{R}^{″}\right)+8\frac{R{R^{\prime }}^2{R}^{″}}{{(1+{R^{\prime }}^2)}^{3/2}}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\mathit{Ma}}^2[-\frac{56}{15}+\frac{2}{225}\frac{gR}{{\overline{W}}^2}\left(4\mathit{Re}-\frac{w_s}{\overline{W}}\mathit{Re}+2\frac{w_s}{\overline{W}}\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{225}\frac{w_s}{\overline{W}}\left(254-11\frac{w_s}{\overline{W}}+\frac{796}{\sqrt{1+{R^{\prime }}^2}}-\frac{199}{\sqrt{1+{R^{\prime }}^2}}\frac{w_s}{\overline{W}}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{2}{225}\frac{1}{\mathit{Re}}\frac{w_s}{\overline{W}}\left(-210\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}+11\frac{w_s}{\overline{W}}\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}+199\frac{w_s}{\overline{W}}\frac{R^{\prime }}{1+{R^{\prime }}^2}\right)\end{array}$$

(36)

Note that substituting Ma = 0 and $\frac{dR}{dz}=0$ leads to $\frac{W}{\overline{W}}=2-\frac{w_s}{\overline{W}}$, which agrees with Equation (A9) for w = w^poi, as it should, and that substitution of Ma = 0 and w_s = 0 for $\frac{dR}{dz}\ne 0$ in the above solution gives Forrester and Young’s [1] result for incompressible no-slip flow.

4. Equation for Density

In order to plot the velocity profile obtained in the previous section, the explicit solution for ρ(z) has to be found, since Re and Ma depend on ρ(z), due to their local nature. To achieve this, the ideal gas equation of state Equation (20) can be used to replace pressure terms with density in Equation (A2), while constant flow rate can be used to replace local Re and Ma numbers with upstream values and ρ(z) terms. Specifically, constant flow rate implies (see Equation (30)):

$$\overline{W}=\frac{{\overline{W}}_0{\rho }_0{R}_0^2}{\rho R^2}$$

(37)

where the zero subscript indicates constant upstream values in the unconstricted portion of the cylinder. Thus:

$$\begin{array}{l}\mathit{Re}=\frac{\rho \overline{W}R}{\mu }=\frac{\rho {\overline{W}}_0{\rho }_0{R}_0^{2}R}{\rho R^2\mu }\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=R{e}_0\frac{R_0}{R}\end{array}$$

(38)

and:

$$\begin{array}{l}\mathit{Ma}=\frac{\overline{W}}{\sqrt{\frac{k_{B}T}{m}}}=\frac{{\overline{W}}_0{\rho }_0{R}_0^2}{\rho R^2\sqrt{\frac{k_{B}T}{m}}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\mathit{Ma}}_0\frac{{\rho }_0}{\rho }{\left(\frac{R_0}{R}\right)}^2\end{array}$$

(39)

Lastly, the dimensionless slip velocity can be written as:

$$\frac{w_s}{\overline{W}}=\frac{w_s}{{\overline{W}}_0}{\left(\frac{R}{R_0}\right)}^2\frac{\rho }{{\rho }_0}$$

(40)

It follows that the pressure equation can be written in terms of ρ(z) using the equation of state, Equation (20), giving:

$$\begin{array}{l}\frac{-R_0^2}{\mu {\overline{W}}_0}{\left(\frac{R}{R_0}\right)}^4\frac{\rho }{{\rho }_0}\frac{k_{B}T}{m}\frac{d\rho }{dz}[1-\frac{388}{225}{\mathit{Ma}}_0^2{\left(\frac{{\rho }_0}{\rho }\right)}^2{\left(\frac{R_0}{R}\right)}^4+\frac{97}{225}\frac{w_s}{\overline{W}}{\mathit{Ma}}_0^2-\frac{{\rho }_0}{\rho }{\left(\frac{R_0}{R}\right)}^2\\ -\frac{194}{225}\frac{w_s}{{\overline{W}}_0}\frac{R^{\prime }}{\sqrt{1+{R^{\prime }}^2}}\frac{{\mathit{Ma}}_0^2}{\mathit{R}{\mathit{e}}_0}\frac{{\rho }_0}{\rho }\frac{R_0}{R}]=\frac{388}{225}{R}^{\prime }\mathit{Re}+\frac{\mathit{gR}}{{\overline{\mathit{W}}}^{\mathit{2}}}\mathit{Re}\\ -8-\frac{8}{25}\frac{w_s}{\overline{W}}-\frac{97}{225}\frac{w_s^2}{{\overline{W}}^2}{R}^{\prime }\mathit{Re}+\frac{w_s}{\overline{W}}\frac{1}{\sqrt{1+{R^{\prime }}^2}}\left(\frac{208}{25}-\frac{194}{75}\frac{d}{dz}(R{R}^{\prime })\right)\\ +\frac{97}{751}+\frac{R^{\prime }}{1+{R^{\prime }}^2}\frac{w_s^2}{{\overline{W}}^2}\mathit{Re}-\frac{388}{75}R{R}^{\prime }\frac{w_s}{\overline{W}}\frac{d}{dz}{(1+{R^{\prime }}^2)}^{-1/2}\end{array}$$

(41)

where Re and Ma must be written in terms of Re₀, Ma₀ and ρ as given by Equations (38)–(40).

5. Flow Geometry

In order to be able to consider an asymmetric stenosis, the radius is taken to have the idealized polynomial form:

$$R\left(z\right)=\{\begin{array}{ll}{R}_0,\hfill & \text{for}\hspace{0.17em}z\le z_1\hfill \\ a{z}^3+b{z}^2+cz+d,\hfill & \text{for}\hspace{0.17em}{z}_1\le z\le z_2\hfill \\ e{z}^3+f{z}^2+gz+h,\hfill & \text{for}\hspace{0.17em}{z}_2\le z\le z_3\hfill \\ R_0\hfill & \text{for}\hspace{0.17em}z\ge z_3\hfill \end{array}$$

(42)

where z₂ = z₁ + l₁ and z₃ = z₂ + l₂.

Imposing that R(z) be continuously differentiable and that R(z₂) = R₀ – δ and R′(z₂) = 0 requires:

$$a=\frac{2\delta }{l_1^3}$$

(43)

$$b=-\frac{3\delta \left(2z_1+l_1\right)}{l_1^3}$$

(44)

$$c=\frac{6\delta z_1\left(z_1+l_1\right)}{l_1^3}$$

(45)

$$d=-\frac{2\delta z_1^3+3\delta z_1^2{l}_1-R_0{l}_1^3}{l_1^3}$$

(46)

$$e=-\frac{2\delta }{l_2^3}$$

(47)

$$f=\frac{3\delta \left(2z_1+2l_1+l_2\right)}{l_2^3}$$

(48)

$$g=-\frac{6\delta \left(z_1+l_1\right)\left(z_1+l_1+l_2\right)}{l_2^3}$$

(49)

$$h=\frac{3\delta l_2\left(z_1^2+l_1^2\right)+6\delta z_1{l}_1\left(z_1+l_1+l_2\right)+2\delta \left(z_1^3+l_1^3\right)+\left(R_0-\delta \right)l_2^3}{l_2^3}$$

(50)

The resulting axisymmetric flow domain is shown in Figure 1. As can be seen from the figure, by construction, δ controls the severity of the constriction, while l₁ can be used to create the asymmetry about the z₂ location.

For all results that follow, R₀ = 10.5, z₁ = 600.5 and l₁ + l₂ = 30.

6. Numerical Results and Discussion

For all (dimensionless) MPC simulations that follow, there were approximately N = 8.5 million particles of unit mass m = 1 in the system, Δx = 1 = Δy =Δz; there were 1, 200 cells in the z direction and 25 cells in the x and y directions, respectively. The time step was taken to be Δt = 1 and k_BT = 1, together with n = 20. For the cumulative average, the averaging started after 5, 000 time steps and was performed for 35, 000 time steps thereafter. The initial system was set up with x and y velocities drawn from a Maxwellian velocity distribution, and z velocity drawn from the steady velocity profile of flow through a cylinder of fixed radius R₀.

A length of 1, 200 cells in the z-direction was chosen so as to ensure that periodic boundary conditions are valid. For this cylinder length, the velocity settled back to the expected parabolic profile in an unconstricted cylinder prior to reaching the exit for all constrictions considered here. In addition, since the velocity and density were found to be affected upstream in some simulations, starting the constriction at z = 600.5 ensured that there was a region upstream for which this effect was not present. Although some constrictions did not require a length of 1,200, this length was fixed for all simulations, so as to ensure that the most severe constriction with the highest Reynolds number would satisfy the periodic boundary condition.

The initial velocity distribution in the z direction was chosen, so as to reduce the simulation time. Test simulations (not reported here) were performed using a Maxwellian velocity distribution in all three directions as the initial state. The system maintained the Maxwellian velocity distribution in the x and y directions, and on average, the expected z velocity distribution that was later chosen as the initial state. In this way, the system reached equilibrium earlier, and the cumulative averaging could start after 5,000 time steps in all cases considered.

Simulations were done using serial code on an Intel Xeon X5482 3.2 GHz machine with 8 GB RAM. Typical run times were 3–4 days.

To obtain the required upstream values for ρ₀, the particle-based numerical results were averaged over the centerline density values for z ∈ [0, 100], and a best parabolic fit to the cross-section at z = 100.5 gave rise to the values for W̄₀ and w_s, as provided in

Table 1. Parameter values used in the analytical solution in Figure 2 for comparison with the particle-based method for compressible no-slip flow (λ = 0, w_s = 0).

g	ρ0	ρequil	$\frac{{\rho }_0-{\rho }^{\mathit{equil}}}{{\rho }_0}$	W̄0	Re0
0.005	20.025321	20.55929025	−0.0267	0.168938215102975	4.126
0.01	20.187479	20.73938113	−0.0273	0.338610045766591	8.277
0.02	20.9408427	21.57604487	−0.0303	0.683610724092841	16.770

and

Table 2. Parameter values used in the analytical solution in Figure 3 for comparison with the particle-based method for compressible flow with slip (λ = 0.5).

g	ρ0	ρequil	$\frac{{\rho }_0-{\rho }^{\mathit{equil}}}{{\rho }_0}$	ws	W̄0	Re0
0.005	12.3602273	12.70014982	−0.0275	0.0358641	0.164985484616928	3.784
0.01	12.4565002	12.57105274	−0.0092	0.0733492	0.331580364858540	7.616
0.02	12.8636085	12.53377480	+0.0256	0.145642	0.670483591369728	15.479

. These values were then used to determine the density from numerical integration of Equation (41).

It can be seen in Figure 4 that the bounce-back rule (MPC-BB, λ = 0) correctly leads to the expected zero velocity at the wall, while slip is clearly present in the MPC-LIT(λ = 0.5) case.

The differential equation for density was found to have a stable positive steady state, ρ^equil, that differed slightly from the ρ₀ determined from the MPC results. The values have been added to the Tables, as well as the relative errors from ρ₀. The density equation was solved numerically using the fourth-order Runge–Kutta scheme with Δz = 0.001 using MAPLE. Since the geometry is a piecewise defined function, the equation was solved one piece at a time, and instead of imposing ρ₀ as an initial condition at z = 0, ρ^equil was used. The differential equation was then solved on [0, z₁] with the value at z₁ becoming the initial condition for the differential equation on [z₁, z₂], and so on. In this way, the numerical solution was found for z ∈ [0, 1200]. Since the system has a steady state, the density settled back to the equilibrium value downstream of the constriction, thus ensuring that periodic boundary conditions are obtained in the analysis allowing comparison with the MPC results.

6.1. Compressible No-Slip Flow

In Figure 2a, a comparison of the theoretically-predicted centerline density arising from the numerical solution of Equation (41) is made with the particle-based MPC density results in the no-slip case. It can be seen that although there are some discrepancies between the predicted density curves and those obtained from the MPC simulations, both predict a density increase through the constriction, and the best agreement is found for the lowest Reynolds number considered (g = 0.005 curves). Worth noting in Table 1 is the increase in ρ₀ as Re₀ increases, which is consisted with the increase in ρ^equil.

Using the theoretically-predicted density curves in the centerline velocity expression (36) gives rise to the theoretically-predicted centerline velocity curves in Figure 2b. It can be seen that the theoretically-predicted centerline velocity agrees fairly well with the MPC result for g = 0.005,but as the Reynolds number increases, the agreement worsens. Worth noting is the appearance of a dip in the centerline velocity in both the theoretically-predicted and MPC results as a result of the constriction for the largest Reynolds number considered (g = 0.02).

6.2. Compressible Flow with Slip

For compressible flow with slip at the wall, relevant parameter values arising from the theoretical and numerical results are shown in Table 2. Theoretical scaled centerline densities and centerline velocities are compared to MPC results in Figure 3. It can be seen in (a) of the figure that there is some discrepancy between the theoretically predicted and MPC density results, but that the agreement is somewhat better than in the no-slip case. Likely due to the better agreement between the density curves, the scaled centerline velocities agree better, as well, and the dip for the largest Reynolds number (g = 0.02)is slightly overestimated by the theoretical predictions, contrary to the no-slip case.

Worth noting here is that, although the density curves seem to match better in the slip case, glancing at Table 2, ρ₀ is found to increase as the Reynolds number increases, while the reverse is predicted with ρ^equil.

6.3. Effect of the Severity of the Constriction

For the smallest Reynolds number considered (g = 0.005), the severity of the constriction is varied for both slip (λ = 0.5) and no-slip (λ = 0) flow. Corresponding parameter values are given in

Table 3. Parameter values used in the analytical solution in Figure 5 for comparison with particle-based method for compressible flow through constrictions of varying degrees.

δ	λ	ρ0	ρequil	ws	W̄0	Re0
0.5	0	20.025321	20.55929025	0	0.168938215102975	4.126
0.5	0.5	12.3602273	12.70014982	0.0358642	0.164985484616928	3.784
1.5	0	19.9307477	20.43003857	0	0.168324883949003	4.109
1.5	0.5	12.3319318	12.67373252	0.0379518	0.165379390657969	3.792
2	0	19.8767493	20.32769045	0	0.167506930369402	4.088
2	0.5	12.3103817	12.61661932	0.0356837	0.163995726106273	3.759

, and resulting scaled centerline velocity plots are shown in Figure 5. It can be seen that there is relatively good agreement between the theoretically-predicted curves and those from the MPC results for the mildest constriction (δ = 0.5) and that there is some discrepancy as the constriction becomes more severe. The appearance of a dip in the scaled centerline velocity for the more severe constrictions is captured in the slip case, while the decrease in scaled centerline velocity upstream of the constriction is found in the MPC no-slip results, but not in the theoretical predictions. On these same no-slip plots, the theoretical results predict a lower scaled centerline velocity in the post-constriction region, while MPC results do not show this feature.

6.4. Effect of Increasing Slip

Increasing the slip parameter, λ, and, thus, the wall slip, w_s, leads to Figure 6. Parameter values used in the analytical velocity profiles are provided in

Table 4. Parameter values used in the analytical solution in Figure 6 for comparison with the particle-based method for compressible flow through a constriction with δ = 0.5, g = 0.005, l₁ = 20 and variable slip parameter values.

λ	ρ0	ρequil	ws	W̄0	Re0
0	20.025321	20.55929025	0	0.168938215102975	4.126
0.2	14.9570845	14.90941499	0.00852200	0.155209851730363	3.663
0.4	12.6323301	12.66736550	0.0242417	0.155669931290677	3.583
0.5	12.3602273	12.70014982	0.0358642	0.164985484616928	3.784

. There is very good agreement between the analytical and the numerical results as the slip is varied, and the equilibrium density values from the theoretical predictions agree well with the centerline densities obtained in the MPC results.

6.5. Contour Plot Comparison

Figure 7 shows contour plots for the scaled centerline velocity for both the analytical and numerical particle-based method results for a constriction with δ = 2, g = 0.005, l₁ = 20 and λ = 0.5. For the analytical results, the values of the last row of Table 3 were used.

7. Discussion and Conclusions

An approximate analytical solution for the density, and for the axial velocity distribution, in an asymmetric constriction have been developed and compared to the numerical solution of a particle-based system governed by the multiparticle collision (MPC) dynamics. The solutions in all cases correspond to compressible flow with slip at the cylinder wall. Reynolds numbers varied from approximately four to 17.

Analysis of results revealed that increasing the Reynolds number in a fixed geometry leads to the appearance of a dip in the scaled centerline velocity in the entry region of the constriction, together with more pronounced flow acceleration following the location of maximum constriction. This is true with and without slip. In addition, as the Reynolds number increases, there is an increase in scaled centerline density, ρ₀, in all cases considered, except in the analytical results with slip that predict a decrease in centerline density (ρ^equil) instead. As the severity of the constriction increases, both slip and no-slip results show acceleration through the constriction, although the analytical and MPC results agree best for the mildest constriction (δ = 0.5) considered. Consistent with theory and MPC is the appearance of a dip in the scaled centerline velocity in the post-constriction region that is more pronounced as the severity of the constriction increases. This dip is, however, missing from the no-slip MPC results, which, instead, show a dip in the upstream section that is not captured in the theoretical predictions. Lastly, increasing slip has the effect of leading to faster flow through the constriction with the appearance of a dip in the post-constriction region that is consistent with the MPC results.

Since many key features compare well between the theoretically predicted results and those obtained by MPC, it is expected that improvements in the theory will lead to even better agreement in the constriction region and thereafter. In particular, an approximation was made for ${\int }_0^{R}r{w}^{2}dr$, which led to some errors in the pressure equation and all equations in the subsequent analysis. In Figure 8, plots of W and W_approx = 2W̄ − w_s can be found for the constrictions considered in Figures 2 and 3. It can be seen that relationship Equation (A9) is true for the smallest constriction considered and fails to hold for the higher Reynolds numbers, more so for the no-slip case in (a). This is likely a key reason as to why the agreement between MPC and theory is worse for larger Reynolds numbers. Furthermore, all quadratic (dP/dz)² and second-order d²P/dz² terms were dropped in the analysis, which likely led to some errors, as well. It would be interesting to explore whether or not keeping such terms in the analysis leads to significant improvements over what was found here, and this is currently under investigation. An additional source of discrepancy between the results could be the use of a thermostat in the MPC simulations that is applied uniformly, rather than locally, and whether or not using a local thermostat leads to better agreement is currently under investigation. A discussion on the use of thermostats in MPC simulations has been given in [43,44], and it would be interesting to see whether or not changing the thermostat in the simulations can lead to better agreement with the theoretical results.

In summary, an analytical solution for the flow of a compressible Newtonian fluid with slip at the wall was developed and found to compare fairly well to a numerical solution for a particle-based fluid governed by MPC in mild constrictions with low Reynolds numbers. Various Reynolds numbers, Mach numbers, wall slip values and flow geometries were considered in the analysis for asymmetric flow domains.