Two-dimensional off-lattice Boltzmann model for van der Waals fluids with variable temperature

TL;DR

This paper introduces a novel two-dimensional off-lattice Boltzmann model for van der Waals fluids with variable temperature, demonstrating high accuracy and stability in simulating phase separation and thermal flows.

Contribution

It presents a new off-lattice Boltzmann approach with fourth-order accuracy for thermal liquid-vapor systems, including variable temperature effects and phase separation.

Findings

Achieves at least fourth order convergence in shear and wave propagation.

Maintains small spurious velocities (<1%) even at high temperature ratios.

Preserves Galilean invariance up to second order.

Abstract

We develop a two-dimensional Lattice Boltzmann model for liquid-vapour systems with variable temperature. Our model is based on a single particle distribution function expanded with respect to the full-range Hermite polynomials. In order to ensure the recovery of the hydrodynamic equations for thermal flows, we use a fourth order expansion together with a set of momentum vectors with 25 elements whose Cartesian projections are the roots of the Hermite polynomial of order Q = 5. Since these vectors are off-lattice, a fifth-order projection scheme is used to evolve the corresponding set of distribution functions. A fourth order scheme employing a 49 point stencil is used to compute the gradient operators in the force term that ensures the liquid-vapour phase separation and diffuse reflection boundary conditions are used on the walls. We demonstrate at least fourth order convergence with respect to the lattice spacing in the contexts of shear and longitudinal wave propagation through the van der Waals fluid. For the planar interface, fourth order convergence can be seen at small enough lattice spacings, while the effect of the spurious velocity on the temperature profile is found to be smaller than 1.0%, even when T w ' 0.7 T c . We further validate our scheme by considering the Laplace pressure test. Galilean invariance is shown to be preserved up to second order with respect to the background velocity. We further investigate the liquid-vapour phase separation between two parallel walls kept at a constant temperature T w smaller than the critical temperature T c and discuss the main features of this process.

Tables

Table 1. Table 1: The roots ℛ q subscript ℛ 𝑞 {\cal{R}}_{q} of the full-range Hermite polynomial of order ℓ = 5 ℓ 5 \ell=5 and their associated weights 𝒲 q subscript 𝒲 𝑞 {\cal{W}}_{q} [ 11 ] .

$q$		${ℛ}_q$		${𝒲}_q$
1		$-\sqrt{5+\sqrt{10}}$		$\left(7-2\sqrt{10}\right)/60$
2		$-\sqrt{5-\sqrt{10}}$		$\left(7+2\sqrt{10}\right)/60$
3		$0$		$8/15$
4		$\sqrt{5-\sqrt{10}}$		$\left(7+2\sqrt{10}\right)/60$
5		$\sqrt{5+\sqrt{10}}$		$\left(7-2\sqrt{10}\right)/60$

Table 2. Table 2: The values of the weighting factors ω ¯ q subscript ¯ 𝜔 𝑞 \overline{\omega}_{q} ( 28 ) when one, two or all three of the indicators of smoothness σ q subscript 𝜎 𝑞 \sigma_{q} ( q = 1 , 2 , 3 𝑞 1 2 3 q=1,2,3 ) have vanishing values.

	${\overline{\omega }}_1$	${\overline{\omega }}_2$	${\overline{\omega }}_3$
${\varsigma }_1={\varsigma }_2={\varsigma }_3=0$	$0.1$	$0.6$	$0.3$
${\varsigma }_2={\varsigma }_3=0$	$0$	$2/3$	$1/3$
${\varsigma }_3={\varsigma }_1=0$	$1/4$	$0$	$3/4$
${\varsigma }_1={\varsigma }_2=0$	$1/7$	$6/7$	$0$
${\varsigma }_1=0$	$1$	$0$	$0$
${\varsigma }_2=0$	$0$	$1$	$0$
${\varsigma }_3=0$	$0$	$0$	$1$

Table 3. Table 3: Summary of the values for the surface tension γ pl subscript 𝛾 pl \gamma_{\rm pl} and γ Lap subscript 𝛾 Lap \gamma_{\rm Lap} obtained by evaluating numerically Eq. ( 51 ) for the planar interface profiles shown in Fig.  4 (b) and using the Laplace pressure test, as shown in Fig.  7 (c).

$\sigma$	${\gamma }_{\mathrm{pl}}$	${\gamma }_{\mathrm{Lap}}$
$5\times {10}^{-5}$	$3.386731\times {10}^{-3}$	$3.40689\times {10}^{-3}$
${10}^{-4}$	$4.883771\times {10}^{-3}$	$4.86615\times {10}^{-3}$

Table 4. Table 4: Exponents of δ ​ x = 1 / N x 𝛿 𝑥 1 subscript 𝑁 𝑥 \delta x=1/N_{x} for the relative error [ ( η app / η 0 ) − 1 ] delimited-[] subscript 𝜂 app subscript 𝜂 0 1 [(\eta_{\rm app}/\eta_{0})-1] and for the L 2 subscript 𝐿 2 L_{2} norm ( 59 ) for the ideal gas and vapour and liquid phases of the van der Waals fluid, for various longitudinal velocities, in the context of the damping of shear waves.

$V$	$\eta$ order			$L_2$ order
	Ideal gas	Vapour	Liquid	Ideal gas	Vapour	Liquid
$0.0$	$4.927334$	$4.927341$	$4.927331$	$4.907067$	$4.907078$	$4.907064$
$0.5$	$4.914673$	$4.914681$	$4.914674$	$4.892658$	$4.892671$	$4.892660$
$1.0$	$4.893496$	$4.893502$	$4.893497$	$4.864938$	$4.864947$	$4.864940$
$2.0$	$4.876409$	$4.876413$	$4.876410$	$4.831781$	$4.831786$	$4.831781$

Table 5. Table 5: Exponents of V 𝑉 V for the differences between the values measured at longitudinal velocity V 𝑉 V and at rest V = 0 𝑉 0 V=0 of the relative apparent shear viscosity [ ( η app ; V − η app ; V = 0 ) / η 0 ] delimited-[] subscript 𝜂 app V subscript 𝜂 app V 0 subscript 𝜂 0 [(\eta_{\rm app;V}-\eta_{\rm app;V=0})/\eta_{0}] ( η 𝜂 \eta order) and L 2 subscript 𝐿 2 L_{2} norm ( L 2 ; V − L 2 ; V = 0 ) subscript 𝐿 2 𝑉 subscript 𝐿 2 𝑉 0 (L_{2;V}-L_{2;V=0}) ( L 2 subscript 𝐿 2 L_{2} order), for the ideal gas and the vapour and liquid phases of the van der Waals fluid, in the context of the damping of shear waves.

$N_x$	Ideal gas		Vapour		Liquid
	$\eta$ order	$L_2$ order	$\eta$ order	$L_2$ order	$\eta$ order	$L_2$ order
$20$	$1.848778$	$1.848780$	$1.848780$	$1.845612$	$1.845614$	$1.845614$
$30$	$2.003266$	$2.003267$	$2.003280$	$2.003027$	$2.003027$	$2.003041$
$40$	$1.978353$	$1.978370$	$1.978403$	$1.978272$	$1.978290$	$1.978323$
$50$	$1.987273$	$1.987361$	$1.987420$	$1.987260$	$1.987352$	$1.987413$
$60$	$1.982981$	$1.983316$	$1.983342$	$1.983050$	$1.983397$	$1.983421$

Table 6. Table 6: Exponents in the context of the damping of a longitudinal wave of δ ​ x = 1 / 𝔑 x 𝛿 𝑥 1 subscript 𝔑 𝑥 \delta x=1/\mathfrak{N}_{x} for ( ν a ; app / ν a ; lin ) − 1 subscript 𝜈 𝑎 app subscript 𝜈 𝑎 lin 1 (\nu_{a;{\rm app}}/\nu_{a;{\rm lin}})-1 , ( c s ; app / c s ; lin ) − 1 subscript 𝑐 𝑠 app subscript 𝑐 𝑠 lin 1 (c_{s;{\rm app}}/c_{s;{\rm lin}})-1 and the L 2 subscript 𝐿 2 L_{2} norm, for various values of the background velocity V 𝑉 V .

$V$	${\nu }_a$ order			$L_2$ order			$c_s$ order
	Ideal gas	Vapour	Liquid	Ideal gas	Vapour	Liquid	Ideal gas	Vapour	Liquid
$0.0$	$4.576$	$4.579$	$4.311$	$4.480$	$4.452$	$4.235$	6.510	6.488	6.983
$0.5$	$4.565$	$4.571$	$4.294$	$4.469$	$4.439$	$4.221$	5.960	5.860	7.423
$1.0$	$4.533$	$4.540$	$4.258$	$4.433$	$4.407$	$4.178$	5.906	5.803	7.018
$2.0$	$4.554$	$4.533$	$4.481$	$4.095$	$4.004$	$3.781$	1.515	0.942	1.479

Table 7. Table 7: Exponents in the context of the damping of a longitudinal wave of V 𝑉 V for the differences ( ν a ; app ; V − ν a ; app ; 0 ) / ν a ; lin subscript 𝜈 𝑎 app 𝑉 subscript 𝜈 𝑎 app 0 subscript 𝜈 𝑎 lin (\nu_{a;{\rm app};V}-\nu_{a;{\rm app};0})/{\nu_{a;{\rm lin}}} ; ( L 2 ) V − ( L 2 ) 0 subscript subscript 𝐿 2 𝑉 subscript subscript 𝐿 2 0 (L_{2})_{V}-(L_{2})_{0} ; and ( c s ; app ; V − c s ; app ; 0 ) / c s ; lin subscript 𝑐 𝑠 app 𝑉 subscript 𝑐 𝑠 app 0 subscript 𝑐 𝑠 lin (c_{s;{\rm app};V}-c_{s;{\rm app};0})/c_{s;{\rm lin}} , where the quantities bearing the subscript V 𝑉 V are obtained by numerically fitting Eq. ( 79 ) to the data obtained when the background longitudinal velocity has value V 𝑉 V (the subscript 0 0 denotes the laboratory frame, in which V = 0 𝑉 0 V=0 ). The quantities bearing the subscript lin lin {\rm lin} refer to the solution of the linearised equations derived in Eq. ( 70 ).

${𝔑}_x$	${\nu }_a$ order			$L_2$ order			$c_s$ order
	Ideal gas	Vapour	Liquid	Ideal gas	Vapour	Liquid	Ideal gas	Vapour	Liquid
$20$	$1.951$	$1.955$	$1.944$	$1.962$	$1.970$	$1.928$	1.807	1.796	2.164
$30$	$1.902$	$1.902$	$1.903$	$1.874$	$0.133$	$1.884$	1.940	1.939	1.937
$40$	$1.889$	$1.884$	$1.867$	$1.822$	$1.982$	$1.891$	1.940	1.939	1.770
$50$	$1.767$	$1.326$	$2.076$	$1.840$	$1.965$	$1.897$	1.940	1.939	1.940
$60$	$1.992$	$1.973$	$1.961$	$1.865$	$1.961$	$1.900$	1.940	1.939	1.940

Equations

$${\tilde{n}}_R=\frac{N_A}{{\tilde{V}}_m({\tilde{T}}_R)}=\frac{N_A}{V_{m;R}({\tilde{\ell }}_R{)}^3}\mathrm{and}n=\frac{V_{m;R}}{V_{m;T}}.$$

$$p^w=\frac{3nT}{3-n}-\frac{9}{8}{n}^2.$$

$${\mathcal{W}}_q=\frac{Q!}{[H_{Q+1}({\mathcal{R}}_q){]}^2}.$$

$$\omega (p)=\frac{1}{\sqrt{2\pi }}{e}^{-p^2/2},$$

$${\int }_{-\infty }^{+\infty }dp\omega (p)H_{\ell }(p)H_{{\ell }^{\prime }}(p)=\ell !{\delta }_{\ell ,{\ell }^{\prime }}.$$

$${\partial }_t{f}_{\mathbf{\kappa }}+\frac{1}{m}{\mathbf{p}}_{\mathbf{\kappa }}\cdot \nabla f_{\mathbf{\kappa }}+\mathbf{F}\cdot ({\nabla }_{\mathbf{p}}f{)}_{\mathbf{\kappa }}=-\frac{1}{\tau }[f_{\mathbf{\kappa }}-f_{\mathbf{\kappa }}^{eq}],$$

$$f_{\mathbf{\kappa }}^{eq}=n\prod _{\alpha =1}^D{g}_{k_{\alpha }},$$

$$g_{k_{\alpha }}\equiv g_{k_{\alpha }}(u_{\alpha },T)=w_{k_{\alpha }}\sum _{\ell =0}^N{H}_{\ell }(p_{k_{\alpha }})\sum _{s=0}^{\lfloor \ell /2\rfloor }\frac{(mT-1{)}^s(m{u}_{\alpha }{)}^{\ell -2s}}{2^{s}s!(\ell -2s)!},$$

$$\mathbf{F}=\frac{1}{n}\nabla (p^i-p^w)+\sigma \nabla (\Delta n).$$

$$f(\mathbf{x},\mathbf{p},t)=\frac{1}{\sqrt{2\pi }}{e}^{-p_{\alpha }^2/2}\sum _{\ell =0}^{\infty }\frac{1}{\ell !}{\mathcal{F}}_{\alpha ,\ell }(\mathbf{x},p_{\overline{\alpha }},t)H_{\ell }(p_{\alpha }),$$

$${\mathcal{F}}_{\alpha ,\ell }(\mathbf{x},p_{\overline{\alpha }},t)$$

$$\overline{\alpha }$$

$${\partial }_{p_{\alpha }}f(\mathbf{x},\mathbf{p},t)=-\frac{1}{\sqrt{2\pi }}{e}^{-p_{\alpha }^2/2}\sum _{\ell =0}^{\infty }\frac{1}{\ell !}{\mathcal{F}}_{\alpha ,\ell }(\mathbf{x},p_{\overline{\alpha }},t)H_{\ell +1}(p_{\alpha }).$$

$$({\partial }_{p_{\alpha }}f{)}_{\mathbf{\kappa }}\equiv ({\partial }_{p_{\alpha }}f{)}_{\mathbf{\kappa }}(\mathbf{x},t)$$

$${\mathcal{F}}_{\alpha ,\ell ;k_{\overline{\alpha }}}(\mathbf{x},t)$$

$${\partial }_{t}n+\nabla \cdot (n\mathbf{u})=0,$$

$$\rho \left[{\partial }_t{u}_{\alpha }+(\mathbf{u}\cdot \nabla )u_{\alpha }\right]=n{F}_{\alpha }-{\partial }_{\alpha }{p}^i-{\partial }_{\beta }{\Pi }_{\alpha \beta },$$

$$n\left[{\partial }_{t}T+(\mathbf{u}\cdot \nabla )T\right]+{\partial }_{\alpha }{q}_{\alpha }+p^i(\nabla \cdot \mathbf{u})+{\Pi }_{\alpha \beta }{\partial }_{\alpha }{u}_{\beta }=0,$$

$${\Pi }_{\alpha \beta }+p^i{\delta }_{\alpha \beta }=\int d^{2}pf\frac{{\xi }_{\alpha }{\xi }_{\beta }}{m},q_{\alpha }=\int d^{2}pf\frac{{\mathbf{\xi }}^2}{2m}\frac{{\xi }_{\alpha }}{m}.$$

$$\rho \left[{\partial }_t{u}_{\alpha }+(\mathbf{u}\cdot \nabla )u_{\alpha }\right]=n\sigma \nabla (\Delta n)-{\partial }_{\alpha }{p}^w-{\partial }_{\beta }{\Pi }_{\alpha \beta }.$$

$${\Pi }_{\alpha \beta }=-\eta ({\partial }_{\alpha }{u}_{\beta }+{\partial }_{\beta }{u}_{\alpha }-{\delta }_{\alpha \beta }\nabla \cdot \mathbf{u}),q_{\alpha }=-{\kappa }_T{\nabla }_{\alpha }T,$$

$$\eta =\tau nT,{\kappa }_T=\frac{2}{m}\tau nT.$$

$$\mathrm{Pr}=1,$$

$${\partial }_t{f}_{\mathbf{\kappa }}=L[f_{\mathbf{\kappa }}],L[f_{\mathbf{\kappa }}]=-\frac{1}{m}{\mathbf{p}}_{\mathbf{\kappa }}\cdot \nabla f_{\mathbf{\kappa }}-\mathbf{F}\cdot ({\nabla }_{\mathbf{p}}f{)}_{\mathbf{\kappa }}-\frac{1}{\tau }[f_{\mathbf{\kappa }}-f_{\mathbf{\kappa }}^{eq}].$$

$$f_{\mathbf{\kappa }}^{(1)}(\mathbf{x},t)=$$

$$f_{\mathbf{\kappa }}^{(2)}(\mathbf{x},t)=$$

$$f_{\mathbf{\kappa }}(\mathbf{x},t+\delta t)=$$

$$\frac{1}{m}{\mathbf{p}}_{\mathbf{\kappa }}\cdot \nabla f_{\mathbf{\kappa }}=\frac{1}{m}{p}_{k_1}{\partial }_x{f}_{\mathbf{\kappa }}+\frac{1}{m}{p}_{k_2}{\partial }_y{f}_{\mathbf{\kappa }},$$

$${\left(\frac{1}{m}\mathbf{p}\cdot \nabla f\right)}_{\mathbf{\kappa };i,j}=\frac{{\mathcal{F}}_{\mathbf{\kappa };i+1/2,j}-{\mathcal{F}}_{\mathbf{\kappa };i-1/2,j}}{x_{i+1/2,j}-x_{i-1/2,j}}+\frac{{\mathcal{F}}_{\mathbf{\kappa };i,j+1/2}-{\mathcal{F}}_{\mathbf{\kappa };i,j-1/2,}}{y_{i,j+1/2}-y_{i,j-1/2}}$$

$${\mathcal{F}}_{i+1/2}={\overline{\omega }}_1{\mathcal{F}}_{i+1/2}^1+{\overline{\omega }}_2{\mathcal{F}}_{i+1/2}^2+{\overline{\omega }}_3{\mathcal{F}}_{i+1/2}^3,$$

$${\mathcal{F}}_{i+1/2}^1=$$

$${\mathcal{F}}_{i+1/2}^2=$$

$${\mathcal{F}}_{i+1/2}^3=$$

$${\overline{\omega }}_q=\frac{{\tilde{\omega }}_q}{{\tilde{\omega }}_1+{\tilde{\omega }}_2+{\tilde{\omega }}_3},{\tilde{\omega }}_q=\frac{{\delta }_q}{{\varsigma }_q^2}.$$

$${\delta }_1=1/10,{\delta }_2=6/10,{\delta }_3=3/10,$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

Lattice Boltzmann (LB) models with variable temperature are known since at least two decades [1, 2, 3] and their development is still in progress today. Basically, the thermal LB models belong to one of the following families [4, 5]: multi-speed models [6, 7, 8, 9, 10, 11, 12, 13, 14, 15], double distribution function models [16, 17, 18, 19, 20, 21] and hybrid models [22, 23, 24, 25, 26, 27, 28, 29]. Such models are currently applied to investigate physical and engineering processes involving heat transfer with or without phase change, as well as micro- and nano-scale flow phenomena. The diversity of these applications are confirmed by the rich literature related to LB models and by two series of regular conferences [30, 31].

The purpose of this paper is to explore the capabilities of a minimal LB model with variable temperature used to simulate the behaviour of a two-dimensional ($2D$) phase-separating fluid which obeys the Van der Waals equation of state. The model is constructed using the Gauss-Hermite quadrature of order $Q=5$ and has $K=Q\times Q=25$ velocity vectors. Since the roots of the fifth order Hermite polynomial are irrational numbers, the ensuing velocity set is off-lattice. For this reason, finite difference techniques must be employed in order to obtain the numerical solution of the evolution equations in the LB model. In this paper, we employ the third order total variation diminishing (TVD) Runge-Kutta (RK-3) time stepping procedure, together with the fifth-order weighted essentially non-oscillatory (WENO-5) scheme for the advection.

This paper is structured as follows. Section 2 introduces the off-lattice Boltzmann model, with a brief discussion of our non-dimensionalization convention, as well as of the momentum space discretization procedure, the expansion of the equilibrium distribution and the implementation of the forcing term which ensures the recovery of the van der Waals equation of state. The finite difference schemes RK-3 and WENO-5 are introduced in Sec. 3, together with a discussion of the implementation of the boundary conditions. Sec. 4 presents our validation tests (interface width, phase diagram) for the case of a planar interface (when the system is assumed to be homogeneous along the direction parallel to the walls). The Laplace pressure test is performed for circular gas bubbles in Sec. 5. We discuss the transport coefficients appearing in our model and its Galilean invariance in Sec. 6, in the context of the damping of shear and longitudinal waves in a periodic one-dimensional domain, as well as for a droplet in a constant velocity background flow. Section 7 presents our simulation results for the phase separation process induced in a van der Waals fluid at the critical point enclosed between parallel plates which are cooled suddenly. Our conclusions are summarized in Sec. 8.

2 Description of the model

2.1 Non-dimensionalized quantities

In LB models, all quantities of interest are expressed in non-dimesional form. For convenience, in this subsection the tilde ($\widetilde{\hskip 5.0pt}$) symbol over a letter which denotes a physical (measurable) quantity makes the difference between its dimensional form $\widetilde{A}$ and the non-dimensionalized form $A$. These two forms are related through $\widetilde{A}=A\,\widetilde{A}_{R}$, where $\widetilde{A}_{R}$ is the corresponding reference quantity. The non-dimensionalization procedure of the Boltzmann equation [9] amounts to defining four basic reference quantities, namely the particle number density $\widetilde{n}_{R}$, the length $\widetilde{\ell}_{R}$, the mass $\widetilde{m}_{R}$ and the energy $\widetilde{e}_{R}=\widetilde{k}_{B}\widetilde{T}_{R}$, where $\widetilde{k}_{B}$ is the Boltzmann constant and $\widetilde{T}_{R}$ is the reference temperature. In our LB model, the reference length $\widetilde{\ell}_{R}$ is some characteristic length of the fluid system, which may be, e.g., the width of the flow channel. For the single component van der Waals fluid considered in this paper, the mass of its particles is the natural choice for $\widetilde{m}_{R}$, while the properties of the fluid at the critical point are the natural choice for the reference temperature and the reference density, i.e. $\widetilde{T}_{R}=\widetilde{T}_{c}$ and $\widetilde{n}_{R}=\widetilde{n}_{c}$. The values of $\widetilde{T}_{c}$ and $\widetilde{n}_{c}$ defining the critical point can be considered free parameters of the van der Waals model.

The reference values for other physical quantities in the LB model are derived from the above-mentioned basic quantities. In particular, we get the reference speed $\widetilde{c}_{R}=\sqrt{\widetilde{k}_{B}\widetilde{T}_{R}/{\widetilde{m}_{R}}}$ and the reference time $\widetilde{t}_{R}=\widetilde{\ell}_{R}/\widetilde{c}_{R}$ [9, 32]. We choose the reference pressure to be the pressure of the ideal gas at $\widetilde{n}=\widetilde{n}_{c}$ and $\widetilde{T}=\widetilde{T}_{c}$, namely $\widetilde{p}_{R}=\widetilde{n}_{R}\,\widetilde{k}_{B}\widetilde{T}_{R}$. The reference density can also be written in terms of the molar volume $\widetilde{V}_{m}(\widetilde{T}_{R})$ at the reference temperature as follows:

[TABLE]

Here $N_{A}$ is the Avogadro number, while $V_{m;R}$ and $V_{m;T}$ are the nondimensionalized molar volumes at temperatures $\widetilde{T}_{R}$ and $\widetilde{T}$, respectively. With this choice, the non-dimensionalized form of the van der Waals equation reads [12, 32]

[TABLE]

Finally, we note that in our model, the lattice spacing is $\widetilde{\delta s}={\widetilde{\ell}}_{R}/{\mathfrak{N}}$ and hence its non-dimensionalized value is $\delta s=1/{\mathfrak{N}}$, where ${\mathfrak{N}}$ is an integer number.

2.2 Discretization of the momentum space,

evolution equations and the equilibrium distribution functions

In LB models, the non-dimensionalized values of the fluid particle number density $n\equiv n({\bm{x}},t)$, velocity ${\bm{u}}\equiv\bm{u}({\bm{x}},t)$ and temperature $T\equiv T({\bm{x}},t)$, defined in the nodes ${\bm{x}}$ of a lattice ${\mathcal{L}}$, are retrieved at time $t$ through the calculation of the moments (up to second order) of the single particle distribution function $f\equiv f({\bm{x}},{\bm{p}},t)$ defined in the point $({\bm{x}},{\bm{p}})$ of the phase space [4, 11, 33]. Current multispeed LB models use the Cartesian coordinate system in the $D$-dimensional momentum space and the moments of $f({\bm{x}},{\bm{p}},t)$ are computed using a convenient quadrature. For this purpose, $f({\bm{x}},{\bm{p}},t)$ is approximated by its expansion $f^{N}({\bm{x}},{\bm{p}},t)$ up to a certain order $N$ with respect to a set of orthogonal polynomials, e.g., the full-range Hermite polynomials defined on each Cartesian axis of the momentum space [4, 11, 33, 34, 35, 36]. The resulting quadrature points form a discrete vector set $\{{\bm{p}}_{\bm{\kappa}}\equiv(p_{k_{1}},p_{k_{2}},\ldots p_{k_{D}})\}$ in the momentum space (${\bm{\kappa}}\equiv\{k_{1},k_{2},\ldots k_{D}\}$ is a set of integer indices and $p_{k_{\alpha}}$, $1\leq\alpha\leq D$ is the projection of the vector ${\bm{p}}_{\bm{\kappa}}$ on the Cartesian axis $\alpha$). As a result of the application of the Gauss-Hermite quadrature method, in the LB model the fluid system is described by the set of functions $f_{{\bm{\kappa}}}\equiv f_{{\bm{\kappa}}}({\bm{x}},t)=f^{N}({\bm{x}},{\bm{p}}_{\bm{\kappa}},t)$, defined in the nodes ${\bm{x}}$ of a lattice ${\mathcal{L}}$.

In the $D$-dimensional LB model where the full-range Gauss-Hermite quadrature of order $Q$ is used on each Cartesian axis, we have $1\leq k_{\alpha}\leq Q$ for all $\alpha$, $1\leq\alpha\leq D$, and hence the momentum set $\{{\bm{p}}_{{\bm{\kappa}}}\}$ has $K=Q^{D}$ elements. The order $Q$ of the quadrature should satisfy the condition $Q\geq N+1$ in order to retrieve all the moments of $f({\bm{x}},{\bm{p}},t)$ up to order $N$ [11, 34, 35, 37]. Although the number $K$ of the quadrature points can be reduced by very elaborated pruning techniques by sacrificing some higher order moments of the distribution function [11, 38, 39], we will not consider such models in this paper.

When the Bhatnagar-Gross-Krook (BGK) collision term is used in a LB model with variable temperature, the moments of the distribution function $f({\bm{x}},{\bm{p}},t)$ up to order $N=4$ are needed in order to get the evolution equations of the macroscopic fields at the Navier - Stokes - Fourier level [4, 7, 11, 40]. Thus, the minimum number of the momentum vectors in the two-dimensional ($2D$) thermal LB model based on the full-range Gauss-Hermite quadrature that ensures all the moments of $f({\bm{x}},{\bm{p}},t)$ up to order $N=4$ is $K=(N+1)^{2}=25$. In this $2D$ LB model, the full-range Gauss-Hermite quadrature of order $Q=5$ is used on each Cartesian axis. The quadrature points, namely ${\bm{p}}_{\bm{\kappa}}\equiv(p_{k_{1}},p_{k_{2}})$, $1\leq k_{1},k_{2}\leq Q$, are constructed using the direct product rule [11, 35, 36]. For each $\alpha\in\{1,2\}$, the Cartesian projections $p_{k_{\alpha}}$ belong to the set $\{{\cal{R}}_{q}\}$, $1\leq q\leq Q=5$, of the roots of the full-range Hermite polynomial $H_{5}(p)$ [11, 35, 41]. For convenience, Table 1 shows the elements of this set, as well as their associated weights ${\cal{W}}_{q}$ given by [35, 36, 42, 43, 44]

[TABLE]

To avoid confusion, we recall that the full range Hermite polynomials $H_{\ell}(p)$ used in this paper are the so-called probabilistic Hermite polynomials, which are orthogonal with respect to the weight function

[TABLE]

and their orthogonality relation reads [42]

[TABLE]

As usual in the current LB models involving the BGK collision term [11], the non-dimensionalized form of the evolution equation of the functions $f_{\bm{\kappa}}$ for a single-component fluid is:

[TABLE]

where ${\bm{F}}$ is the force acting on a particle of mass $m$ and $\tau$ is the relaxation time. Even though $m=1$ according to the nondimensionalization conventions discussed in Subsec. 2.1, we keep $m$ explicit in all equations in order to avoid confusion. For simplicity, in this paper we assume that the relaxation time $\tau$ is constant. The Cartesian components $(\partial_{p_{\alpha}}f)_{\bm{\kappa}}$, $\alpha\in\{1,2\}$, of the elements in the discrete vector set $\{(\nabla_{\bm{p}}f)_{\bm{\kappa}}\}$ that replace the momentum gradient $\nabla_{\bm{p}}f$ in the Boltzmann equation, will be detailed in Subsec. 2.3.

The equilibrium functions $f^{eq}_{\bm{\kappa}}\equiv f^{eq}({\bm{x}},{\bm{p}}_{\bm{\kappa}},t)$ are given by [35, 36]:

[TABLE]

and $\lfloor\ell/2\rfloor$ is the integer part of $\ell/2$. In the above, it is understood that $g_{k_{\alpha}}\equiv g_{k_{\alpha}}(\bm{x},t)$ since $u_{\alpha}$ and $T$ depend on $\bm{x}$ and $t$. To each $p_{k_{\alpha}}\in\{{\cal{R}}_{q}\}$, $k_{\alpha}=1,\,2,\ldots Q$, $\alpha\in\{1,2\}$, there is an associated weight $w_{k_{\alpha}}\in\{{\cal{W}}_{q}\}$, given by Eq. (3) and we will use the notation $w_{\bm{\kappa}}\equiv w_{k_{1}k_{2}}=w_{k_{1}}w_{k_{2}}$.

2.3 Force term

The following expression of the force $\bm{F}\equiv{\bm{F}}({\bm{x}},t)$, which appears in Eq.(5), was used in order to simulate the evolution of a van der Waals fluid [4, 32, 45, 46, 47, 48]:

[TABLE]

In this expression, the parameter $\sigma$ controls the surface tension, $p^{i}=nT$ is the non-dimensionalized ideal gas pressure and the non-dimensionalised form of the van der Waals pressure $p^{w}$ is given in Eq. (1). The spatial gradients appearing in Eq. (7) are computed using 49-point stencils which are given in Refs. [49, 50]. For the reader’s convenience, these stencils are summarized in Sec. 3.3.

To account for the Cartesian components $(\partial_{p_{\alpha}}f)_{\bm{\kappa}}$, $\alpha\in\{1,2\}$ of $\{(\nabla_{\bm{p}}f)_{\bm{\kappa}}\}$, which appear in Eq.(5), we first expand the single particle distribution function $f({\bm{x}},{\bm{p}},t)$ with respect to the full-range Hermite polynomials $H_{\ell}(p_{\alpha})$ defined on the Cartesian axis $\alpha\in\{1,\,2\}$ of the momentum space, to get [36]:

[TABLE]

where

[TABLE]

Using the recurrence property of the Hermite polynomials $\partial_{x}[e^{-x^{2}/2}H_{\ell}(x)]=-e^{-x^{2}/2}H_{\ell+1}(x)$ [11, 36] to compute the derivative with respect to $p_{\alpha}$ of $f({\bm{x}},{\bm{p}},t)$ given in Eq. (8), we get

[TABLE]

After truncation of this expresion up to order $N$, the application of the discretisation procedure in the momentum space gives:

[TABLE]

Note that the sum in Eq.(14) above runs up to $\ell=N-1$ since $H_{N+1}(p_{k_{\alpha}})=0$ for all $k_{\alpha}=1,\,2,\ldots Q=N+1$, $\alpha\in\{1,2\}$.

2.4 Macroscopic equations

Multiplying the Boltzmann equation (5) with the collision invariants $1$, $\bm{p}$ and $\bm{p}^{2}/2m$ and integrating over the momentum space yields the following macroscopic equations:

[TABLE]

where $\Pi_{\alpha\beta}$ is the viscous part of the stress tensor and $q_{\alpha}$ is the heat flux. This quantities are defined in terms of the peculiar momentum $\bm{\xi}=\bm{p}-m\bm{u}$ as follows:

[TABLE]

The force (7) has the effect of replacing the ideal gas pressure $p^{i}$ in the momentum equation (16b) with the van der Waals pressure $p^{w}$, while also adding a surface tension term:

[TABLE]

The above modification to the momentum equation is sufficient to induce spontaneous phase separation when the temperature $T$ of the fluid is smaller than the critical temperature $T_{c}$.

Furthermore, a Chapman-Enskog analysis shows that, at first order, the viscous stress tensor and the heat flux are given by [40]:

[TABLE]

where the dynamic viscosity $\eta$ and heat conductivity $\kappa_{T}$ have the following expressions:

[TABLE]

The ensuing Prandtl number ${\rm Pr}=c_{p}\eta/\kappa_{T}$ ($c_{p}=2/m$ is the specific heat for a two-dimensional monatomic gas) is fixed in the BGK model at:

[TABLE]

while the hard sphere model predicts that ${\rm Pr}=2/3$.

Given the above mentions there are two remarks we want to highlight: first, since the phase separation mechanism is induced through the use of a body force, the pressure appearing in the energy equation (16c) is still the ideal pressure, instead of the van der Waals pressure [51]; and second, the value of ${\rm Pr}$ (21) is fixed at $1$. The energy equation could be altered such that the ideal fluid pressure is replaced by the van der Waals pressure by employing the modified Boltzmann (i.e. Enskog) equation [52]. Furthermore, there are various methods to alter the value of ${\rm Pr}$, of which we mention the Shakhov [40] and the MRT [53, 54] models. These possible enhancements are the subject of forthcoming work. In this paper, we are interested to perform a first exploration of the capabilities of the single particle distribution function LB model based on Gauss-Hermite quadratures to simulate liquid-vapour thermal flows and, for simplicity, we only considered the simple form of both the body force term and of the collision term in Eq. (5).

3 Numerical scheme and boundary conditions

3.1 Time stepping

In this paper, the time stepping is implemented using the explicit third order total variation diminishing (TVD) Runge-Kutta (RK-3) time marching procedure [55, 56, 57, 58] associated to the fifth-order weighted essentially non-oscillatory (WENO-5) scheme [59, 60] for computing the advection. In order to implement the time stepping algorithm, it is convenient to cast the Boltzmann equation (5) in the following form:

[TABLE]

The third-order TVD Runge-Kutta integrator gives the following algorithm for computing the values of the distribution functions $f_{\bm{\kappa}}$ at time $t+\delta t$:

[TABLE]

3.2 Advection scheme

The advection term which appears in Eq. (22) above, namely

[TABLE]

is computed using the WENO-5 scheme [59, 60] along each Cartesian coordinate. Assuming that the flow domain is discretized using $1\leq i\leq\mathfrak{N}_{1}$ and $1\leq j\leq\mathfrak{N}_{2}$ nodes on the $x$ and $y$ axes, respectively, Eq. (24) becomes:

[TABLE]

where $\mathcal{F}_{\bm{\kappa};i+1/2,j}$ represents the flux of $f$ advected with velocity $p_{k_{1}}/m$ through the vertical interface between the cells centered on $\bm{x}_{i,j}$ and $\bm{x}_{i+1,j}$. Similarly, $\mathcal{F}_{\bm{\kappa};i,j+1/2}$ represents the flux of $f$ advected with velocity $p_{k_{2}}/m$ through the horizontal interface between the cells centered on $\bm{x}_{i,j}$ and $\bm{x}_{i,j+1}$. The construction of these fluxes is summarized below only for the horizontal direction and under the assumption of a positive advection velocity $p_{k_{1}}/m>0$. In this case, the flux $\mathcal{F}_{\bm{\kappa};i+1/2,j}$ can be computed using the following expression [59]:

[TABLE]

where for brevity, the velocity and vertical indices $\bm{\kappa}=\{k_{1},k_{2}\}$ and $j$ were omitted. The interpolating functions $\mathcal{F}^{q}_{i+1/2}$ ($q=1,2,3$) are given by:

[TABLE]

The weighting factors $\overline{\omega}_{q}$ appearing in Eq. (26) are given by:

[TABLE]

where the ideal weights $\delta_{q}$ have the following values:

[TABLE]

The indicators of smoothness $\varsigma_{q}$ can be computed as follows:

[TABLE]

The computation of the weighting factors $\overline{\omega}_{q}$ (28) implies the division between the ideal weights $\delta_{q}$ (29) and the indicators of smoothness $\varsigma_{q}$ (30). In order to avoid division by [math] when either one, two or all three of the indicators of smoothness vanish, it is customary to modify Eq. (30) by adding a small quantity $\varepsilon\simeq 10^{-6}$ to $\varsigma_{q}$. According to Ref. [56], $\varepsilon$ is a dimensionful quantity and its effect on the WENO-5 scheme depends on the typical magnitude of the advected function $f$. It can be seen from Tab. 1 that the ratio between the largest weight $(8/15)^{2}$ and the smallest weight $[(7-2\sqrt{10})/60]^{2}$ is $\sim 2250$, i.e. the set of discrete distributions typically spans three orders of magnitude. Under these circumstances, we follow Refs. [61, 62] and compute the weighting factors $\overline{\omega}_{q}$ directly using Tab. 2 in the limiting cases when any of the indicators of smoothness vanishes.

3.3 High order stencils for the gradient and gradient of the Laplacian

The WENO-5 scheme described in Subsec. 3.2 is of fifth order with respect to the lattice spacing for smooth functions [60]. The smallest square covering all nodes involved in updating a given lattice site comprises $7\times 7=49$ lattice sites. It is therefore natural to seek the computation of the gradient and gradient of the Laplacian appearing in Eq. (7) using the $49$-point stencils described in Refs. [49, 50]. For the reader’s convenience, these stencils are summarized below.

Let $\Delta p=p^{i}-p^{w}$ be the difference between the ideal and van der Waals pressures. Following the discretization of the fluid domain using $\mathfrak{N}_{1}\times\mathfrak{N}_{2}$ equal sized square cells of sides $\delta s$, the function $\Delta p$ is replaced by a set of time-dependent functions $\Delta p_{i,j}$ ($1\leq i\leq\mathfrak{N}_{1}$, $1\leq j\leq\mathfrak{N}_{2}$). In order to compute the first term in Eq. (7) corresponding to the cell $(i,j)$, the gradient of $\Delta p$ can be obtained using the following procedure [49]:

[TABLE]

where $l,m\in\{0,\pm 1,\pm 2,\pm 3\}$. The weights $w^{{\left|l\right|},{\left|m\right|}}$ are symmetric with respect to $l$ and $m$ ($w^{{\left|l\right|},{\left|m\right|}}=w^{{\left|m\right|},{\left|l\right|}}$), having the following values:

[TABLE]

while $w^{0,0}=w^{2,0}=w^{3,2}=w^{3,3}=0$. The resulting scheme is $6$th order accurate and $8$th order isotropic for smooth functions [49].

The second term in Eq. (7) involves the gradient of the Laplacian of the density $n$, which is replaced using the notation introduced above by a set of time-dependent functions $n_{i,j}$ ($1\leq i\leq\mathfrak{N}_{1}$, $1\leq j\leq\mathfrak{N}_{2}$). The computation of the gradient of the Laplacian is performed using the following procedure:

[TABLE]

The convention in the above is that $S_{0,0}$ is the central matrix element (i.e., $S_{1,2}=-c_{7}$). The coefficients $c_{k}$ have the following values:

[TABLE]

while $c_{2}=c_{3}=c_{5}=c_{9}=0$. Eq. (33) is $4$th order accurate and $6$th order isotropic for smooth functios [50].

3.4 Boundary conditions

In this paper, we consider the phase separation in a van der Waals fluid placed between two parallel walls which are perpendicular to the $x$ axis. The flow is always assumed to be periodic along the $y$ axis. As already mentioned in Subsec. 3.3, the flow domain is discretized using a $2D$ square lattice with ${\mathfrak{N}}_{1}\times{\mathfrak{N}}_{2}$ nodes. Let ${\bm{x}}=\delta s[(i-\frac{1}{2})\bm{e}_{1}+(j-\frac{1}{2})\bm{e}_{2}]$, $1\leq i\leq{\mathfrak{N}}_{1}$, $1\leq j\leq{\mathfrak{N}}_{2}$, be the position vectors of the nodes in this lattice. The discussion in this section focuses on the implementation of the specular and diffuse reflection boundary conditions for the distribution functions $f_{\bm{\kappa}}=f_{k_{1},k_{2}}$, as well as for the macroscopic fields involved in the computation of the force term (7).

3.4.1 Specular boundary conditions

During the validation tests considered in Secs. 4 and 5, concerning a plane interface and the Laplace pressure test, the final configuration is considered to be symmetric with respect to the channel centerline, such that the simulation domain can be reduced by implementing specular reflection along the symmetry planes. In particular, let us consider that the center of the channel is located at $(i,j)=(1/2,1/2)$. In order to perform the advection of the distribution function $f_{\bm{\kappa};i,j}=f_{k_{1},k_{2};i,j}$ using Eq. (26), the value of $f_{\bm{\kappa};i,j}$ must be defined below the bottom fluid domain boundary, where $j\in\{-2,-1,0\}$ and $1\leq i\leq\mathfrak{N}_{1}$, as well as to the left of the fluid domain where $i\in\{-2,-1,0\}$ and $1\leq j\leq\mathfrak{N}_{2}$. The specular reflection concept can be implemented along the bottom horizontal boundary as follows:

[TABLE]

On the left vertical boundary, the following procedure can be employed:

[TABLE]

The values $\widetilde{k}_{1}$ and $\widetilde{k}_{2}$ are defined with respect to $k_{1}$ and $k_{2}$, such that the corresponding velocity component is reverted, i.e.:

[TABLE]

3.4.2 Diffuse reflection boundary conditions

Let us now consider the implementation of the diffuse reflection boundary conditions. For definiteness, we refer to the right wall, which is located at $i=\mathfrak{N}_{1}+\frac{1}{2}$. According to the diffuse reflection concept, the distribution function of the fluid particles that return from a plane wall is the Maxwell-Boltzmann equilibrium distribution function corresponding to the wall velocity ${\bm{u}}_{w}$ and temperature $T_{w}$ [35, 40, 63]. This amounts to setting the flux $\mathcal{F}_{k_{1},k_{2};\mathfrak{N}_{1}+\frac{1}{2},j}$ (26) through the interface between the fluid and the wall, located at $i=\mathfrak{N}_{1}+\frac{1}{2}$, to the following value:

[TABLE]

where $f^{eq}_{w;\bm{\kappa}}$ is computed using Eq. (6) by setting $u_{x}=0$, $u_{y}=u_{w}$, $T=T_{w}$ and $n=n_{w}$, where $u_{w}$ represents the vertical velocity of the wall and the wall particle number density $n_{w}$ will be determined below in Eq. (42).

The flux in Eq. (38) can be achieved in the context of the WENO-5 scheme discussed in Subsec. 3.2 by fixing the distribution functions in the ghost nodes at $i\in\{\mathfrak{N}_{1}+1,\mathfrak{N}_{1}+2,\mathfrak{N}_{1}+3\}$ ($1\leq j\leq\mathfrak{N}_{2}$) as follows [62]:

[TABLE]

In order to compute the fluxes of the distributions corresponding to particles travelling towards the wall, two ghost nodes are required inside the wall. The distributions in these ghost nodes are computed using a quadratic extrapolation from the fluid domain, as follows:

[TABLE]

The value of $n_{w}$ in the expression of $f^{eq}_{\bm{\kappa};w}$ is found for each value of $j$ by requiring the total flux of particles to vanish at the wall:

[TABLE]

such that $n_{w}$ can be computed using:

[TABLE]

where $g_{k_{\alpha}}(u_{\alpha},T)$ is defined in Eq. (6b).

3.4.3 Boundary conditions for the macroscopic fields

As discussed in Subsec. 3.3, the computation of the force term (7) is performed using $7\times 7$ stencils. In order to employ these stencils near the fluid domain boundary, the macroscopic fields ($\Delta p=p^{i}-p^{w}$ and $n$) are reflected with respect to the coordinate axes, for both the specular and the diffuse reflection boundary conditions. This reflection has to be performed at the corners of the ghost domain, i.e. for the bottom left corner, where a double reflection occurs as follows:

[TABLE]

where $1\leq i,j\leq 3$.

4 Planar interface

In this section, the capabilities of our thermal models are discussed in the simple case of a planar interface. Complete homogeneity is assumed along the vertical ($y$) direction, such that the flow domain becomes essentially one-dimensional. The number of grid points in this case is $\mathfrak{N}_{1}\times 1$ and the spatial derivatives are computed only along the horizontal direction. In particular, the gradient of $\Delta p=p^{i}-p^{w}$ (31) reduces to:

[TABLE]

Similarly, the stencil (33) for the computation of $\nabla(\Delta n)$ reduces to:

[TABLE]

The analysis presented in this section concerns only the stationary state, in which we assume that the gas phase occupies the central half of the channel, while the liquid phase is confined to the vicinity of the walls. The stationary state is assumed to be symmetric with respect to the center of the channel. Assuming that the time derivatives in Eqs. (16) vanish, Eq. (16a) shows that $u_{x}=0$, while Eq. (16c) shows that $\partial_{x}q_{x}=0$, since the viscous part of the stress tensor $\Pi_{\alpha\beta}$ (19) vanishes. Since the heat flux must vanish at the center of the domain due to symmetry, $\partial_{x}q_{x}=0$ implies $q_{x}=0$ throughout the fluid domain. Finally, Fourier’s law (19) shows that $T=T_{w}$ everywhere inside the fluid domain. The interface shape can be found by solving the stationary limit of Eq. (18):

[TABLE]

An approximate solution due to Wagner and Pooley [64] for the interface profile in the right half of the channel is:

[TABLE]

where $n_{g}$ and $n_{l}$ are the gas and liquid densities, $\xi$ is the interface width and $x_{0}$ is the interface position. This solution loses validity when $T$ is smaller than $1$.

In order to save computational time, only the right half of the channel is simulated, while specular reflection boundary conditions are imposed at the left boundary, as discussed in Sec. 3.4.1. The wall is located at $i_{w}=\mathfrak{N}_{x}+\frac{1}{2}$, where $x_{i_{w}}=0.5$, such that the number of nodes is $\mathfrak{N}_{x}=\mathfrak{N}/2$ and the lattice spacing is $\delta s=1/\mathfrak{N}=1/2\mathfrak{N}_{x}$. The system is initialized with a fluid at constant temperature $T=T_{w}<1$ with the density profile given in Eq. (47), where $n_{g}$ and $n_{l}$ are the liquid and gas densities at $T=T_{w}$ predicted through the well-known Maxwell construction (also known as the equal area construction rule), a procedure described in Refs. [65, 66, 67]. The interface width $\xi$ in Eq. (47) is approximated using the following expression:

[TABLE]

The interface position is set to $x_{0}=0.25$.

Since Eq. (47) is not an exact solution, after the initialization the fluid undergoes an interface adjustment which causes the temperature to rise inside the fluid domain. Due to the extraction of the heat through the lateral walls, the fluid temperature near the wall remains close to $T_{w}$. Around the liquid-gas interface, where the fluid density is not constant, heat generation occurs due to the non-vanishing spurious velocity, which is known to plague multiphase simulations [12, 32, 68, 69, 70, 71, 72, 73]. The stability and accuracy of our simulations is directly improved when the magnitude of these effects is reduced. Thus, Subsec. 4.1 is devoted to the analysis with respect to the lattice spacing and time step of the maximum temperature difference $\Delta T=T-T_{w}$, as well as of the maximum value of the spurious velocity observed in the stationary state. This test is performed at $T_{w}=0.8$, which is considered the working temperature in this paper. After choosing a suitable grid, an analysis of the robustness of our simulations is performed in Subsec. 4.2 by considering various values of $T_{w}$. This analysis is important in order to highlight the validity domain of our simulations. A further validation test is performed in Subsec. 4.3, where the width of the planar interface is discussed. Finally, the phase diagram is discussed in Subsec. 4.4.

4.1 Grid convergence

To find the convergence order of the numerical scheme employed in this paper, we performed a series of simulations with constant time step $\delta t=2\times 10^{-4}$ at $T_{w}=0.8$ and $\sigma=10^{-4}$, for various values of the lattice spacing $\delta s$. Figure 1 shows the general decrease of the temperature difference $\Delta T=T-T_{w}$ and of the spurious velocity $u_{x}$ when the lattice spacing $\delta s$ is decreased. The half-channel profiles of $\Delta T$ and $u_{x}$ are shown in Figures 1(a) and (b) for various values of $\delta s$. It can be seen that both $\Delta T$ and $u_{x}$ exhibit strong oscillations at large lattice spacings (i.e., when $\delta s\gtrsim 1/242$), which are suddenly damped when $\delta s$ decreases under a certain threshold value (i.e., when $\delta s\lesssim 1/244$). A more quantitative analysis of the $\delta s$ dependence of these spurious effects can be made at the level of the maximum temperature differece $(\Delta T)_{\rm max}=T_{\rm max}-T_{w}$ and maximum absolute value of the velocity $|u_{x}|_{\rm max}$. Contrary to expectations, Figs. 1(c) and (d) reveal two exponents. The first corresponds to large values of $\delta s$, when the corresponding profiles are plagued by high amplitude oscillations, and has the values $\sim 3.42$ for $(\Delta T)_{\rm max}$ and $\sim 1.88$ for $|u_{x}|_{\rm max}$. The second exponent refers to the case when the oscillations magnitude is small, having the values $\sim 4.80$ for $(\Delta T)_{\rm max}$ and $\sim 4.14$ for $u_{x;{\rm max}}$.

We thus conclude that, in order to perform accurate simulations at $T_{w}=0.8$ and $\sigma=10^{-4}$, the lattice spacing should be decreased below $1/244$. Unless otherwise stated, we will employ $\delta s=1/320$ for the rest of the simulations presented in this paper.

Before ending this subsection, it is worth mentioning that the convergence with respect to the time step is much less instructive. This is because, for a fixed value of the lattice spacing, the time step is constrained via the CFL condition:

[TABLE]

In particular, for the fifth quadrature order model employed in this paper, $p_{x;\rm max}\simeq 2.86$, such that for $\delta s=1/320$, the time step is constrained via $\delta t<10^{-3}$. Already at this value, the error due to the time integration is smaller than the error introduced by the spatial discretization, such that further decreasing the time step does not seem to significantly improve the numerical results and the convergence test in this particular case is inconclusive. For the rest of this paper, we employ $\delta t=2\times 10^{-4}$.

4.2 Stationary profiles at various temperatures

Figure 2 shows the density, temperature and velocity profiles at various values of the wall temperature $T_{w}$ when $\delta s=1/320$, $\delta t=2\times 10^{-4}$, $\sigma=10^{-4}$ and $\tau=5\times 10^{-3}$. The magnitude of the spurious velocities is less than $10^{-4}$, while the maximum temperature difference is smaller than $5\times 10^{-4}$ even at $T_{w}=0.74$. It can also be seen that for lower values of $T_{w}$, the density gradient along the interface, and hence the amplitude of the spurious velocity and the temperature difference increase, as already observed in Refs. [12, 68].

To better assess the range within which our models can be reliably used for simulations, a series of computer runs were performed with $\sigma=10^{-4}$ by varying the wall temperature between $0.5$ and $0.99$. Figures 3(a) and (b) show the profiles of $\Delta T$ and $u_{x}$ for various values of the temperature $T_{w}$, obtained for $\delta s=1/320$ and $\delta t=2\times 10^{-4}$. A sudden decrease in the magnitudes of $\Delta T$ and $u_{x}$ can be observed when the wall temperature is increased from $T_{w}=0.725$ to $T_{w}=0.726$. This sudden change of magnitude is visible also in Figs. 3(c) and (d), where the maximum values $(\Delta T)_{\rm max}$ and $u_{x;{\rm max}}$ are represented with respect to $1-T_{w}$. These figures show two sets of simulation results, the first corresponding to $\delta s=1/320$ and $\delta t=2\times 10^{-4}$, while the second corresponds to $\delta s=1/512$ and $\delta t=10^{-4}$. Two apparently disjoint regimes can be observed, corresponding to the cases when the fluctuations of the amplitudes are large (small $T_{w}$) or small (large $T_{w}$). As expected, at fixed $T_{w}$, the values of $(\Delta T)_{\rm max}$ and $u_{x;{\rm max}}$ decrease when the lattice spacing is decreased. It can be seen that the point $T_{w}=0.8$ lies inside the region of smaller errors, thus we conclude that the grid spacing $\delta s=1/320$ and the time step $\delta t=2\times 10^{-4}$ are suitable for the simulation of one-dimensional thermal phase separation between parallel plates having the temperature $T_{w}=0.8$.

4.3 Interface width test

In order to study the properties of the interface between the gas and liquid phases, the density profile is investigated in the stationary state, for various values of the wall temperature $T_{w}$ and of the surface tension parameter $\sigma$. In particular, we aim to characterize the interface shape using the approximate formula (47). The gas density $n_{g}$, liquid density $n_{l}$, interface location $x_{0}$ and interface width $\xi$ are obtained by performing a four-parameter nonlinear fit of Eq. (47). The first two plots in Fig. 4 represent the nondimensionalized density $\nu=(2n-n_{l}-n_{g})/(n_{l}-n_{g})$, for: (a) various values of the wall temperature $T_{w}$ and $\sigma=10^{-4}$; (b) various values of $\sigma$ at $T_{w}=0.8$. The values of $n_{g}$ and $n_{l}$ are determined separately for each data set as described above. The numerical results are compared with the approximate formula (47). The legend gives the ratio between the interface width $\xi$ obtained using the nonlinear fit and the value given in Eq. (48). It can be seen that the ratio approaches $1$ as $T_{w}$ approaches the critical temperature.

In order to assess the validity of Eqs. (47) and (48) away from the critical point, Fig. 4(c) shows the relative difference $1-\xi/\xi_{w}$ with respect to the distance from the critical point, measured by $1-T_{w}$. It can be seen that the formula loses validity as the departure from the critical point is increased. In particular, for $T_{w}=0.8$ the relative departure from the predicted value is $7\%$.

Finally, Fig. 4(d) tests the linear dependence of $\xi$ on $\sqrt{\sigma}$ for $T_{w}=0.8$ and $0.9$. This dependence is tested in two steps, as follows. First, the value of $\xi$ is obtained via the four-parameter numerical fit described above. Next, the values of $\xi$ corresponding to various values of $\sigma$ for the same value of $T_{w}$ are used to perform a fit of the law $\xi=\alpha\sqrt{\sigma}$. Excellent agreement is found, while the values of $\alpha$ found at $T_{w}=0.9$ and $0.8$ are within $3.5\%$ and $7\%$ departure from the value predicted through Eq. (48).

4.4 Phase diagram

Figure 5 shows the liquid-vapour phase diagram, as recovered with our model using $\sigma=10^{-4}$ and two sets of values for the lattice spacing and time step, namely $(\delta s,\delta t)\in\{(1/320,2\times 10^{-4}),(1/512,10^{-4})\}$. The values of the density are collected from the first lattice node near the wall for the liquid phase ($n_{l}$) and from the center of the channel for the vapour phase ($n_{g}$). Very good agreement is observed between the density values obtained with our model and those obtained using Eq. (1) via the Maxwell construction. The details of the Maxwell construction are given in Refs. [65, 66, 67].

5 Laplace pressure test

In the case of circular droplets or bubbles, the pressure inside the interface is larger than the pressure outside of it. In this section, we consider the system described in Fig. 6, where a gas bubble of radius $R_{0}$ located at the center of the channel is immersed inside the liquid phase. At initial time, the fluid is assumed to be everywhere in thermal equlibrium ($f=f^{eq}$) at the wall temperature ($T=T_{w}$). The density is initialized according to Eq. (47), where $n_{g}$ and $n_{l}$ are the gas and liquid densities obtained via the Maxwell construction [65, 66, 67] at $T=T_{w}$. The interface width $\xi$ is computed using the approximate formula (48), while the coordinate $x$ is replaced by the distance $R=\sqrt{x^{2}+y^{2}}$ measured from the center of the channel. The initial interface location is at $x_{0}=R_{0}\in\{0.15,0.2,0.25,0.3,0.35\}$.

Assuming that $R_{0}$ is sufficiently large for the bubble to be stable, thermodynamic equlibrium is reached in the stationary state, such that $T=T_{w}$ throughout the domain. The pressure difference between the inside and outside of the interface satisfies:

[TABLE]

where $\gamma$ is the surface tension and $R_{f}$ is the bubble radius after the interface is stabilized.

Since we are interested only in the stationary state which we assume to be symmetric with respect to the horizontal and vertical lines that intersect at the channel center, the computational demand can be reduced by implementing the specular reflection boundary conditions described in Sec. 3.4.1 on the bottom, left and top domain boundaries, as indicated in Fig. 6. The channel center is located at $x=y=0$, while the wall and the top boundary are located at $x=0.5$ and $y=0.5$, respectively. The time step is set to $\delta t=2\times 10^{-4}$ and the lattice spacing is $\delta s=1/320$, such that the simulations are performed on a square domain of size $160\times 160$. For this test case, all simulations were performed with $T_{w}=0.8$ and $\tau=5\times 10^{-3}$.

The density and non-ideal (Van der Waals) pressure are shown in Figs. 7 (a) and (b), respectively, for $\sigma=5\times 10^{-5}$. Good agreement can be observed between the fitted value for $\xi$ ($\sim 0.935\xi_{w}$) and that obtained in Fig. 4(b) for the planar interface ($\sim 0.937\xi_{w}$). The surface tension is computed by multiplying the fitted value $R_{f}$ for the location of the interface by the difference $p^{w}_{\rm in}-p^{w}_{\rm out}$ between the van der Waals pressure (1) inside ($p^{w}_{\rm in}$) and outside ($p^{w}_{\rm out}$) of the bubble. Figure 7(c) represents this difference with respect to $R_{f}^{-1}$, for two values of $\sigma$. A linear fit gives the value of the surface tension $\gamma_{\rm Lap}$.

Alternatively, the surface tension can be computed in the context of a planar interface using the following formula:

[TABLE]

The integration domain is understood to cross only one interface. For this purpose, the profiles presented in Fig. 4(b) can be used to obtain the value $\gamma_{\rm pl}$ for the planar interface for the values of $\sigma$ considered in Fig. 7(c). The numerical results obtained using the Laplace pressure test and the planar interface are summarized in Tab. 3. The relative error $\gamma_{\rm Lap}/\gamma_{\rm pl}-1$ is below $1\%$ for both $\sigma=5\times 10^{-5}$ and $\sigma=10^{-4}$.

6 Transport coefficients, sound speed and Galilean invariance

We further investigate the features of our models in several contexts. First, we demonstrate the ability of our models to correctly recover the transport coefficients of the fluid that we are simulating in the context of the damping of transversal (shear) and longitudinal (sound) waves. Galilean invariance of LB models has been discussed in many studies throughout the past two decades [24, 74, 75, 76, 77, 78, 79, 80]. In this section, we evaluate the sensitivity of our models to Galilei transformations, by considering the wave damping problems at non-vanishing background fluid velocities. We conclude this section by presenting a study of the evolution of a gas bubble enclosed between parallel walls kept at constant temperature for various values of the background velocity. In this section, we use the tilde ($\widetilde{\hskip 5.0pt}$) to denote time-dependent amplitudes. This notation should not be confused with the one introduced in Subsec. 2.1 for dimensional quantities.

6.1 Shear waves

A popular benchmark of lattice Boltzmann models is the damping of shear waves [69, 71, 75, 79]. The system is considered to be homogeneous along the $y$ direction. At initial time, the system is considered to be in local thermal equilibrium at constant density and temperature, while the velocity along the $y$ axis is initialised according to:

[TABLE]

where $k=2\pi/L$ is the wavenumber of the shear wave, $L$ is its wavelength and we use $u_{0}=10^{-3}$. Without loss of generality, we set $L=1$ and choose a number of $\mathfrak{N}_{x}$ cells to discretise the system along the $x$ direction. The coordinate of the center of cell $i$ ($1\leq i\leq\mathfrak{N}_{x}$) is

[TABLE]

such that at initial time, we set $u_{y;i}(t=0)=u_{0}\sin kx_{i}$.

Considering that $u_{0}$ is a small quantity, the continuity, Navier-Stokes and temperature equations (16) reduce to:

[TABLE]

while $n(t)=n_{0}={\rm const}$ and $T(t)=T_{0}={\rm const}$. Assuming that for $t>0$, $u_{y}=\widetilde{u}(t)\sin kx$, Eq. (54) yields

[TABLE]

where the subscript ${\rm lin}$ refers to the analytic solution derived in the linearised limit of the macroscopic equations (16). According to the Galilean invariance of the theory, the solution (55) should be valid also when seen by an observer travelling towards negative values of $x$ with a constant velocity $V$, according to $X(t)=-Vt$, expressed in the laboratory frame. In the inertial frame of this observer, the transverse velocity profile becomes:

[TABLE]

where $\widetilde{u}(t)$ is given in Eq. (55) and $u_{x}=V$ everywhere in the fluid. In order to recover the amplitude $\widetilde{u}(t)$ during our simulations, we use

[TABLE]

Throughout this section, we set the background temperature to $T_{0}=0.8$, $\tau=5\times 10^{-4}$ in order to ensure that our simulations lie in the hydrodynamic regime, $\delta t=2\times 10^{-4}$ and $u_{0}=10^{-3}$ in order to ensure the validity of the linearisation ansatz which leads to Eq. (55). The time variable is discretised using $N_{t}=2\times 10^{5}$ values $t_{n}=n\delta t$ ($1\leq n\leq N_{t}$) in addition to the initial time $t_{0}=0$. We further consider a division with respect to $\Delta t=1000\delta t$, giving rise to $S=200$ values $\hat{t}_{s}=s\Delta t=1000s\times\delta t$ ($1\leq s\leq S$). For each value of $s$, the quantity $\widetilde{u}_{s}$ is computed using a discrete analogue of Eq. (57):

[TABLE]

where $u_{y;s;i}$ is the value of $u_{y;i}$ at time $t=\hat{t}_{s}=1000s\times\delta t=0.2s$.

In order to perform quantitative analyses, the values $\widetilde{u}_{s}$ are used to perform a numerical fit of Eq. (55) which allows the parameter $\nu_{\rm app}$ to be extracted, using which the apparent viscosity can be computed via $\eta_{\rm app}=\rho_{0}\nu_{\rm app}/k^{2}$. The second type of quantitative analysis concerns the $L_{2}$ norm of the relative difference between $\widetilde{u}_{s}$ and the expected value $\widetilde{u}(t)$ (55), which we compute using the trapezoidal rule:

[TABLE]

where $\mathfrak{f}_{s}=0.5$ when $s=0$ and $s=S$ and $\mathfrak{f}_{s}=1$ for $1\leq s<S$.

We consider three batches of simulations. The first corresponds to the case of the ideal gas at unit density ($n_{0}=n_{0;{\rm ideal}}=1$), when the forcing term in Eq. (5) is not taken into account. The second and third batches correspond to the cases of the van der Waals vapour ($n_{0}=n_{0;{\rm g}}\simeq 0.2397$) and liquid ($n_{0}=n_{0;{\rm l}}\simeq 1.933$) phases, respectively. For each simulation batch, we consider discretisations with $\mathfrak{N}_{x}=20$, $30$, $40$, $50$ and $60$ points. For each value of $\mathfrak{N}_{x}$, we consider velocities $V$ ranging from [math] (laboratory frame) to $2$, with a step of $0.1$. Since in this problem, the gradients of the density and temperature (and hence, of the pressure) are expected to vanish, the numerical results for these three media are very similar.

Figure 8 shows in the top panels the typical time dependence of the amplitude $\widetilde{u}(t)$, while in the bottom panels, the convergence tests are presented, as discussed below. For simplicity, only the results for the first batch of simulations are shown (the case of the ideal gas).

In panel (a) of Fig. 8, the longitudinal velocity is $V=0$ (laboratory frame) and the domain is discretised using various number of nodes. It can be seen that already at $\mathfrak{N}_{x}=30$, a reasonable agreement is found compared to the analytic formula (55). In panel (b), various values of the overall longitudinal velocity $V$ are considered, while keeping $\mathfrak{N}_{x}=20$ in order to enhance the differences between the various numerical results. It can be seen that the results deteriorate as $V$ is increased.

In the bottom panels of Fig. 8, the error in $[(\eta_{\rm app}/\eta_{0})-1]$ and in the $L_{2}$ norm computed using Eq. (59) are presented.

In panel (c) of Fig. 8, the longitudinal velocity is set to $V=0$ and the number of nodes $\mathfrak{N}_{x}$ is varied. A numerical fit of $[(\eta_{\rm app}/\eta_{0})-1]$ and $L_{2}$ as functions of $a(\delta x)^{\gamma}$ gives values of $\gamma$ close to $5$, confirming that the WENO-5 scheme employed in this paper is fifth-order accurate, as also shown in Ref. [60]. The analysis discussed above is performed for overall longitudinal velocities $V=0$, $0.5$, $1$ and $2$ for the ideal gas and for the vapour and liquid phases of the van der Waals fluid and the results are reported in Tab. 4. It can be seen that the decrease of the exponents with $V$ is insignificant (less than $2\%$ difference between the cases $V=0$ and $V=2$). As expected, the exponents at fixed values of $V$ are very similar for the three media considered herein.

Finally, panel (d) of Fig. 8 measures the effects of increasing the longitudinal velocity on the numerical viscosity and the $L_{2}$ norm. The number of nodes is kept fixed at $\mathfrak{N}_{x}=30$. Since at $V=0$, the error compared to the analytic estimates is finite, the influence of $V$ can be isolated by considering the numerical results for ${\eta_{\rm app;V}}/\eta_{0}$ and $L_{2;V}$ obtained at finite $V$ relative to their values when $V=0$. A numerical fit of the scaling law $aV^{\gamma}$ shows that the differences $(\eta_{{\rm app};V}-\eta_{{\rm app};0})/\eta_{0}$ and $L_{2;V}-L_{2;0}$ grow with exponent $\gamma\simeq 2$. Since this scaling holds only for small values of $V$, the fits are performed for $0<V\leq 0.5$. Further results for this test are shown in Tab. 5, where the number of grid points $\mathfrak{N}_{x}$ is varied from $20$ to $60$. It can be seen that the exponent $\gamma$ is very close to $2$ for all tested cases. As expected, the difference between the results obtained for the ideal gas and the vapour and liquid phases of the van der Waals fluid is negligible.

6.2 Longitudinal waves

We now turn to another important problem in fluid dynamics which concerns the study of longitudinal waves [80, 81, 82, 83, 84, 85, 86, 87]. The propagation of longitudinal waves induces fluctuations in the macroscopic properties of the fluid, the amplitudes of which decay due to viscous and thermal dissipation. Considering that the wave propagates through a background state characterised by $n_{0}$, $T_{0}$ and $u_{x;0}=V$, the macroscopic quantities can be written as:

[TABLE]

The background velocity $V$ is introduced above to enable the solution of the longitudinal wave problem to be considered in any Galilean frame in motion along the wave’s direction of propagation. For simplicity, we set $V=0$ in the following, since $V$ can be restored at the end of the calculation using the principle of Galilean invariance.

Taking the limit when the wave amplitudes $\delta n$, $\delta T$ and $\delta u$ are small, the macroscopic equations (16) can be linearised as follows:

[TABLE]

We remind the reader that the heat capacity at constant volume is $C_{V}=1$ in our two-dimensional framework, corresponding to an adiabatic index $\gamma=2$. In the above, $\rho_{0}=mn_{0}$, $\sigma$ is the surface tension parameter, while $\eta_{0}=\tau n_{0}T_{0}$ and $\kappa_{T;0}=2\tau n_{0}T_{0}/m$ are the shear viscosity and thermal conductivity for the BGK model discussed in Sec. 2.4. The pressures $P_{\eta}$ and $P_{\kappa}$ appearing in the Navier-Stokes and heat equations are left unspecified in order to allow the same framework to be applied for the ideal and van der Waals fluids. Their values at $n=n_{0}$ and $T=T_{0}$ are denoted by $P_{\eta;0}$ and $P_{\kappa;0}$, while the perturbation $\delta P_{\eta}=P_{\eta}-P_{\eta;0}$ can be written as:

[TABLE]

Since all the relations in Eq. (61) are linear and homogeneous with respect to the perturbation amplitudes, a harmonic decomposition can be made. Let $k=2\pi/L$ be the wave number of a longitudinal wave with wavelength $L$. In the laboratory frame ($V=0$), the following ansatz can be made:

[TABLE]

where the amplitudes $\widetilde{\delta u}$, $\widetilde{\delta n}$ and $\widetilde{\delta T}$ depend only on time. The analysis of the wave damping can be made at the level of independent modes by writing these amplitudes as follows:

[TABLE]

where the coefficients $\delta u_{\nu}$, $\delta n_{\nu}$ and $\delta T_{\nu}$ are constant.

Substituting Eqs. (63) and (64) into Eq. (61) gives:

[TABLE]

Setting $\delta u_{\nu}=0$ yields the trivial solution $\delta n_{\nu}=\delta T_{\nu}=0$. Thus, non-trivial dynamics occur only when:

[TABLE]

The above equation is cubic with respect to $\nu$, thus it admits at least one real solution, which corresponds to the thermal mode $\nu_{t}$. The other two roots $\nu_{\pm}$, corresponding to the acoustic modes, must be complex in order to allow the wave to propagate. Writing $\nu_{\pm}=\nu_{a}\pm i\nu_{s}$, we see that $\nu_{a}$ induces acoustic dissipation, while $\nu_{s}=kc_{s;0}$ is related to the speed of sound $c_{s;0}$ at the background parameters. In order to derive expressions for the allowed values of $\nu$, we remember that $\eta_{0}$ and $\kappa_{T;0}$ are proportional to $\tau$, which is assumed to be a small number in order for the flow to remain within the hydrodynamic regime. Thus, we seek solutions of the form:

[TABLE]

where the imaginary unit $i$ was inserted in front of the leading order term. Inserting Eq. (67) into Eq. (66) gives, to orders $\tau^{0}$ and $\tau^{1}$,

[TABLE]

where we have identified the speed of sound as:

[TABLE]

Setting $\nu_{0}=0$ yields the thermal mode $\nu_{t}=\tau\nu_{1}$, while the cases with $\nu_{0}\neq 0$ correspond to the acoustic modes $\nu_{\pm}=\nu_{a}\pm i\nu_{s}$:

[TABLE]

In the case of the ideal gas, we have $P_{\eta;0}=P_{\kappa;0}=p^{i}_{0}=n_{0}T_{0}$, hence

[TABLE]

The above expressions can be seen to coincide with those deduced in Eq. (39) of Ref. [80] by setting $\gamma=2$, ${\rm Pr}=1$, $\kappa=\kappa_{T}/\rho_{0}c_{p}$, $c_{p}=2/m$ and $\kappa_{T}=2\tau n_{0}T_{0}/m$.

The behaviour of the true van der Waals fluid can be recovered by setting $P_{\eta;0}=P_{\kappa;0}=p^{w}$. For the model used in this paper, only $P_{\eta;0}=p^{w}$ is fulfilled, while $P_{\kappa;0}=p^{i}$ is approximated through the ideal gas pressure. This latter approximation affects the speed of sound $c_{s}$, as well as the damping coefficients $\nu_{a}$ and $\nu_{t}$. Figure 9 presents the dependence of $c_{s}$, $\nu_{a}$ and $\nu_{t}$ with respect to $n_{0}$ for the ideal gas, for the true van der Waals fluid and for our model. The most important feature for the study of phase separation phenomena is the correct recovery of the spinodal curve, which is given by the solutions of $\nu_{t}=0$ for various values of $T$. It can be seen from Eq. (70) that the spinodal curve is determined by solving

[TABLE]

The surface tension term reduces the breadth of the spinodal region. For $k=2\pi/L$ and $\sigma=10^{-4}$, the contribution of this term is negligible. However, there is always a minimal wavelength under which spinodal decomposition cannot occur, namely $L_{\rm min}=2\pi/k_{\rm min}$, where

[TABLE]

At $\sigma=10^{-4}$ and $T=0.8$, $L_{\rm min}\simeq 0.0966$, while the interface width predicted through Eq. (48) is $\xi_{w}\simeq 0.021$. For large wavelengths, the surface tension can be neglected and the equation $\partial_{n}p^{w}=0$ predicts that spinodal decomposition can occur for densities between $0.521$ and $1.574$. Another interesting feature of the real van der Waals fluid is that the speed of sound becomes imaginary for densities between $0.771$ and $1.631$ (the values are computed for $T_{0}=0.8$ and $\sigma=10^{-4}$), which lie within the spinodal curve. In our model, the speed of sound remains real for all values of $n$. In the regions outside the spinodal curve, the qualitative behaviour of the speed of sound in our model and in the true van der Waals model is similar. There are however discrepancies between the value of the speed of sound in the vapour and liquid phases for the true van der Waals fluid ($c_{s;{\rm vapour}}\simeq 1.0289$ and $c_{s;{\rm liquid}}\simeq 1.479$) compared to those occuring in our model ($c_{s;{\rm vapour}}\simeq 1.130$ and $c_{s;{\rm liquid}}\simeq 2.056$).

Let us now write the solution of the longitudinal wave problem. Considering the coefficients $\delta u_{t}$ and $\delta u_{\pm}$ of the $\nu_{t}$ and $\nu_{\pm}$ modes as independent variables, the solutions $\widetilde{\delta u}$, $\widetilde{\delta n}$ and $\widetilde{\delta T}$ can be written as:

[TABLE]

where Eq. (65) can be used to obtain:

[TABLE]

In the case when the density and temperature perturbations vanish at initial time (i.e., $\delta n_{0}=\delta T_{0}=0$), the constants $\delta u_{c}$, $\delta u_{t}$ and $\delta u_{s}$ are given up to $\tau^{2}$ through:

[TABLE]

In this case, $\delta n_{t}\simeq 0$ and $\delta P_{t}\simeq 0$, such that the contribution of the purely evanescent mode $\nu_{t}$ is negligible. The full solution for the velocity amplitude can be written up to $O(\tau^{2})$ as:

[TABLE]

We now present our simulation results. As in Subsec. 6.1, we perform three batches of simulations: the first is for the ideal gas with $n_{0}=n_{\rm ideal}=1$; the second and third are for the vapour ($n_{0}=n_{g}=0.2396669$) and liquid ($n_{0}=n_{l}=1.932703$) phases of the van der Waals fluid with $\sigma=10^{-4}$, respectively. All simulations are performed at $T=0.8$ and we consider the wavelength fixed at $L=1$ ($k=2\pi/L$). The relaxation time is fixed at $\tau=5\times 10^{-4}$ while the time step is $\delta t=2\times 10^{-4}$. The initial wave amplitudes are $\delta u_{0}=10^{-3}$, $\delta n_{0}=\delta T_{0}=0$ and the evolution of the velocity amplitude $\widetilde{\delta u}$ is given in Eq. (77). In each simulation batch, we use between $N_{x}=20$ and $N_{x}=60$ grid points and the background velocity along the $x$ axis is varied between $V=0$ and $V=2$.

In all cases, we perform $N_{t}=50000$ iterations, up to $t_{f}=10$. The values of the amplitude of the velocity are stored at intervals $\Delta t=100\delta t$ and are labelled $\widetilde{\delta u}_{s}$ ($0\leq s\leq S=500$). The procedure for computing $\widetilde{\delta u}_{s}$ is:

[TABLE]

where $\hat{t}_{s}=s\times\Delta t=100s\times\delta t$. The quantitative analysis is performed at the level of three quantities, namely: the acoustic damping coefficient $\nu_{a}$; the sound speed $c_{s}$; and the $L_{2}$ norm. In order to extract $\nu_{a}$ and $c_{s}$ from the numerical data, Eq. (77) is written as:

[TABLE]

The parameters $\nu_{a}$, $c_{s}$ and $\mathcal{S}$ are obtained by performing a three-parameter fit of Eq. (79) for the case of the van der Waals fluid, while in the case of the ideal gas, $\mathcal{S}$ is set to [math] and the fit is performed using only two free parameters. The other parameters are $\delta u_{0}=10^{-3}$ and $k=2\pi$. The $L_{2}$ norm is computed as follows:

[TABLE]

where $\widetilde{\delta u}_{\rm lin}$ is the solution of the linearised hydrodynamic equation derived in Eq. (77), evaluated at $t=\hat{t}_{s}$.

Figure 10 shows the typical evolution of the amplitude $\widetilde{u}(t)$, as obtained using our numerical method, compared to the analytic solution (77) of the linearised hydrodynamics equations, for the cases of the ideal gas and of the vapour and liquid phases of the van der Waals fluid at $T=0.8$. It can be seen that our numerical results are well overlapped with the analytic solution even when $\mathfrak{N}_{x}=20$ and $V=2$, thus demonstrating the capabilities of the numerical scheme and the degree of Galilean invariance of our implementation. Figures 11 and 12 describe the typical procedure that we employed for the quantitative analyses discussed below.

In Fig. 11, the relative error of the numerically obtained values for $\nu_{a}$ and $c_{s}$, computed with respect to their analytic expectations in Eq. (70), as well as $L_{2}$, are represented with respect to $\delta x=1/\mathfrak{N}_{x}$ for the case when $V=0$. It can be seen that the errors decrease with $\delta x$ only for $\delta x\gtrsim 0.04$. Further decreasing $\delta x$ allows these quantities to stabilize at values which have a relative difference compared to Eq. (70) of about $10^{-3}$. This difference is comparable to both $\delta u_{0}$ and $\tau$, while Eqs. (70) and (77) are valid only at linear order in these quantities. Thus, the deviatation of the numerical results from the results obtained in the linearised regime is consistent with the assumptions employed in deriving the analytic solution. Thus, in order to extract the order of the numerical scheme, a numerical fit of the function $a(\delta s)^{\gamma}$ is performed on the relative errors, but only for $\delta s\geq 1/30$, thus avoiding the effects of the plateau region which appears at smaller values of $\delta s$. The exponents of the above mentioned fits are given in the legend. Further analysis was performed by considering a selection of values for $V$, between $V=0$ and $V=2$ and the results are summarised in Tab. 6. The values of the exponents corresponding to $\nu_{d}$ and $L_{2}$ are generally confined between $4$ and $5$. The exponent corresponding to $c_{s}$ presents wider variations, having values larger than $5$ for $V\lesssim 1$, and decreasing below $2$ at $V=2$. It should be noted that the relative error in obtaining $c_{s}$ is several orders of magnitude below $(\nu_{a;{\rm app}}/\nu_{a;{\rm lin}}-1)$ and $L_{2}$. We conclude that our numerical results generally support that our numerical scheme is at least of order $4$ for small and moderate values of $V$.

Figure 12 summarises the procedure that we used in order to determine the order at which changing the background longitudinal velocity $V$ affects our numerical results. Thus, we considered the relative errors of the differences of the numerically determined quantities $\nu_{a}$, $c_{s}$ and $L_{2}$ for a given value of $V$ and their values obtained when $V=0$. It can be seen that for sufficiently small values of $V$, a power law can be observed, the exponent of which indicates an almost quadratic dependence of these differences on $V$. The results presented in Fig. 12 are restricted to the case when $\mathfrak{N}_{x}=60$. Further results for various values of $\mathfrak{N}_{x}$ between $20$ and $60$ are summarised in Tab. 7. These results are generally supportive of the nearly quadratic dependence of the relative errors on $V$.

6.3 Laplace pressure test of a moving bubble

In this subsection we present a study of the evolution of a gas bubble enclosed between parallel walls kept at constant temperature for various values of the background velocity parallel to the system walls. We initialised the system with a bubble centered on $x=0,y=0$ using formula (47), temperature $T=T_{w}=0.8$, fluid velocity $u_{y}=V$ and wall velocity $u_{y}^{w}=V$. We test the Galilean invariance of our model by tracking the value $\gamma(V,t)$ of the surface tension evaluated using the Laplace law (50), as well as of the vertical coordinate $y_{c}(V,t)$ of the bubble center, evaluated using formula (47), for background velocities $V=0.01,0.02,0.05$ and $0.1$.

In Fig.14 (a) we present the evolution of the relative error $\epsilon_{\gamma}(V,t)=|\gamma(V,t)/\gamma(0,t)-1|$, where $\gamma(0,t)$ is the value of the surface tension for the stationary bubble. It can be observed that this error is well below $1\%$ for all the velocities considered. The numerical fit of the average values $\overline{\epsilon}_{\gamma}(V)$, over the interval $20\leq t\leq 30$, against the function $aV^{b}$, as seen in Fig. 14(b), yields a value of $b$ very close to $2$. In Fig. 14(c) we plot the evolution of the relative deviation $\epsilon_{y}(V,t)=[y_{c}(V,t)/Vt]-1$ of the bubble center. In Fig. 14(d) we plot the deviation $\Delta y_{c}(V,t)=y_{c}(V,t)-Vt$ with respect to background velocity at $t=30$. A numerical fit gives a nearly linear dependence of the deviation $\Delta y_{c}(V,30)$ on the background velocity.

In order to better illustrate the deviations of the bubble center in time, we present in Fig. 15 the density isocontour corresponding to $n=1.0$ for $V=0.1$ after one ($t=10$), two ($t=20$) and three ($t=30$) cycles. It can be seen that the bubble constantly lags behind the background flow, such that the contours corresponding to successive cycles do not overlap. We attribute this effect to the spurious currents which are always present at the interface.

7 Phase separation dynamics between heat extracting parallel plates

We now consider the phase separation in a van der Waals fluid placed between two parallel walls. The fluid between the walls is in isothermal conditions at the critical temperature $T=1$ and the density field is initialized with small fluctuations $n_{i,j}=1+\delta n_{i,j}$, where the indices $i$ and $j$ identify the grid node, while $\delta n_{i,j}$ is a random point-dependent number satisfying $-0.01\leq\delta n_{i,j}\leq 0.01$. At initial time, the temperature of the walls is suddenly decreased to $T_{w}<1$ (in this section, we only consider $T_{w}=0.8$), thus inducing the phase separation process by cooling the system via the gradual extraction of heat through the diffuse reflecting walls. Figure 16 shows the initial setup of the problem. The temperature of the walls is kept constant throughout the simulation. The simulation domain is comprised of $2\mathfrak{N}\times\mathfrak{N}$ nodes. The lattice spacing is set to $\delta s=1/\mathfrak{N}$, where $\mathfrak{N}=320$, such that the left and right walls are located at $x=-1$ and $x=1$, respectively. The domain has unit vertical span. Diffuse reflection boundary conditions are implemented along the walls, which are located at $i=\frac{1}{2}$ and $i=2\mathfrak{N}+\frac{1}{2}$ ($1\leq j\leq\mathfrak{N}$). Periodic boundary conditions apply along the top and bottom domain boundaries, where $j=\frac{1}{2}$ and $j=\mathfrak{N}+\frac{1}{2}$ ($1\leq i\leq 2\mathfrak{N}$).

Figure 17 shows the evolution of the liquid - vapor separation process on a $640\times 320$ lattice. At $t=0$, the fluid temperature was $T=T_{c}=1.0$ and the wall temperature was set to $T_{w}=0.8$. The simulation was conducted with the time step and lattice spacing $\delta t=2\times 10^{-4}$ and $\delta s=1/320$, respectively, while $\sigma=10^{-4}$ and $\tau=5\times 10^{-3}$. At early times, one observes the liquid deposition on the cold walls. As the bulk temperature decreases, further parallel bands of low and high density appear near the walls. Afterwards, these bands break into individual droplets due to the action of surface tension. The formation of liquid droplets in the central region of the channel is observed at later stages of the simulation. This happens because the temperature in the center of the channel decreases during the heat extraction through the walls, but always remains higher than the wall temperature, as seen in the right column of Figure 17 . This feature was observed also when investigating the liquid-vapour phase separation with a different thermal LB model [12]. Moreover, the right column of Figure 17 revealed that the local maxima in the temperature profiles are located in the liquid-vapour interface regions. In these regions, where large density gradients are present, there is an unphysical heat generation process, due to the spurious velocity, a numerical effect that plagues the LB models [5, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97]. This numerical effect has been succesfully reduced in this paper by using the fifth order weighted essentially non-oscillatory (WENO-5) numerical scheme.

8 Conclusion

A single particle distribution function thermal lattice Boltzmann model based on the full-range Gauss-Hermite quadrature of order $Q=5$ was tested by simulating the liquid-vapour phase separation in a van der Waals fluid bounded by two parallel walls. The Van der Waals force term was implement using 49-point stencils. We validated our thermal model by considering the plane interface problem, the Laplace pressure test and by comparing the phase separation results against the Maxwell construction results. Good agreement was obtained for temperatures as low as $0.72\,T_{c}$. We also present a discussion on transport coefficients, sound speed and Galilean invariance.

Starting from an initial state in which the fluid is at the critical temperature $T_{c}$, with random density fluctuations of at most $1\%$ around the critical density, we investigated the phase separation between two walls kept at a constant temperature $T_{w}<T_{c}$. Our simulations show that the condensation starts in the vicinity of the walls. As the bulk temperature decreases, liquid droplets develop further in the channel. We observed that in the stationary state, spurious currents persist, which affect the temperature field through spurious heating at the vapour-liquid interface. To avoid the large relative errors induced by the numerical effects, it is necessary to use high order schemes such as the fifth order weighted essentially non-oscillatory (WENO-5) scheme with sufficiently small values of the lattice spacing $\delta s$ and the time step $\delta t$, e.g. $(\delta s,\delta t)=(1/160,2\times 10^{-4})$. With these parameters, the non-dimensionalised magnitude of the spurious velocity is below $6\times 10^{-5}$, while the fluid temperature profile is within the range of $1.0\%$ above the wall temperature $T_{w}$, even when $T_{w}=0.72\,T_{c}$.

Acknowledgments

This work is supported by a grant from the Romanian National Authority for Scientific Research, CNCS-UEFISCDI, project number PN-II-ID-PCE-2011-3-0516. The authors are indebted to Adrian Horga for invaluable insight regarding the development of our CUDA code.