An Introduction to Computational
Fluctuating Hydrodynamics

Our first example is the stochastic heat equation (SHE) for the transport of internal energy by thermal diffusion. This section presents a short derivation, for details see [2]. We start with the continuity equation for energy

$$\frac{\partial}{\partial t}(\rho e)=-\nabla\cdot\boldsymbol{Q}$$

(1)

where $\rho$ is the (constant) mass density and $e$ is the specific internal energy. The heat flux, $\boldsymbol{Q}$, can be separated into deterministic and stochastic contributions (denoted by $\overline{\boldsymbol{Q}}$ and $\widetilde{\boldsymbol{Q}}$, respectively) and is written in Onsager form as

$$\boldsymbol{Q}=\overline{\boldsymbol{Q}}+\widetilde{\boldsymbol{Q}}=L\boldsymbol{X}+\widetilde{\boldsymbol{Q}}$$

(2)

where $L$ is the Onsager coefficient and $\boldsymbol{X}$ is the thermodynamic force. Fluctuation-dissipation gives that the stochastic flux has zero mean and correlation

$$\langle\widetilde{\boldsymbol{Q}}(\mathbf{r},t)\widetilde{\boldsymbol{Q}}(\mathbf{r}^{\prime},t^{\prime})\rangle=2k_{B}L~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime})$$

(3)

where $k_{B}$ is the Boltzmann constant and $\delta(\cdot)$ denotes the Dirac delta function. Readers familiar with the Langevin equation for Brownian motion [3] will notice the similarities between that stochastic ordinary differential equation (ODE) and the above stochastic PDE.

From non-equilibrium thermodynamics [2, 4], the deterministic total rate of entropy change in a region $\Omega$ is

$$\frac{d\overline{S}}{dt}=\int_{\Omega}\frac{d\overline{s}}{dt}\;d\mathbf{r}=\int_{\Omega}\boldsymbol{X}\cdot\boldsymbol{J}\;d\mathbf{r}-\left[\frac{\overline{\boldsymbol{Q}}}{T}\right]_{\partial\Omega}$$

(4)

where $\boldsymbol{J}$ is the thermodynamic flux. The first term is the internal entropy production and the second is the change due to heat flow at the boundary. Using the Gibbs equation, $d\overline{u}=\rho d\overline{e}=Td\overline{s}$, and Eq. (1) gives

$$\int_{\Omega}\frac{d\overline{s}}{dt}\;d\mathbf{r}=\int_{\Omega}\frac{\rho}{T}\frac{d\overline{e}}{dt}\;d\mathbf{r}=-\int_{\Omega}\frac{1}{T}\nabla\cdot\overline{\boldsymbol{Q}}\;d\mathbf{r}=\int_{\Omega}\frac{1}{T}\nabla\cdot(\lambda\nabla T)\;d\mathbf{r}.$$

(5)

The last equality uses the phenomenological Fourier law for heat flow, $\overline{\boldsymbol{Q}}=-\lambda\nabla T$, where $\lambda$ is the thermal conductivity. Applying Green’s first identity gives

$\displaystyle\int_{\Omega}\frac{1}{T}\nabla\cdot(\lambda\nabla T)\;d\mathbf{r}$

$\displaystyle=\left[-\frac{\overline{\boldsymbol{Q}}}{T}\right]_{\partial\Omega}-\int_{\Omega}\nabla\left(\frac{1}{T}\right)\cdot(\lambda\nabla T)\;d\mathbf{r}.$

(6)

The entropy expressions (4) and (5) can be equated, and so the chain rule gives

$\displaystyle\int_{\Omega}\boldsymbol{X}\cdot\boldsymbol{J}\;d\mathbf{r}$

$\displaystyle=-\int_{\Omega}\nabla\left(\frac{1}{T}\right)\cdot(\lambda\nabla T)\;d\mathbf{r}=\int_{\Omega}\nabla\left(\frac{1}{T}\right)\cdot\left(\lambda T^{2}\;\nabla\frac{1}{T}\right)\;d\mathbf{r}.$

(7)

From this equality, $\boldsymbol{X}\cdot\boldsymbol{J}=\boldsymbol{X}\cdot(L\boldsymbol{X})$, and Eq. (2), we can identity the thermodynamic force $\boldsymbol{X}$, the thermodynamics flux $\boldsymbol{J}$, and the Onsager coefficient $L$ to be [4, 5]

$$\boldsymbol{X}=\nabla\left(\frac{1}{T}\right)\;\mathrm{,}\;\;\;\boldsymbol{J}=\lambda T^{2}\;\nabla\left(\frac{1}{T}\right)=\overline{\boldsymbol{Q}}\;,\;\;\;\mathrm{and}\;\;\;L=\lambda T^{2},$$

(8)

respectively. The noise correlation is then given by¹¹1Formally, the temperature in the noise should be the deterministic value but in CFHD the instantaneous temperature is typically used.

$$\langle\widetilde{\boldsymbol{Q}}(\mathbf{r},t)\widetilde{\boldsymbol{Q}}(\mathbf{r}^{\prime},t^{\prime})\rangle=2k_{B}\lambda T^{2}~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime}).$$

(9)

Finally, by writing $e=c_{V}T$, where $c_{V}$ is the specific heat at constant volume, the SHE is given as

$$\rho c_{V}\frac{\partial T}{\partial t}=\nabla\cdot(\lambda\nabla T+\sqrt{2k_{B}\lambda T^{2}}~{}\widetilde{\boldsymbol{Z}})$$

(10)

where $\widetilde{\boldsymbol{Z}}$ is a Gaussian white noise vector field uncorrelated in space and time, that is,

$$\langle\widetilde{\boldsymbol{Z}}_{s}(\mathbf{r},t)\widetilde{\boldsymbol{Z}}_{s^{\prime}}(\mathbf{r}^{\prime},t^{\prime})\rangle=~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime})~{}\delta^{\mathrm{Kr}}_{s,s^{\prime}}$$

(11)

where $s\in\{x,y,z\}$ is a Cartesian coordinate direction and $\delta^{\mathrm{Kr}}$ is Kronecker delta.

For simplicity we focus on the one-dimensional case (e.g., thin rod) and write the SHE in compact form as

	$\displaystyle\partial_{t}T$	$\displaystyle=\partial_{x}\left(\kappa\partial_{x}T+\alpha T\tilde{Z}\right)$		(12)
		$\displaystyle=\kappa\partial_{x}^{2}\,T+\alpha\partial_{x}\,T\tilde{Z}$		(13)

where $\kappa=\lambda/\rho c_{V}$ and $\alpha=\sqrt{2k_{B}\lambda}/\rho c_{V}$ are taken to be constants. In the one-dimensional case the noise variance is

$$\langle\tilde{Z}(x,t)\tilde{Z}(x^{\prime},t^{\prime})\rangle=\frac{1}{\mathcal{A}}~{}\delta(x-x^{\prime})~{}\delta(t-t^{\prime})$$

(14)

where $\mathcal{A}$ is the cross-sectional area of the system.

Equation (13) can be numerically integrated in time to produce sample trajectories for $T(x,t)$. We discretize space and time as $x_{i}=i\Delta x$ and $t^{n}=n\Delta t$ so temperature is $T_{i}^{n}=T(x_{i},t^{n})$. We use forward difference in time

$$\partial_{t}f\Rightarrow\frac{f_{i}^{n+1}-f_{i}^{n}}{\Delta t}$$

(15)

and centered differences in space

$$\partial_{x}f\Rightarrow\frac{f_{i+{\frac{1}{2}}}^{n}-f_{i-{\frac{1}{2}}}^{n}}{\Delta x}\;,\qquad\partial_{x}^{2}f\Rightarrow\frac{f_{i+1}^{n}-2f_{i}^{n}+f_{i-1}^{n}}{\Delta x^{2}}\;.$$

(16)

Note that for the spatial index the integer values are at cell centers while half-integer values are at cell faces.

With this we have the forward Euler (FE) scheme

	$\displaystyle T_{i}^{n+1}=T_{i}^{n}$	$\displaystyle+\frac{\kappa\Delta t}{\Delta x^{2}}\left(T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}\right)$
		$\displaystyle+\frac{\alpha\Delta t}{\Delta x}\left(T_{i+{\frac{1}{2}}}^{n}\tilde{Z}_{i+{\frac{1}{2}}}^{n}-T_{i-{\frac{1}{2}}}^{n}\tilde{Z}_{i-{\frac{1}{2}}}^{n}\right).$		(17)

For the temperature at a face we can use the arithmetic average, $T_{i+{\frac{1}{2}}}^{n}=(T_{i+1}^{n}+T_{i}^{n})/2$. Similarly, the standard Predictor-Corrector (PC) scheme is

	$\displaystyle T_{i}^{*}=T_{i}^{n}$	$\displaystyle+\frac{\kappa\Delta t}{\Delta x^{2}}\left(T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}\right)$
		$\displaystyle+\frac{\alpha\Delta t}{\Delta x}\left(T_{i+{\frac{1}{2}}}^{n}\tilde{Z}_{i+{\frac{1}{2}}}^{n}-T_{i-{\frac{1}{2}}}^{n}\tilde{Z}_{i-{\frac{1}{2}}}^{n}\right)$		(18)

and

	$\displaystyle T_{i}^{n+1}=\frac{1}{2}\Big{[}T_{i}^{n}+T_{i}^{*}$	$\displaystyle+\frac{\kappa\Delta t}{\Delta x^{2}}\left(T_{i+1}^{*}-2T_{i}^{*}+T_{i-1}^{*}\right)$
		$\displaystyle+\frac{\alpha\Delta t}{\Delta x}\left(T_{i+{\frac{1}{2}}}^{*}\tilde{Z}_{i+{\frac{1}{2}}}^{n}-T_{i-{\frac{1}{2}}}^{*}\tilde{Z}_{i-{\frac{1}{2}}}^{n}\right)\Big{]}$		(19)

where $T_{i}^{*}$, computed in the predictor step, is used in the corrector step.²²2Random numbers used in predictor and corrector steps may be the same or different; see [6].

The discretized noise has variance

$$\langle\tilde{Z}_{i+{\frac{1}{2}}}^{n}\tilde{Z}_{j+{\frac{1}{2}}}^{m}\rangle=\frac{1}{\mathcal{A}}~{}\frac{\delta^{\mathrm{Kr}}_{i,j}}{\Delta x}~{}\frac{\delta^{\mathrm{Kr}}_{n,m}}{\Delta t}=\frac{\delta^{\mathrm{Kr}}_{i,j}}{\Delta V}~{}\frac{\delta^{\mathrm{Kr}}_{n,m}}{\Delta t}$$

(20)

where $\Delta V$ is the volume of a grid point cell. Numerically this noise is generated as

$$\tilde{Z}_{i+{\frac{1}{2}}}^{n}=\frac{1}{\sqrt{\Delta V\Delta t}}~{}\mathfrak{N}_{i+{\frac{1}{2}}}^{n}$$

(21)

where $\mathfrak{N}$ are independent, normal (Gaussian) distributed psuedo-random numbers.

The Python program, StochasticHeat, illustrates the numerical methods described in the previous section for the stochastic heat equation (SHE). The program is listed in the appendix and the most up-to-date version can be downloaded as a Jupyter Notebook from github.com/AlejGarcia/IntroFHD. The simulation computes trajectories for the one-dimensional SHE and analyzes their statistical properties (e.g., variance of temperature). The reader is encouraged to download and run the program since results are obtained in about 5 minutes on a laptop.

At the top of the program there are various global options and parameters:

• •

  Periodic or Dirichlet boundary conditions

• •

  Thermodynamic equilibrium or constant temperature gradient

• •

  Initialize with or without temperature perturbations

• •

  Run simulation with or without thermal fluctuations

• •

  Use Forward Euler (FE) or Predictor-Corrector (PC) scheme

• •

  Number of grid cells, $N$

• •

  Number of time steps, $M$

Other parameters (e.g., system length, $\ell$) can also be changed from their default values within the program. Simulation results in this section are from runs using $N=32$ grid cells, running for either $M=2\times 10^{6}$ (equilibrium cases) or $M=2\times 10^{7}$ (temperature gradient case) time steps. Default values were used for the other parameters (see program listing in Appendix).

The program computes time-averaged quantities, for example the average temperature in a cell is

$$\langle T_{i}\rangle=\frac{1}{M_{\mathrm{s}}}~{}\sum_{n=n_{0}+1}^{n_{0}+M_{\mathrm{s}}}T_{i}^{n}$$

(22)

where $M_{\mathrm{s}}$ is the number of statistical samples; the system “relaxes” for the first $n_{0}$ time steps before sampling begins. At thermodynamic equilibrium the variance and static correlation of temperature fluctuations ($\delta T_{i}^{n}=T_{i}^{n}-\langle T_{i}\rangle$) are predicted by statistical mechanics to be [7]

$$\langle\delta T_{i}^{2}\rangle=\frac{k_{B}\langle T_{i}\rangle^{2}}{\rho c_{V}\Delta V}\qquad\mathrm{and}\qquad\langle\delta T_{i}\delta T_{j}\rangle=\langle\delta T_{i}^{2}\rangle~{}\delta^{\mathrm{Kr}}_{i,j}$$

(23)

with the average temperature, $\langle T_{i}\rangle$, equal to the equilibrium temperature, $T_{\mathrm{eq}}$. Figures 1 and 2 show that the simulation results for Dirichlet boundaries (i.e., fixed temperature at the boundaries) are in good agreement with theory, especially for the PC scheme.

Figure 3 shows that temperature correlations are slightly different for periodic boundary conditions. This is because, due to conservation of energy, $\sum_{i}\rho c_{V}\langle\delta T_{i}\rangle=0$, the correlation is shifted by an additive constant such that $\sum_{i}\langle\delta T_{i}\delta T_{j}\rangle=0$. This is a common feature in systems with conserved variables, such as density fluctuations in closed systems.

The assessment of numerical schemes in CFHD is best done by measuring the static structure factor for fluctuations. For the SHE, we define the unnormalized structure factor, $S_{k}=\langle\hat{T}_{k}\hat{T}^{*}_{k}\rangle$, where

$$\hat{T}_{k}=\sum_{j=0}^{N-1}T_{j}\exp(-2\pi\mathrm{i}\,jk/N)$$

(24)

is the discrete Fourier transform with wavenumber index $k=0,1,\ldots,N-1$, and $\hat{T}^{*}_{k}$ is its complex conjugate transpose.³³3In CFHD the structure factor is often normalized to unity using the equilibrium variance. At thermodynamic equilibrium the discretized structure factor is

$$S_{k}=\frac{k_{B}T_{\mathrm{eq}}^{2}}{\rho c_{V}}~{}N.$$

(25)

Figure 4 shows that the FE scheme has significant error for large $k$ while the PC scheme is accurate for all $k$. An analysis of discretization errors for various numerical schemes is presented in [6]; for FE the error in $S_{k}$ is $O(\kappa\Delta t\,k^{2})$ while for PC it is $O(\kappa^{2}\Delta t^{2}k^{4})$.

Finally we consider a non-equilibrium scenario, specifically a linear temperature gradient imposed by Dirichlet boundary conditions. In this case there is a weak, long-range correlation of temperature fluctuations [8]:

$$\langle\delta T_{i}\delta T_{j}\rangle=\frac{k_{B}T_{\mathrm{eq}}^{2}}{\rho c_{V}\Delta V}\delta^{\mathrm{Kr}}_{i,j}+\frac{k_{B}(\nabla T)^{2}}{\rho c_{V}\mathcal{A}\ell}\times\left\{\begin{array}[]{cc}x_{i}(\ell-x_{j})&\mathrm{if}~{}x_{i}<x_{j},\\ x_{j}(\ell-x_{i})&\mathrm{otherwise}.\end{array}\right.$$

(26)

Figure 5 shows that the simulation results are in good agreement with this theoretical result.

In fluctuating hydrodynamics there are several stochastic diffusion equations that are closely related to the stochastic heat equation. For example, the diffusion of species mass is described by a similar stochastic PDE with the main difference being the functional form of the stochastic flux. As above we start with

$$\partial_{t}(\rho c)=-\nabla\cdot\boldsymbol{F}$$

(27)

where $c$ is the species mass fraction and the species flux is

$$\boldsymbol{F}=\overline{\boldsymbol{F}}+\widetilde{\boldsymbol{F}}=\mathcal{L}\boldsymbol{\mathcal{X}}+\widetilde{\boldsymbol{F}}.$$

(28)

Assuming that the system is isothermal (i.e., no Soret effect), irreversible thermodynamics tells us that the thermodynamic force is [4, 5]

$$\boldsymbol{\mathcal{X}}=-\nabla\left(\frac{\mu}{T}\right)=-\frac{\nabla\mu}{T}$$

(29)

where $\mu$ is the chemical potential. For ideal solutions

$$\mu=\mu_{0}(T)+\frac{k_{B}T}{m}\ln c\qquad\mathrm{so}\qquad\nabla\mu=\frac{k_{B}T}{mc}~{}\nabla c$$

(30)

where $m$ is the particle mass. Since the phenomenological law for species flow is Fick’s law, $\overline{\boldsymbol{F}}=-\rho D\nabla c$, then the Onsager coefficient is

$$\mathcal{L}=\frac{m\rho D\,c}{k_{B}}\;.$$

(31)

The noise correlation is

$$\langle\widetilde{\boldsymbol{F}}(\mathbf{r},t)\widetilde{\boldsymbol{F}}(\mathbf{r}^{\prime},t^{\prime})\rangle=2m\rho D\,c~{}\delta(\mathbf{r}-\mathbf{r}^{\prime})~{}\delta(t-t^{\prime})$$

(32)

so the stochastic species diffusion equation is

$$\partial_{t}\,c=\nabla\cdot\left(D\nabla c+\sqrt{\frac{2D\,c}{n_{0}}}~{}\widetilde{\boldsymbol{Z}}\right)$$

(33)

where $n_{0}=\rho/m$ is the number density of the pure state ($c=1$). This may also be written in terms of species number density, $n=c\,n_{0}$, as

$$\partial_{t}\,n=\nabla\cdot(D\nabla n+\sqrt{2D\,n}~{}\widetilde{\boldsymbol{Z}}).$$

(34)

This stochastic PDE is also known as the Dean-Kawasaki equation. Recall that the stochastic heat equation is

$$\partial_{t}\,T=\nabla\cdot\left(\kappa\nabla T+\sqrt{\frac{2\kappa k_{B}\,T^{2}}{\rho c_{V}}}~{}\widetilde{\boldsymbol{Z}}\right).$$

(35)

Note that the deterministic parts of the two equations above are functionally similar while the stochastic parts are distinctly different.

The stochastic species diffusion equation can also be derived by coarse-graining the discrete random walk model for diffusion. Additional examples of discrete models that have associated stochastic diffusion equations are: the excluded random walk model; the “train model” for the diffusion of momentum [9, 10]; and the “Knudsen chain” model for gas effusion [11].

The general form of the fluctuating hydrodynamic equations for a set of state variables $\boldsymbol{U}$ is

$$\partial_{t}\boldsymbol{U}=-\nabla\cdot(\boldsymbol{\mathcal{F}}_{H}(\boldsymbol{U})+\boldsymbol{\mathcal{F}}_{D}(\boldsymbol{U})+\boldsymbol{\mathcal{N}}(\boldsymbol{U})\widetilde{\boldsymbol{Z}})$$

(36)

where $\boldsymbol{\mathcal{F}}_{H}$ and $\boldsymbol{\mathcal{F}}_{D}$ are the hyperbolic and diffusive fluxes with $\boldsymbol{\mathcal{N}}$ being the noise amplitude. For example, for the stochastic species diffusion equation $\boldsymbol{U}=n$, $\boldsymbol{\mathcal{F}}_{H}=0$, $\boldsymbol{\mathcal{F}}_{D}=-D\nabla n$, and $\boldsymbol{\mathcal{N}}=\sqrt{2Dn}$. In general, the fluctuation-dissipation theorem relates $\boldsymbol{\mathcal{F}}_{D}$ to $\boldsymbol{\mathcal{N}}$ and tells us that the stochastic term is independent of the hyperbolic term.

This brief introduction to CFHD has focused on stochastic diffusion equations due to the pedagogical benefit of their simplicity. We close with an outline of other fluctuating hydrodynamic equations including references describing finite volume schemes for solving them.

• •

  Stochastic Burger’s equation [12]:

  $$\partial_{t}\,c=-\partial_{x}\,(ac(1-c)-D\,\partial_{x}c+\sqrt{2Dc(1-c)}\,\tilde{Z}).$$

  (37)

  Note that for $c\ll 1$ this resembles the stochastic species diffusion equation with the addition of a hyperbolic flux, $\boldsymbol{F}_{H}=ac(1-c)$.

• •

  Stochastic “train model” equations [9, 10]:

  	$\displaystyle\partial_{t}\,\rho$	$\displaystyle=-\partial_{x}(-D\partial_{x}\rho+\sqrt{2mD\rho}\,\tilde{Z}),$		(38)
  	$\displaystyle\partial_{t}\,(\rho u)$	$\displaystyle=-\partial_{x}(-D\partial_{x}(\rho u)+\sqrt{2mD\rho u^{2}}\,\tilde{Z}).$		(39)

  These equations are a simple model for the transport of momentum density, $\rho u$, in a fluid.

• •

  Reaction-diffusion systems [13]:
  The reaction-diffusion system has $N_{s}$ species diffusing and undergoing $N_{r}$ reactions. By denoting the number density of species $s$ by $n_{s}$, the equations of fluctuating reaction-diffusion can be written formally as the stochastic PDEs

  $$\frac{\partial n_{s}}{\partial t}=\nabla\cdot\left(D_{s}\nabla n_{s}+\sqrt{2D_{s}n_{s}}\boldsymbol{\mathcal{Z}}_{s}^{(D)}\right)+\sum_{r=1}^{N_{r}}\nu_{sr}\left(a_{r}({\bf n})+\sqrt{a_{r}({\bf n})}\mathcal{Z}_{r}^{(R)}\right)$$

  (40)

  where $D_{s}$ is the diffusion coefficient of species $s$, $a_{r}({\bf n})$ is the propensity function indicating the rate of reaction $r$, and $\nu_{sr}$ is the stoichiometric coefficient of species $s$ in reaction $r$.

• •

  Compressible Navier-Stokes equations (single species fluid) [14, 15]:

  	$\displaystyle\partial_{t}\,\rho=-\nabla\cdot\left(\rho\boldsymbol{u}\right),$		(41)
  	$\displaystyle\partial_{t}\,(\rho\boldsymbol{u})=-\nabla\cdot\left[(\rho\boldsymbol{u}\otimes\boldsymbol{u}+p\mathbb{I})+\overline{\Pi}+\widetilde{\Pi}\right],$		(42)
  	$\displaystyle\partial_{t}\,(\rho E)=-\nabla\cdot\left[\left(\rho E+p\right)\boldsymbol{u}+\overline{\boldsymbol{Q}}+\widetilde{\boldsymbol{Q}}+\left(\overline{\Pi}+\widetilde{\Pi}\right)\cdot\boldsymbol{u}\right],$		(43)

  where $\boldsymbol{u}$, $p$, and $\Pi=\overline{\Pi}+\widetilde{\Pi}$ are fluid velocity, pressure, and stress tensor; total specific energy is $E=e+\frac{1}{2}u^{2}$. Note that (43) reduces to the stochastic heat equation when $\boldsymbol{u}=0$. See references for explicit expression of the deterministic and stochastic fluxes.

• •

  Compressible Navier-Stokes equations (multi-species fluid) [15, 16]:
  Same as above but replace (41) with

  $$\partial_{t}\,\rho_{k}=-\nabla\cdot(\rho_{k}\boldsymbol{u}+\overline{\boldsymbol{F}}_{k}+\widetilde{\boldsymbol{F}}_{k})$$

  (44)

  where $\rho_{k}$ is the mass density for species $k$ and $\rho=\sum_{k}\rho_{k}$. For ideal mixtures the diffusive flux, $\boldsymbol{F}_{k}$, is similar to that of the species diffusion equation described in the previous section.

• •

  Hydrodynamics plus chemical reactions [17, 18]:
  Add source terms to the right-hand side of (44) for the deterministic and stochastic rates of reactions, $\overline{\Omega}_{k}$ and $\widetilde{\Omega}_{k}$. These can be formulated from the chemical Langevin equation [19].

• •

  Incompressible / low Mach hydrodynamic equations [20, 21, 22, 23]:
  The multi-species methodology can be extended to model isothermal mixtures of miscible incompressible liquids. Incompressibility of the fluids leads to a constrained evolution, similar to the incompressible Navier-Stokes equations, given by

  	$\displaystyle\partial_{t}\,(\rho\boldsymbol{u})=-\nabla\cdot\left[\rho\boldsymbol{u}\otimes\boldsymbol{u}+\overline{\Pi}+\widetilde{\Pi}\right]-\nabla\pi,$		(45)
  	$\displaystyle\partial_{t}\,\rho_{k}=-\nabla\cdot(\rho_{k}\boldsymbol{u}+\overline{\boldsymbol{F}}_{k}+\widetilde{\boldsymbol{F}}_{k}),$		(46)
  	$\displaystyle\nabla\cdot\boldsymbol{u}=-\nabla\cdot\left(\sum_{k}\frac{\overline{\boldsymbol{F}}_{k}+\widetilde{\boldsymbol{F}}_{k}}{\bar{\rho}_{k}}\right).$		(47)

  Here, $\pi$ is a perturbational pressure and $\bar{\rho}_{k}$ is the pure component density for species $k$. The resulting system retains in influence of fluctuations on mixing but is considerable more computationally efficient, particularly for liquids with a high sound speed.

• •

  Stochastic Poisson-Nernst-Planck equations [24, 25]:
  The incompressible version of CFHD for multi-species fluids can model electrolyte solutions by including charged species (ions). The chemical potential becomes the electrochemical potential with the electrical mobility given by the Nernst–Einstein relation and the diffusive flux by the Nernst-Planck equation. For scales comparable to the Debye length, the electric field is obtained by solving the Poisson equation. At scales larger than the Debye length, the fluid is electro-neutral. Electro-neutrality can be imposed as a constraint by solving a modified elliptic equation to compute the electric potential.

• •

  Multi-phase fluids [26, 27, 28, 29]:
  By adding non-local gradient terms to the free energy, one can incorporate interfacial tension in CFHD model to treat multiphase systems. This can be done for both single component systems (e.g., liquid/vapor) and for multicomponent systems (e.g., oil/water). For multicomponent systems, the low Mach number version of the methodology has also be derived. These types of systems introduce additional higher order operators into the momentum equations and, for multicomponent systems, the species transport equations. See the references for specific forms of these equations.