On the entropic property of
the Ellipsoidal Statistical model
with the Prandtl number below 2/3

In the frame of the mutually orthogonal eigenvectors of $\Theta$, the $\mathcal{T}$ as well as $\Theta$ is diagonalized:

where $\mathrm{Id}$ is the $3\times 3$ identity matrix. Since the three diagonal components multiplied by $\lambda$ are eigenvalues of $\mathcal{T}$, the condition for $\mathcal{T}$ to be positive definite is simply reduced to

in the present parameter range $-1\leq\nu<-1/2$. Note that, since $0\leq\epsilon<1$ by construction, the constraint (2) is effective only when the right-hand side of the last inequality is less than unity, namely

is fulfilled. The above $\theta_{\mathcal{T}}$ varies monotonically from $1/2$ to $1$ when $\nu$ varies from $-1$ to $-1/2$.

After straightforward manipulations, the difference $D$ of $\det\mathcal{T}$ with $\det\Theta$ is written as

Here, we need to care the restriction in (ii) of the last paragraph of Sec. 2 that $D=0$ holds only when $\Theta$ is a scalar multiple of the identity matrix, which is identical to $\epsilon=0$ in the present notation. Hence the condition for $D\geq 0$ with this restriction is reduced to that the inside of the curly brackets is positive, i.e.,

for $\epsilon\neq 0$ in the range $-1\leq\nu<-1/2$. The left-hand side is non-negative because $0\leq\theta\leq 1$, while the sign of the right-hand side depends on that of the factor $F(\theta;\nu)\equiv[2\nu+1+(1-\nu)\theta][(2\nu+1)\theta+1-\nu]$. Let the smaller and larger zeros of $F$ with fixed $\nu$ be $\theta_{a}$ and $\theta_{b}$, i.e., $\theta_{a}=-(2\nu+1)/(1-\nu)$ and $\theta_{b}=-(1-\nu)/(2\nu+1)$. Then, it is readily seen that $0<\theta_{a}<1<\theta_{b}$ and that $F\geq 0$ for $\theta_{a}\leq\theta\leq\theta_{b}$ and $F<0$ elsewhere. Hence, we have the following condition:

We shall merge the conditions (2) and (6) by taking account of the original range of $\epsilon$, i.e., $0\leq\epsilon<1$.

Let us start with the observation that $\theta_{\mathcal{T}}-\theta_{a}=3(\nu+1)/(1-\nu)\geq 0$, so that $\theta_{a}\leq\theta_{\mathcal{T}}<1$ for $-1\leq\nu<-1/2$. As described before, the condition (2) becomes effective only when $\theta>\theta_{\mathcal{T}}$. Our first concern is the size relation between $3/[(1-\nu)(1+\theta)]$ in (2) and $\mathcal{F}(\theta;\nu)$ in (6) in the range $\theta_{\mathcal{T}}<\theta\leq 1$. Then, checking the size relation is reduced to checking the sign of the following quantity

because all of $(1+\theta)$, $(1-\nu)$, and $F$ are positive. Here, $\Delta\gtrless 0$ implies $3/[(1-\nu)(1+\theta)]\gtrless\mathcal{F}$. The simple substitution of $F$ shows, after some manipulations, that

Remind that $\theta_{\mathcal{T}}<1$ for $-1\leq\nu<-1/2$. Thus $\Delta>0$, i.e., $\mathcal{F}<3/[(1-\nu)(1+\theta)]$ for $\theta_{\mathcal{T}}<\theta\leq 1$. Our next concern is the size relation between $\mathcal{F}$ and unity in the range $\theta_{a}<\theta\leq\theta_{\mathcal{T}}$, because $3/[(1-\nu)(1+\theta)]$ is over or equal to unity in this range. Then, because $\Delta\leq 0$ for $\theta_{a}<\theta\leq\theta_{\mathcal{T}}$, we have $1\leq 3/[(1-\nu)(1+\theta)]\leq\mathcal{F}$, so that there is no restriction on $\epsilon$ in the concerned range of $\theta$.

To summarize, we have arrived at the following statement of criterion.

Let $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ be positive real eigenvalues of the symmetric positive definite matrix $\Theta$ defined by (1d). Let $\epsilon=1-\lambda_{1}/\lambda_{3}$ and $\theta=(1-\lambda_{2}/\lambda_{3})/(1-\lambda_{1}/\lambda_{3})$, where $0\leq\epsilon<1$ and $0\leq\theta\leq 1$ by construction. Then, for $-1\leq\nu<-1/2$, which corresponds to $1/2\leq\mathrm{Pr}<2/3$, the following criterion must be fulfilled in order for the H theorem to remain valid:

where $\mathcal{F}$ is the one defined by (6).

Consider a normal shock wave standing at rest against a uniform flow of a monatomic gas. Let $X_{1}$ be the spatial coordinate along the flow direction whose origin is located at the center of the shock.²²2The center of the shock is defined by the position where the density $\rho$ takes the average of the densities at the upstream and downstream infinities. Let $\rho_{-}$, $T_{-}$, and $(u_{-},0,0)$ be the density, temperature, and flow velocity of the uniform equilibrium state at upstream infinity, while $\rho_{+}$, $T_{+}$, and $(u_{+},0,0)$ be those at downstream infinity. Then, under the spatially uniform and the null flow-velocity assumption in the $X_{2}$- and $X_{3}$-directions, the problem is formulated as

where the upstream and the downstream quantities are related by the Rankine–Hugoniot conditions [12]:

Here $\gamma=5/3$ and $M_{-}(\geq 1)$ is the Mach number of the flow at upstream infinity, i.e., $M_{-}=u_{-}/\sqrt{\gamma RT_{-}}$.

The problem can be reduced to that of the marginal VDFs [5]. After the reduction, it is solved numerically by the standard iterative finite-difference method. The numerical accuracy is carefully checked enough to assure the main conclusions drawn from the present computations. Here we omit the details of the numerical method itself, grid resolution, etc. and simply present the obtained results.

It is readily seen that $\Theta$ is diagonal³³3Thanks to the assumption made just before (9a), the problem is rotationally invariant around the $X_{1}$-axis. The off-diagonal components are thus null. in the present coordinates system, and its diagonal components are nothing else than the eigenvalues $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ [or $(1-\epsilon)\lambda$, $(1-\epsilon\theta)\lambda$, and $\lambda$]. Firstly, the computations have been performed for $M_{-}=2$, $5$, $10$, 20, $40$, and $60$ under the setting $\mathrm{Pr}=2/3$. Figure 1 shows the structure of the shock wave in some typical cases. In all the cases, $\Theta_{11}$ is found to be the largest among the diagonal components. Because the problem is isotropic in the directions normal to $X_{1}$-direction, $\Theta_{22}=\Theta_{33}$ and thus $\theta\equiv 1$. Then, the condition (8) is simplified as

or equivalently

using the relation $\nu=1-1/\mathrm{Pr}$. Note that $\mathcal{S}(\mathrm{Pr})$ monotonically increases from $0$ to $1$ as $\mathrm{Pr}$ increases from $1/2$ to $2/3$. The obtained profiles of $\epsilon$, which is given by $1-\Theta_{33}/\Theta_{11}$ here, inside the shock are also shown by dash-dotted lines in Figure 1. With the aid of (12), the lower bound of the acceptable Prandtl number for each Mach number is suggested from the figure, the result of which is shown in Table 1.

With these results in mind, further computations have been carried out for different Prandtl numbers as well. The maximum values of $\epsilon$ in individual computations are plotted against $\mathrm{Pr}$ for each Mach number in Figure 2. It is seen that they depend only weakly on $\mathrm{Pr}$ for every fixed value of $M_{-}$. The intersections of solid lines and the curve $\epsilon=\mathcal{S}(\mathrm{Pr})$ are indicated by the symbol $\Box$. Their abscissas indicate the actual lower bounds of acceptable Prandtl number for the individual $M_{-}$. They are all slightly smaller than the values of $\mathrm{Pr}_{*}$ in Table 1 for those $M_{-}$. It is seen that in the present test the condition (8) [or (12)] is tolerant; even for $M_{-}=60$, $(3/2)\mathrm{Pr}$ can be down to $\simeq 0.83$, i.e., a reduction of $\mathrm{Pr}$ by around $17\%$ from $2/3$ is acceptable.

Incidentally, even for $\mathrm{Pr}=1/2$ (or $\nu=-1$), where the condition (12) is violated, stable numerical computations have been carried out. It would be partially explained by that the condition (2) for $\mathcal{T}$ to be positive definite is reduced to $0\leq\epsilon<3/4$ and is never violated in those computations (see Figure 2).

Consider a monatomic rarefied gas in the gap between two parallel plates which are kept at a common uniform temperature $T_{0}$. The confined gas is in the equilibrium state at rest with the temperature $T_{0}$ and the density $\rho_{0}$. Set the $X_{1}$-direction to be normal to the plates, with its origin at the middle of the gap. Let the gap width be $L$. At the instance $t=0$, the two plates suddenly start sliding with a constant speed $V$ in the directions opposite to each other: the plate at $X_{1}=L/2$ moves in the $X_{2}$-direction, while that at $X_{1}=-L/2$ in the $-X_{2}$-direction. Assuming the diffuse reflection boundary condition on the plates, the problem is formulated as

This is the time-dependent Couette flow problem.

As the previous test case, the problem can be reduced to that of the marginal VDFs [5]. After the reduction, it is solved numerically by the standard finite-difference method implicit in time. The numerical accuracy is carefully checked enough to assure the main conclusions drawn from the present computations. We again omit all the details of the numerical computations. In the present test, $\Theta$ is not diagonal in the specified Cartesian coordinates, and the three eigenvalues are given by $\Theta_{33}$ and $\lambda_{\pm}\equiv[\Theta_{11}+\Theta_{22}\pm\sqrt{(\Theta_{11}-\Theta_{22})^{2}+4\Theta_{12}^{2}}]/2$.

The computations have been carried out for various values of $V$, the reference Knudsen number, and the Prandtl number. In the present test, $\epsilon$ is found to be always $1-\lambda_{-}/\lambda_{+}$ and to take its maximum value on the plates at the initial stage of the time evolution. Consequently, that value is found to depend only on the plate speed, neither on the Knudsen number nor the Prandtl number, because at the initial stage the collision integral is not effective. These observations imply that the most severe situation is identical to the initial stage of the Rayleigh problem [2, 13] and is represented by the combination of the initial data for incident molecules and the boundary condition for reflected molecules on the plates. That is, the most severe test can be conducted by the following VDF:

Note that this is identical to the steady free molecular solution of the Couette flow with one plate sliding with speed $V$ and the other plate at rest. The above reduction of test, which we call a sliding plate test, allows us to assess the condition of our interest in detail.

For the above VDF (14), $\epsilon$ and $\theta$ are readily obtained as

where $\hat{V}=V/\sqrt{2RT_{0}}$. Note that (a) $\theta>1/2$ from the above, while $0<\theta_{a}\leq 1/2$ in the range $-1\leq\nu<-1/2$ or $1/2\leq\mathrm{Pr}<2/3$. It is also found that $0\leq\epsilon_{S}(\hat{V})<1$. Hence, plugging (15) into (8) yields the following criterion for the H theorem to remain valid in the present test:

where the relation $\nu=1-1/\mathrm{Pr}$ has been used to convert $\nu$ to $\mathrm{Pr}$. The function $\mathcal{F}_{D}$ is plotted by solid lines for various values of $\mathrm{Pr}$ in Figure 3, together with $\epsilon_{S}$ by a dashed line. It is seen from the figure that, by the 20 % (resp. 16 %) reduction of Prandtl number from $2/3$, the condition (16) is satisfied for the sliding speed $\hat{V}\lesssim 1.03$ (resp. 1.59).

Incidentally, in the time-dependent Couette flow test, even when the condition (8) is violated, numerical computations have been carried out stably and consistently, as far as the condition (2) is satisfied. This observation agrees with the shock wave test for $\mathrm{Pr}=1/2$.

Motivated by the sliding plate test, we consider another severe test by the following VDF:

where

The above VDF mimics the gas on the plate at the initial instance of a sudden plate motion in the normal direction with a constant speed $|U|$, by supposing that, before starting the plate motion, the gas is in the equilibrium state at rest with density $\rho_{0}$ and temperature $T_{0}$. Here, $U>0$ implies the pushing plate while $U<0$ implies the pulling plate.

In this case, $\Theta$ is diagonal with $\Theta_{22}=\Theta_{33}$ (thus $\theta=1$) as in the shock wave test, and the diagonal components are given by

where $\hat{U}=U/\sqrt{2RT_{0}}$ and $\hat{\rho}(\hat{U})=(1/2)[\hat{\sigma}(\hat{U})+1+\mathrm{erf}(\hat{U})]$ is the gas density divided by $\rho_{0}$. Note that $\hat{\sigma}>0$ because of its definition. Accordingly, $\Theta_{11}\gtrless RT_{0}$ for $\hat{U}\gtrless 0$ and $\epsilon$ is written as

in the present test. Since $\theta=1$, the condition (8) takes the same form as (12). Substitution of (19) into (12) eventually leads to the following condition:

where

Figure 4 shows the functions $\epsilon_{P}$ and $\hat{\rho}$ respectively in (a) and (b), where the values of $\mathcal{S}(\mathrm{Pr})$ are also shown for reference in (a). The present test assesses the robustness of the desired entropic property for small Prandtl numbers against the anisotropy caused by strong compression (or expansion). In this sense, the type of test is similar to that in the shock wave case. Nevertheless, the observed robustness looks quite different. One possible explanation is the compression rate measured by the density. In the shock wave test, the ratio of density $\rho_{+}/\rho_{-}$ between the upstream and the downstream infinity is at most four even in the limit $M_{-}\to\infty$, see (10). But here, the ratio of the density to its initial value $\hat{\rho}$ or its inverse is beyond this value already around $\hat{U}=1.70$ or $-0.73$. The present test is thus much harder than the shock wave test.