Fluid Motion from the Boltzmann Equation
November 2024
At the macroscopic scale, fluid motion is well described by the Navier-Stokes equations . A set of partial differential equations describing the fluid as a continuum. This is quite a remarkable fact given that fluids ultimately consist of discrete atoms (or molecules) governed by the laws of Newtonian mechanics. It is even possible to mathematically show that the equations of fluid motion emerge from the microscopic motions of the constituent molecules. I found this fact remarkable when I first learned it years ago.
Between the microscopic dynamics of individual atoms & molecules, and the macroscopic fluid dynamics equations sits the so-called ‘mesoscale’, where we partially summarize the microscopic dynamics and consider the distribution of particle speeds in our fluid. This ‘distribution function’ is governed by the Boltzmann equation, which is our jumping off point.
Summarizing the Boltzmann equation to the macroscale leads directly to the fully general equations of fluid motion; the Cauchy momentum equation . However, some of the terms will involve the full distribution function, tying us to the mesoscopic description. The Chapman-Enskog expansion , by decomposing the distribution function, can then be used to move fully to the macroscopic description.
Contents
Notation
Before diving in let’s define some notation that will make our lives much easier:
$$\langle A \rangle := \int_{-\infty}^\infty A(\mathbf q) f(\mathbf q) \,d^3 q$$Directly from its definition as an integral it’s clear that this obeys certain properties we can use in our derivations:
$$\begin{align} \langle \lambda A + B \rangle &= \lambda \langle A \rangle + \langle B \rangle & \text{Linearity} \\ \int A \partial_x f \,d^3 q &= \partial_x \langle A \rangle - \langle \partial_x A \rangle &\text{Product Rule} \\ \end{align}$$where $x$ is any variable that $A$ and $f$ depend on except $q$. We also care about derivatives wrt $q$, but we cannot pull these outside the brackets (integral). In this case we get a slightly different result, namely:
$$\int A \partial_q f \,d^3 q = - \langle \partial_q A \rangle $$Proof
$$ \int A \partial_q f\,d^3 q = \int \partial_q (A f) - f \partial_q A \,d^3 q = -\langle \partial_q A \rangle$$The first integral term is zero because of the fundamental theorem of Calculus , which tells us
$$\int_a^b \partial_q (f(q) A(q)) \, dq = f(b)A(b) - f(a)A(a)$$and the fact that $f$ is a probability distribution, so we have
$$\lim_{|q| \rightarrow \infty} f(q) = 0$$Thus
$$\int_{-\infty}^{\infty} \partial_{q} (f(q) A(q)) \, dq = 0$$This assumes that $f$ goes to zero fast enough that the product $fA$ approaches zero at infinity. This is certainly the case at equilibrium where our $f$ is Gaussian (and so decays exponentially) and our $A$’s are polynomials.
Moments of Boltzmann Equation
The Boltzmann equation is
$$\Omega f = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial \mathbf x}\cdot\mathbf q + \frac{\partial f}{\partial \mathbf q}\cdot \mathbf g$$where \(\mathbf g=\mathbf g(t, \mathbf x)\) is a function of space and time, but not (microscopic) velocity $\mathbf q$. Clearly what we observe at the macroscopic scale emerges from the net configuration and motions of the constituent molecules. That is to say the hydrodynamic variables ($\rho$, $\mathbf u$, etc) can be calculated as averages over the microscopic velocity $\mathbf q$. Therefore taking the same moments of the Boltzmann equation can give us the equations of motion for the macroscopic quantities.
Continuity Equation
To start, we take the zeroth moment of the Boltzmann equation over velocity.
$$\begin{align} \int \Omega \,d\mathbf q &= \int \frac{\partial f}{\partial t} + \frac{\partial f}{\partial \mathbf x} \cdot \mathbf q + \frac{\partial f}{\partial \mathbf q} \cdot \mathbf g \,\, d \mathbf q \\ 0 &= \frac{\partial}{\partial t} \int f d \mathbf q + \frac{\partial}{\partial \mathbf x}\cdot \int f \mathbf q\, d \mathbf q + \Big(\int \frac{\partial f}{\partial \mathbf q} \,d\mathbf q \Big)\cdot\mathbf g \\ &= \frac{\partial \rho}{\partial t} + \nabla\cdot (\rho \mathbf u) \end{align}$$Where we have used the conservation of mass for the collision operator, and some simple rearrangement of the integrals & derivatives. The last term is less obvious but it follows simply from the properties stated above
$$\int \partial_{\mathbf q} f \,d\mathbf q = -\langle \partial_q 1 \rangle = 0$$Thus the final equation we arrive at is
$$\frac{\partial \rho}{\partial t} = - \nabla \cdot (\rho \mathbf u)$$which is the continuity equation for compressible fluid flow in its Eulerian form. The interpretation of this equation is clear, at a given location the change in density is caused by the net imbalance of mass flow in to or out of that point. If we expand the divergence using the chain-rule, we can convert this to the lagrangian form
$$\boxed{ \frac{\mathrm d \rho}{\mathrm d t} = -\rho \nabla \cdot \mathbf u }$$In this particular case I find the Eulerian form slightly more interpretable.
We assumed nothing about the form of $f$ here, and it does not appear in the final equation anywhere. It is fully in terms of the known macroscopic moments $\rho$ and $\mathbf u$. This means that this equation will hold for any situation we consider, whether equilibrium or non-equilibrium.
Cauchy Momentum Equation
Now let’s take the first moment of the Boltzmann equation wrt micro-velocity. This is slightly clearer in component notation, which we will use for the following.
$$\int q_\alpha \Omega f \,d^3 q = \int q_\alpha \partial_t f \,d^3 q + \int q_\alpha q_\beta \partial_{x_\beta} f \,d^3 q + \int q_\alpha g_\beta \partial_{q_\beta} f \,d^3 q $$The procedure is very similar to the zeroth moment. Start by noticing the left-hand side is zero due to conservation of momentum, then consider the terms on the right-hand side. To avoid symbol overload lets calculate each term in sequence
Term 1
$$\begin{align} \int q_\alpha \partial_t f \,d^3 q &= \partial_t \langle q_\alpha \rangle - \cancelto{0}{\langle \partial_t q_\alpha \rangle}= \partial_t (\rho u_\alpha) \end{align}$$Term 2
$$\begin{align} \int q_\alpha q_\beta \partial_{x_\beta} f \,d^3 q = \partial_{x_\beta} \langle q_\alpha q_\beta \rangle := \partial_{x_\beta} \Pi_{\alpha\beta} \end{align}$$Term 3
$$\begin{align} g_\beta \int q_\alpha \partial_{q_\beta} f \,d^3 q &= - \langle \partial_{q_\beta} q_\alpha \rangle g_\beta \\ &= - \langle \delta_{\alpha\beta} \rangle g_\beta \\ &= - \rho \delta_{\alpha\beta} g_\beta \\ &= -\rho g_\alpha \end{align}$$Where we have used the identity $\partial q_\beta / \partial q_\alpha \equiv \delta_{\alpha\beta}$.
Putting all three together gives finally
$$\frac{\partial}{\partial t} (\rho u_\alpha) = - \frac{\partial \Pi_{\alpha\beta}}{\partial x_\beta} + \rho g_\alpha $$or in vector notation
$$ \frac{\partial}{\partial t}(\rho \mathbf u) = - \nabla\cdot \mathbf\Pi + \rho \mathbf g \tag{1}\label{eq:cauchy_momentum} $$This is the Cauchy momentum equation , which describes momentum transport in fluid flows, in its conservation form. The second moment $\Pi_{\alpha\beta}$ is the momentum density flux tensor. We can decompose it as
$$\mathbf \Pi = \rho\mathbf u \mathbf u - \boldsymbol{\sigma}$$where $\boldsymbol{\sigma}$ is the (Cauchy) stress tensor .
This comes from considering the moment in terms of the relative velocity $\mathbf v = \mathbf q - \mathbf u$;
$$\begin{align} q_\alpha q_\beta &= (u_\alpha + v_\alpha)(u_\beta + v_\beta) \\ &= u_\alpha u_\beta + u_\alpha v_\beta + v_\alpha u_\beta + v_\alpha v_\beta \end{align}$$thus the integral in the second term becomes
$$\begin{align} \Pi_{\alpha\beta} = \langle q_\alpha q_\beta \rangle = \langle u_\alpha u_\beta + u_\alpha v_\beta + u_\beta v_\alpha + v_\alpha v_\beta \rangle \\ = u_\alpha u_\beta\langle 1 \rangle + u_\alpha \langle v_\beta \rangle + \langle v_\alpha \rangle u_\beta + \langle v_\alpha v_\beta\rangle \\ := \rho u_\alpha u_\beta - \sigma_{\alpha\beta} \end{align}$$where the first moment of $f$ wrt $\mathbf v$ is zero by construction, and we have defined
$$ \sigma_{\alpha\beta} := - \int v_\alpha v_\beta f\,d^3 q $$This gives us
$$ \frac{\partial (\rho \mathbf u)}{\partial t} + \nabla\cdot(\rho \mathbf u \mathbf u) = \nabla \cdot \boldsymbol{\sigma} + \rho \mathbf g $$We can further transform the left-hand side using the product-rule to pull the density outside the differentials as
$$\frac{\partial (\rho\mathbf u)}{\partial t} + \nabla\cdot(\rho \mathbf u \mathbf u) = \rho\frac{\partial \mathbf u}{\partial t} + \rho(\mathbf u\cdot \mathbf \nabla)\mathbf u = \rho \frac{\mathrm{d} \mathbf u}{\mathrm{d} t} $$So the left-hand side is the density multiplied by the total (aka material) derivative of the velocity.
This lets us rearrange slightly to get the Lagragian form of the Cauchy momentum equation
$$\boxed{ \frac{\mathrm d \mathbf u}{\mathrm d t} = \frac{\nabla \cdot \boldsymbol{\sigma}}{\rho} + \mathbf g }$$Velocity Euler Equation
The Cauchy momentum equation is fully generic. Unfortunately, though, in general we don’t know the full distribution $f$, and thus neither do we know $\boldsymbol{\sigma}$. We can try to close the equations via constituative relations. However lets assume we are at local equilibrium so that $f = f^\mathrm{eq}$. In this case can find a simple solution.
The definition of $\boldsymbol{\sigma}$ is equivalent to minus the covariance matrix of $f$ times a factor $\rho$. The equilibrium distribution is an isotropic multivariate Gaussian, whose (diagonal) covariance matrix is $(k_B T/m)\mathbf I$ - see link. Thus we get
$$\boldsymbol{\sigma} = -\rho \frac{k_B T}{m} \mathbf I$$We can relate this to the pressure using the ideal gas law ;
$$P=n k_B T = \rho \frac{k_B T}{m}$$Where $n=\rho/m$ is the number density of fluid particles. So the stress tensor at local equilibrium is $\sigma_{\alpha\beta} = - P \delta_{\alpha\beta}$ thus
$$ \nabla \cdot \boldsymbol{\sigma} = \partial_{x_\beta} (-\delta_{\alpha\beta} P) = -\partial_{x_\alpha} P = -\nabla P $$giving, finally, the Euler equation for velocity
$$\boxed{ \frac{\mathrm d \mathbf u}{\mathrm d t} = \frac{-\nabla P}{\rho} + \mathbf g }\tag{2}\label{eq:euler_velocity}$$The Euler equation describes the time evolution of the velocity field for an inviscid fluid.
This equation, however, involves the pressure, which in turn is related to both the density and the temperature via the equation of state $P \propto \rho T$ as above. We must go deeper …
Energy Euler Equation
The temperature is related to the internal energy by the so-called equipartion theorem . The internal energy is given by the central second moment of the distribution function, so we take the second moment of the Boltzmann equation:
$$\int v^2 \Omega f \,d^3 q = \int v^2 \partial_t f \,d^3 q + \int v^2 q_\alpha \partial_{x_\alpha} f \,d^3 q + \int v^2 g_\alpha \partial_{q_\alpha} f \,d^3 q$$where $v^2 = |\mathbf v|^2 = v_\beta v_\beta$. Note we are using $\mathbf v$ here, not $\mathbf q$ as above, because we care about the internal energy. Again the left-hand side is zero, this time due to conservation of energy. Now let’s proceed as above and calculate each term in the same manner as for the first moment:
Term 3
Let’s start with the third term because it ends up being zero so we can forget about it hereafter.
$$\begin{align} \int g_\alpha v_\beta v_\beta \partial_{q_\alpha} f \,d^3 q &= - g_\alpha \langle \partial_{q_\alpha} v_\beta v_\beta \rangle \\ &= - g_\alpha \langle \partial_{v_\beta} (v_\beta v_\beta) \partial_{q_\alpha} v_\beta \rangle \\ &= 2 g_\alpha \delta_{\alpha\beta} \langle v_\beta \rangle \\ &= 0 \end{align}$$where again $\partial_{q_\alpha} v_\beta = \delta_{\alpha\beta}$ has been used, in combination with the chain rule, to transform the the right-hand side.
It should not be unexpected that this term is zero; the acceleration $\mathbf g$ acts to move the mean velocity, but it changes all $\mathbf q$ equally and so the second central moment, i.e. the internal-energy, does not change.
Term 1
Here we need to be a bit careful with the $\partial_t$ because although $\mathbf q$ is not a function of time, $\mathbf v = \mathbf q - \mathbf u$ is since it depends on $\mathbf u$.
$$\begin{align} \int v^2 \partial_t f \,d^3 q &= \partial_t\langle v_\alpha v_\alpha \rangle - \cancelto{0}{\langle \partial_t(v_\alpha v_\alpha)\rangle} \\ &= \frac{\partial(2\rho U)}{\partial t} \end{align}$$The last term equals zero because of the chain-rule and $\langle\mathbf v\rangle=0$. The chain rule gives
$$\langle \partial_{t}(v_\alpha v_\alpha) \rangle = \langle\partial_{v_\alpha}(v_\alpha v_\alpha)\partial_{u_\alpha} v_\alpha \partial_{t} u_\alpha \rangle = \langle 2v_\alpha \cdot (-1) \cdot \partial_t u_\alpha\rangle = -2 \partial_t u_\alpha \langle v_\alpha \rangle = 0$$Note we wrote this with $\partial_t$ but it would read identically for $\partial_{x_\beta}$. This is because the dependence on both $t$ and $\mathbf x$ enter only through $\mathbf v$’s dependence on $\mathbf u$. We will use this for the next term.
Term 2
This one is simple enough but quite verbose. Here, again, we need to be cognisant that $\mathbf v$ is a function of $\mathbf x$ via $\mathbf u$.
$$\begin{align} \int v_\beta v_\beta q_\alpha \partial_{x_\alpha} f \,d^3 q &= \partial_{x_\alpha} \langle q_\alpha v_\beta v_\beta \rangle - \langle \partial_{x_\alpha} (q_\alpha v_\beta v_\beta) \rangle \\ &= \partial_{x_\alpha} \langle u_\alpha v_\beta v_\beta + v_\alpha v_\beta v_\beta \rangle - \langle q_\alpha \partial_{x_\alpha} (v_\beta v_\beta)\rangle \\ &= \partial_{x_\alpha} (u_\alpha 2 \rho U) + \partial_{x_\alpha} \langle v_\alpha v_\beta v_\beta \rangle + \langle q_\alpha 2 v_\beta \partial_{x_\alpha} u_\beta \rangle \\ &= 2 \partial_{x_\alpha} (u_\alpha \rho U) + 2 \partial_{x_\alpha} j_\alpha + 2 \langle q_\alpha v_\beta \rangle \partial_{x_\alpha} u_\beta \\ &= 2 \partial_{x_\alpha} (u_\alpha \rho U) + 2 \partial_{x_\alpha} j_\alpha + 2 \sigma_{\alpha\beta} \partial_{x_\alpha} u_\beta \\ \end{align}$$Where we have defined the heat-flux vector $\mathbf j$ as
$$j_\alpha := \frac{1}{2}\langle v_\alpha v_\beta v_\beta \rangle$$and the last term is simply from
$$\langle q_\alpha v_\beta \rangle = \langle u_\alpha v_\beta + v_\alpha v_\beta \rangle = u_\alpha \langle v_\alpha \rangle + \langle v_\alpha v_\beta \rangle =\langle v_\alpha v_\beta \rangle =: \sigma_{\alpha\beta}$$Combining all this gives us
$$\frac{\partial (\rho U)}{\partial t} + \nabla\cdot(\rho U \mathbf u) = - \nabla \cdot \mathbf j + \boldsymbol\sigma : \nabla\mathbf u$$where the colon indicates the double-dot product (aka the Frobenius product). Again, we can pull the density outside the derivatives with the product-rule & continuity equation
$$ \frac{\partial (\rho U)}{\partial t} + \nabla\cdot(\rho U \mathbf u) = \rho\frac{\partial U}{\partial t} + \rho\mathbf u\cdot\nabla U = \rho\frac{\mathrm dU}{\mathrm dt} $$The proof is exactly analagous to that above for the Cauchy momentum equation, so I won’t repeat it. This gives, finally
$$\boxed{ \rho\frac{\mathrm d U}{\mathrm d t} = - \nabla \cdot \mathbf j + \boldsymbol\sigma : \nabla\mathbf u }$$This looks quite daunting but let’s see what happens at equilibrium. Firstly we can notice that $\mathbf j$ has to be zero at equilibrium. Why? Well, because at equlibrium all odd moments of $\mathbf v$ are zero, not just the first moment.
Proof
Firstly, note that the equilibrium distribution can be factored into a product of 1-dimensional distributions (because it’s isotropic) $f^\mathrm{eq}(\mathbf q) = f^\mathrm{eq}_\mathrm{1D}(q_x)f^\mathrm{eq}_\mathrm{1D}(q_y)f^\mathrm{eq}_\mathrm{1D}(q_z)$ which lets us reorder the iterated integrals over the various components of $q$.
Now, consider the case $\alpha\neq\beta$.
$$\langle v_\alpha v_\beta v_\beta \rangle = \langle v_\alpha \rangle _\mathrm{1D} \langle v_\beta v_\beta \rangle _\mathrm{1D} \langle 1 \rangle _\mathrm{1D}^{d-2}= 0$$so $\langle v_\alpha\rangle=0$ collapses the whole integral to zero.
The other case is $\alpha=\beta$, where we can also use the factorizability to give
$$\langle v_\alpha^3\rangle = \langle v_\alpha^3\rangle _\mathrm{1D} \langle 1 \rangle _\mathrm{1D}^{d-1} = 0$$This equals zero because $f^\mathrm{eq}_\mathrm{1D}$ is a Gaussian and all odd central moments of a Gaussian are zero.
Secondly, we already know what $\boldsymbol\sigma$ is at equilibrium from above. Plugging its value in gives us
$$\sigma_{\alpha\beta} \partial_{x_\alpha} u_\beta = -P \delta_{\alpha\beta} \partial_{x_\alpha}u_\beta = -P \nabla\cdot\mathbf u$$So, at equilbrium, we get
$$\boxed{ \frac{\mathrm d U}{\mathrm d t}= - \frac{P}{\rho} \nabla\cdot\mathbf u }$$Temperature Equation
To close the loop on this section, we can use the aforementioned equipartition theorem to convert $U$ to $T$ to let us calculate $P$, or just relate $U$ to $P$ directly.
Units
Lets stop quickly to think about units.
The equipartition theorem tells us that the temperature is related to the average internal energy of the particles in our fluid. From its definition we see the combined quantity $\rho U$ is the total internal energy per unit volume. We know also that $n=\rho/m$ is the number density, that is the particle count per unit volume. Therefore the average energy per particle needed for the equipartition theorem is simply $(\rho U)/n = mU$. For completeness this means $U$ is the internal energy per unit mass.
For a 3-dimensional monatomic gas we have $mU = \frac{3}{2}k_B T$ and $P = \rho k_B T / m$. Converting the above equation to be in terms of $T$ gives
$$ \frac{\mathrm d T}{\mathrm d t} = \frac{2}{3}\frac{m}{k_B\rho}\Big(\boldsymbol\sigma:\nabla\mathbf u - \nabla\cdot\mathbf j\Big) $$which at equilibrium is
$$ \frac{\mathrm d T}{\mathrm d t} = - \frac{2}{3} T \nabla\cdot\mathbf u $$This is the Euler equation for the temperature.
Chapman-Enskog Expansion
We saw in the last section that the moments of the Boltzmann equation give us the macroscopic equations of fluid motion. The most general form of which being the Cauchy momentum equation \eqref{eq:cauchy_momentum}. However, we were only able to say anything about the macroscopic velocity (and higher moments) at local equilibrium, which gave us the Euler equation \eqref{eq:euler_velocity}.
To move beyond the equilibrium case we can use the so-called Chapman-Enskog expansion . The core idea of the expansion is to assume that we are not too far from equilibrium, and so we can perform a perturbative expansion around equilibrium.
When is this a good approximation? If the mean free path ($\lambda$) of a particle is very long compared to the length scale of interest ($L$), then it’s possible for the system to be very far from equilibrium. To see why consider any small volume element, it is unlikely that any particle will have collided inside the volume if $\lambda >\!\!> L$. Conversely, if $\lambda <\!\!< L$ then we can expect the the particles to have collided many times within the volume, and thus the distribution of velocities to be at least approximately at equilibrium.
The ratio of the mean free path $\lambda$ to the characteristic length-scale $L$ is a dimensionless number called the Knudsen number
$$\mathrm{Kn} = \frac{\lambda}{L}$$The approximation is thus valid for small Knudsen numbers.
In fact, the Chapman-Enskog expansion decomposes $f$ into terms with varying orders of $\mathrm{Kn}$. Defining $f^\mathrm{eq}:=f^{(0)}$ we can write the expansion
$$f = f^{(0)} + \varepsilon f^{(1)} + \varepsilon^2 f^{(2)} + \dots$$where $\varepsilon^n$ labels terms which are of the same order in $\mathrm{Kn}$.
Here we are using $\varepsilon$ purely to label terms which are of the same order in $\mathrm{Kn}$. It is a book-keeping device, not something with a numerical value itself. The actual numerical value of $f$ would be given by $f = \sum_i f^{(i)}$, or equivalently by simply setting $\varepsilon=1$ at the end.
The benefit of having $\varepsilon$ is that it allows us to track the ‘size’ of terms as we manipulate the expansion. For example we can see immediately that $(\varepsilon^n f^{(n)})(\varepsilon^m f^{(m)}) = \varepsilon^{n+m}f^{(n)}f^{(m)}$ is $O(\mathrm{Kn}^{n+m})$.
A vague analogy would be how the dual number system is written as $a + \varepsilon b$, except in the Chapman-Enskog expansion we are just dealing with a decomposition of a regular real number into a sum of other real numbers.
An exactly equivalent alternative would be to say that $\varepsilon$ actually has a numerical value of $\mathrm{Kn}$ but that the terms $f^{(n)}$ are correspondingly larger by a factor $\varepsilon^n$ than in the case where it is purely a label.
The second cruicial ingredient of the Chapman-Enskog procedure is a functional Ansatz that $f$ depends on time only through the macroscopic hydrodynamic fields. This means that any temporal dependence of $f$ enters only indirectly through the dependence of the hydrodynamic variables on time; i.e
$$f(t, \mathbf x, \mathbf q) = f(\mathbf x, \mathbf q, \rho(t), \mathbf u(t), T(t))$$This makes physical sense; ultimately we seek an equation for the time dependence of (the moments of) $f$ in terms of the state of the hydrodynamic variables at the current time (this state can obviously include various spatial gradients of the hydrodynamic fields).
Since we are assuming the collision term is dominates the transport terms we write the Boltzmann equation with an extra factor which indicates the transport terms are of order the Knudsen number:
$$\frac{1}{\varepsilon}\Omega f = \frac{\partial f}{\partial t} + \mathbf q \cdot \frac{\partial f}{\partial \mathbf x}$$The additional factor passes at least a cursory sniff-test; if $\varepsilon$ is small then $\varepsilon^{-1}$ will be large and the collsion term will dominate. But how can we just change the Boltzmann equation by introducing a new factor? Well, in truth, we aren’t. Remember that $\varepsilon$ is a label, it has no numerical value. So what we are doing is saying we will only consider situations where it is the case that the transport terms are \(O(\varepsilon)\) compared to the collision term being \(O(1)\). This defines the regimes in which the Chapman-Enskog procedure is valid.
Collision Operator
The collision operator is in general non-linear, which makes the expansion less useful. However, Boltzmann’s original ‘molecular chaos’ collision operator can be expressed as a symmetric bilinear operator. The molecular chaos collision operator reads
$$\Omega(f)(\mathbf q) = \iint h\, \sigma(h, S) \big[ f(\mathbf q^\prime)f(\mathbf q_1^\prime) - f(\mathbf q)f(\mathbf q_1)\big]\,dS \,d^3\mathbf q_1 $$where the primed velocities are the velocities resulting from the collision of two particles with velocities $\mathbf q$ and $\mathbf q_1$, $h = |\mathbf q - \mathbf q_1| = |\mathbf q^\prime - \mathbf q_1^\prime|$, $\sigma$ is the differential cross-section for the collision and $S$ is a solid-angle.
The BGK collision operator commonly used in lattice Boltzmann simulations just is linear (well, affine).
$$\Omega_\mathrm{BKG} (f) = \frac{1}{\tau}(f^\mathrm{eq} - f)$$Thus the analysis would be a little different in that case.
Bilinear Form
We want to express this as $\Omega(f) = J(f, f)$ where $J(f, g)$ is bilinear and symmetric (i.e. $J(f, g) = J(g, f)$). At first glance we can write
$$B(f, g)(\mathbf q) = \iint h\, \sigma(h, S) \big[ f(\mathbf q^\prime)g(\mathbf q_1^\prime) - f(\mathbf q)g(\mathbf q_1)\big]\,dS \,d^3\mathbf q_1 $$this is certainly bilinear and obeys $B(f, f) = \Omega(f)$. However, it is not symmetric under interchange of $f$ and $g$. That’s no trouble though, lets write $J(f, g) = \frac{1}{2}(B(f, g) + B(g, f))$. $J$ is symmetric and clearly retains the properties of bilinearity and $J(f, f) = \Omega(f)$ Writing it out in full we get the symmetric bilinear form of the collision operator as
$$\begin{align} J(f, g)(\mathbf q) = \frac{1}{2}\iint h\, \sigma(h, S) \big[ &f(\mathbf q^\prime)g(\mathbf q_1^\prime) + g(\mathbf q^\prime)f(\mathbf q_1^\prime) \\ &\qquad - f(\mathbf q)g(\mathbf q_1) - g(\mathbf q)f(\mathbf q_1)\big]\,dS \,d^3\mathbf q_1 \end{align}$$Linearized Collision Operator
We can go a step further and linearize $J$ for small deviations away from $f^{(0)}$, which will be useful later. Take the Chapman-Enskog expansion of $f$ and pull out a common factor of $f^{(0)}$:
$$f = f^{(0)}(1 + \varepsilon \phi + \dots) $$where we have defined $f^{(1)} = f^{(0)} \phi$. Now plug the expansion into $J$ and keep only terms up to first order in $\varepsilon$
$$\begin{align} J(f, f) &= J(f^{(0)} + \varepsilon f^{(0)} \phi + \dots,\, f^{(0)} + \varepsilon f^{(0)} \phi + \dots) \\ &= J(f^{(0)}, f^{(0)}) + 2\varepsilon J(f^{(0)}, f^{(0)}\phi) + O(\varepsilon^2)\\ \end{align}$$We know that $J(f^{(0)}, f^{(0)})=0$ since at equilibrium the collision operator does not change $f$, and the factor of two comes from the symmetry of $J$ so $J(f^{(0)}, f^{(1)}) = J(f^{(1)}, f^{(0)})$.
$$\begin{align} 2J(f^{(0)}, f^{(1)}) &= \iint h \sigma(\sigma, S) \big[{f^{(0)}}^\prime{f^{(0)}}^\prime_1 \phi^\prime_1 + {f^{(0)}}^\prime_1{f^{(0)}}^\prime \phi^\prime \\ &\qquad\qquad\qquad\qquad - {f^{(0)}}{f^{(0)}}_1 \phi_1 - {f^{(0)}}_1 {f^{(0)}} \phi \big]\,dS \,d^3\mathbf q_1 \\ &= \iint h \sigma(\sigma, S) {f^{(0)}}{f^{(0)}}_1 \big[\phi^\prime_1 + \phi^\prime - \phi_1 - \phi \big]\,dS \,d^3\mathbf q_1 \end{align}$$The last step is possible because ${f^{(0)}}{f^{(0)}}_1 = {f^{(0)}}^\prime {f^{(0)}}^\prime_1$.
Now, we can pull $f^{(0)}$ outside the integral (because it is not a function of the integration variable) then what’s left in the integral is the linearized collision operator: $$C(\phi) = \iint h \sigma(\sigma, S) {f^{(0)}}_1 \big[\phi^\prime_1 + \phi^\prime - \phi_1 - \phi \big]\,dS \,d^3\mathbf q_1$$So we have
$$2J(f^{(0)}, f^{(1)}) = f^{(0)} C(\phi)$$Using the Expansion
Now lets plug this expansion directly into the Boltzmann equation. For this, we can neglect the force-density term (by setting $\mathbf g = \mathbf 0$) because the resulting term in the Cauchy momentum term is exact for all $f$.
LHS
$$\begin{align} \Omega f &= J(f, f) \\ &=J\Big(\sum_{i=0}^\infty \varepsilon^i f^{(i)}, \sum_{j=0}^\infty \varepsilon^j f^{(j)}\Big) \\ &= \sum_{i=0}^\infty \sum_{j=0}^\infty \varepsilon^{i+j}J(f^{(i)}, f^{(j)}) \\ &= \sum_{n=0}^\infty \varepsilon^n \sum_{m=0}^n J(f^{(m)}, f^{(n-m)}) \end{align}$$In the last step we have collected all the terms where $i+j=n$ together.
RHS
$$\begin{align} \varepsilon \partial_t f + \varepsilon (\partial_{x_\alpha} f) q_\alpha &= \varepsilon \partial_t \Big(\sum_{n} \varepsilon^n f^{(n)}\Big) + \varepsilon \partial_{x_\alpha} \Big(\sum_{n} \varepsilon^n f^{(n)}\Big) q_\alpha \\ &= \sum_{n} \varepsilon^{n+1} \big(\partial_t f^{(n)} + \partial_{x_\alpha} f^{(n)} q_\alpha\big) \end{align}$$Since we have made the functional Ansatz for $f$ we can expand the time derivatives using the chain rule
$$\frac{\partial f^{(n)}}{\partial t} = \frac{\partial \rho}{\partial t}\frac{\partial f^{(n)}}{\partial \rho} + \frac{\partial \mathbf u}{\partial t}\cdot\frac{\partial f^{(n)}}{\partial \mathbf u} + \frac{\partial T}{\partial t}\frac{\partial f^{(n)}}{\partial T}$$Here we must take care with the time derivatives of the hydrodynamic variables. They also involve terms at different orders in $\varepsilon$, so it’s important for us to track these. The time derivatives are given by the continuity, momentum, and temperature equations we derived previously. These equations involve higher moments of $f$ and thus have, in general, non-zero terms involving \(f^{(n)}\) for \(n\geq 1\).
As with \(f\), lets denote the $O(\varepsilon^n)$ part of any quantity with a superscript \(\cdot^{(n)}\). So we have
$$\frac{\partial \rho}{\partial t} = \sum_n \varepsilon^n (\partial_t \rho)^{(n)}$$and similarly for \(\mathbf u\) and \(T\).
In the literature it is common to decompose the time derivative operator itself by writing
$$\frac{\partial}{\partial t} := \sum_n \varepsilon^n \partial_t^{(n)}$$This can seem a bit magical at first, but don’t think of $\partial^{(n)}_t$ as being a deriative itself, rather it simply picks out the $O(\varepsilon^n)$ part of the overall time derivative. This is clear with the notation above
$$\frac{\partial A}{\partial t} = \sum_n \varepsilon^n \partial_t^{(n)} A \equiv \sum_n \varepsilon^n (\partial_t A)^{(n)}$$Since the operator decomposition notation is prevalent in the literature I’ll use it henceforth.
Now, here things get a bit subtle. The local equilibrium distribution is determined entirely by the hydrodynamic variables $\rho$, $\mathbf u$ and $T$. Moreover, by construction, the equilibrium distribution with these moments is what we are expanding about. This means that we must have $\rho = \rho^{(0)}$, $\mathbf u = \mathbf u^{(0)}$ and $T = T^{(0)}$. The same, however, is not true of the higher moments $\mathbf j$ and $\boldsymbol{\sigma}$, for which we must consider the full expansion.
Having made this observation, let’s apply the expanded time derivative notation to the hydrodynamic equations from above, giving
$$\begin{align} \sum_n \varepsilon^n \partial_t^{(n)} \rho &= -\nabla (\rho \mathbf u) \\ \sum_n \varepsilon^n \partial_t^{(n)} \mathbf u &= \frac{1}{\rho}\nabla \cdot\boldsymbol\sigma - (\mathbf u\cdot\nabla)\mathbf u \\ \sum_n \varepsilon^n \partial_t^{(n)} T &= \sum_n \varepsilon^n \Big(\boldsymbol\sigma^{(n)}:\nabla \mathbf u - \nabla\cdot\mathbf j^{(n)}\Big) - \mathbf u\cdot \nabla T \end{align}$$From their definitions as moments of $f$ we know $\mathbf j^{(n)}$ and $\boldsymbol\sigma^{(n)}$ involve only $f^{(n)}$:
$$\begin{align} \mathbf j^{(n)} &= \Big(\int \mathbf v\,v^2 f \,d^3 q\Big)^{(n)} = \Big(\int \mathbf v\,v^2 \sum_m \varepsilon^m f^{(m)} \,d^3 q\Big)^{(n)} \\ &= \Big(\sum_m \varepsilon^m \int \mathbf v\,v^2 f^{(m)} \,d^3 q\Big)^{(n)} =\int \mathbf v\,v^2 f^{(n)} \,d^3 q \end{align}$$The zero order terms give us back exactly the Euler equations. This should be clear because they are exactly the equations we get by assuming we are at equilibrium, which is equivalent to stating \(f = f^{(0)}\).
We will eventually expand the spatial derivative of $f$ in terms of the hydrodynamic field gradients. Note however, that we only need to track orders of \(\varepsilon\) for the time derivatives, not the spatial. This is precisely because we are aiming at an equation in terms of the spatial gradients of the hydrodynamic fields, so once we have these we’re done.
Clearly the next step is to plug our expanded time derivatives into the right-hand side of the Boltzmann equation and crank the handle.
This gives us the (somewhat long!) equation
$$\begin{align} \sum_{n=0}^\infty \varepsilon^n \sum_{m=0}^n J(f^{(m)}, f^{(n-m)}) = &\sum_{n=1}^\infty \varepsilon^{n} \Big(\sum_{m=1}^{n} \big[\partial_\rho f^{(m-1)}\,\partial_t^{(n-m)}\rho + \partial_{\mathbf u} f^{(m-1)} \cdot \partial_t^{(n-m)} \mathbf u \\ &\qquad\qquad\qquad\vphantom{a} + \partial_T f^{(m-1)} \,\partial_t^{(n-m)} T \big] + \mathbf q \cdot \nabla f^{(n)}\Big) \end{align}$$Now we can group all the terms of equal order in $\varepsilon$ and form a series of equations
$$\begin{align} J(f^{(0)}, f^{(0)}) &= 0 & n = 0 \\ \sum_{i=0}^n J(f^{(i)}, f^{(n-i)}) &= \sum_{m=1}^{n} \big[\partial_\rho f^{(m-1)}\,\partial_t^{(n-m)}\rho + \partial_{\mathbf u} f^{(m-1)} \cdot \partial_t^{(n-m)} \mathbf u & \\ &\qquad\qquad\qquad\vphantom{a} + \partial_T f^{(m-1)} \,\partial_t^{(n-m)} T \big] + \mathbf q \cdot \nabla f^{(n)} & n > 0 \end{align}$$The crucial thing here is that terms of order $n$ depend only on the values of $f^{(i)}$ for $i\leq n$. This means we can sequentially solve terms upwards from order zero. Apparently this series converges quickly so it’s usual to only consider the first and second order terms. We already know that to first order we will get $f^{(0)} = f^\mathrm{eq}$ which is a Gaussian :
$$f^{(0)} = \frac{\rho}{\sqrt{(2\pi \theta)^d}}\exp\bigg(- \frac{|\mathbf v|^2}{2\theta}\bigg)$$with $\theta = k_B T / m$.
Second Order CE Expansion
Let’s look at the second order equation, $n=1$. This reads
$$\begin{align} 2 J(f^{(0)}, f^{(1)}) &= \partial_\rho f^{(0)} \partial_t^{(0)} \rho + \partial_{\mathbf u} f^{(0)} \cdot \partial_t^{(0)} \mathbf u + \partial_T f^{(0)} \partial_t^{(0)} T + \mathbf q \cdot \nabla f^{(0)} \tag{3}\label{eq:ce_first_order} \end{align}$$First note that the left-hand side is $f^{(0)} C(\phi)$ from the result for the linearized collision operator above.
Dealing with the right-hand side is not complicated but takes a fair amount of shuffling. To begin, let’s push the spatial derivative down onto the hydrodynamic variables in $f^{(0)}$
$$\nabla f^{(0)} = \partial_\rho f^{(0)} \nabla \rho + (\partial_{\mathbf u} f^{(0)} \cdot \nabla) \mathbf u + \partial_T f^{(0)} \nabla T$$Thus
$$\begin{align} f^{(0)} C(\phi) =\,\, &\partial_\rho f^{(0)} \partial_t^{(0)} \rho + \partial_{\mathbf u} f^{(0)} \cdot \partial_t^{(0)} \mathbf u + \partial_T f^{(0)} \partial_t^{(0)} T \\ \vphantom{x} + \mathbf q \cdot \big[ &\partial_\rho f^{(0)} \nabla \rho + (\partial_{\mathbf u} f^{(0)} \cdot \nabla) \mathbf u + \partial_T f^{(0)} \nabla T \big] \end{align}$$Collecting terms with like pre-factors gives
$$\begin{align} f^{(0)} C(\phi) &= \partial_\rho f^{(0)} (\partial_t^{(0)} \rho + \mathbf q\cdot \nabla \rho) \\ &\vphantom{.} + \partial_{\mathbf u} f^{(0)} \cdot (\partial_t^{(0)} \mathbf u + \mathbf q \cdot (\partial_{\mathbf u} f^{(0)} \cdot \nabla) \mathbf u) \\ &\vphantom{.} + \partial_T f^{(0)} (\partial_t^{(0)} T + \mathbf q \cdot \nabla T) \end{align}$$Interestingly this can be written as
$$f^{(0)}C(\phi) = \partial_\rho f^{(0)} D^{(0)} \rho + \partial_{\mathbf u} f^{(0)} \cdot D^{(0)} \mathbf u + \partial_T f^{(0)} D^{(0)} T$$where $D^{(0)} = \partial^{(0)}_t + \mathbf q \cdot \nabla$ is the zeroth order part of the total derivative (with zero acceleration). This is, in fact, related to a more direct but less explicit way to find the non-equilibrium terms.
Given the known form of $f^{(0)}$ we get the following results for its gradients
$$\begin{align} \frac{1}{f^{(0)}}\frac{\partial f^{(0)}}{\partial \rho} &= \frac{1}{\rho} \\ \frac{1}{f^{(0)}}\frac{\partial f^{(0)}}{\partial \mathbf u} &= \frac{\mathbf v}{\theta} \\ \frac{1}{f^{(0)}}\frac{\partial f^{(0)}}{\partial T} &= \frac{1}{2 T}\Big(\frac{v^2}{\theta} - 3\Big) \\ \end{align}$$Dividing by $f^{(0)}$ and plugging these in gives
$$\begin{align} C(\phi) &= \frac{1}{\rho} \big(\partial_t^{(0)} \rho + \mathbf q \cdot \nabla \rho\big) \\ &\vphantom{.} + \frac{\mathbf v}{\theta} \cdot \big(\partial_t^{(0)} \mathbf u + (\mathbf q \cdot \nabla) \mathbf u\big) \\ &\vphantom{.} + \frac{1}{2T}\Big(\frac{v^2}{\theta} - 3\Big) \big(\partial_t^{(0)} T +\mathbf q \cdot \nabla T\big) \end{align}$$where we have collected together terms with like pre-factors.
Now for the final conceptual step we can substitute in the zeroth order time derivative components. Above we noted that these are given by the Euler equations. Plugging these in gives us
$$\begin{align} C(\phi) &= \frac{1}{\rho} \big(-\nabla \cdot(\rho \mathbf u) + \mathbf q \cdot \nabla \rho\big) \\ &\vphantom{.} + \frac{\mathbf v}{\theta} \cdot \Big(-\frac{\nabla P}{\rho} - (\mathbf u\cdot \nabla)\mathbf u + (\mathbf q \cdot \nabla)\mathbf u\Big) \\ &\vphantom{.} + \frac{1}{2T}\Big(\frac{v^2}{\theta} - 3\Big) \Big(-\frac{2}{3} T \nabla\cdot\mathbf u - \mathbf u \cdot \nabla T + \mathbf q \cdot \nabla T\Big) \end{align}$$At this point, despite the sheer number of terms we are faced with, we have reached a physically interesting point. Note that we now have an equation for $\phi$ where $\phi$ appears only on the left-hand side, and the right-hand side consists entirely of spatial gradients of the hydrodynamic variables. We have successfully converted our first order equation \eqref{eq:ce_first_order}, involving time derivatives, into one using purely information known at the current time. Fantastic.
Before proceeding, let’s crank the simplification handle again for a while; you can view the steps by expanding the following box if you so desire.
Ultimately the simplifcation process yields
$$\boxed{ C(\phi) = \frac{m}{k_B T}\Big(\mathbf v \mathbf v - \frac{v^2}{3}\mathbf{I}\Big):\nabla\mathbf u + \Big(\frac{mv^2}{2k_B T} - \frac{5}{2}\Big)\mathbf{v} \cdot \nabla \mathrm{ln} T }$$This equation (recalling the integral nature of $C$) is an inhomogeneous Fredholm equation of the second kind . According to Harris the solution of this equation is a “truly gruesome task”. For that reason I won’t reproduce it here but it is detailed nicely in the textbook (section 6-5). It can be shown that the solution looks like
$$\phi = \frac{1}{\rho} \Big(\sqrt{\frac{2k_B T}{m}}\,\mathbf a \cdot \nabla \mathrm{ln}\,T + \mathbf{B}:\nabla \mathbf u\Big)$$where $\mathbf a(\mathbf q)$ is a vector and $\mathbf B(\mathbf q)$ is a traceless symmetric tensor, to be determined.
It is, apparently, possible to determine $\mathbf a$ and $\mathbf B$ directly but that is not overly useful. We do not care about $\phi$ directly, but rather its ability to let us calculate the moments of the first order expansion of $f$. For example now we know $\phi$ we can calculate
$$\boldsymbol\sigma^{(1)} = \langle \mathbf v \mathbf v \rangle^{(1)} = \langle \mathbf v \mathbf v \phi \rangle^{(0)} = \int \mathbf v \mathbf v \phi f^{(0)} \,d^3 q$$which results in
$$\boldsymbol\sigma^{(1)} = \mu \big(\nabla \mathbf u + (\nabla \mathbf u)^\mathrm{T} + \frac{2}{3}(\nabla\cdot\mathbf u)\mathbf I\big)$$similarly we can obtain
$$\mathbf j^{(1)} = - k \nabla T$$The values of the factors $\mu$ and $k$ depend on $\mathbf B$ and $\mathbf a$ respectively. They can be calculated, but their exact numerical values depends on the underlying microscopic model (e.g. hard spheres vs power-law repulsive potential vs Lennard-Jones potential ) via the collision cross-section appearing in the collision operator.
Combining the first (Euler equation) and second order terms and plugging into the equations of motion from part 1 yields, finally, the Navier-Stokes equation
$$\boxed{ \rho\frac{\mathrm d \mathbf u}{\mathrm d t} = -\nabla P + \mu \nabla\cdot \big(\nabla \mathbf u + (\nabla \mathbf u)^\mathrm{T} + \frac{2}{3}(\nabla\cdot\mathbf u)\mathbf I\big) + \mathbf g }$$Note that in the case of an incomressible fluid we have $\nabla\cdot\mathbf u=0$ so we get the slightly more familiar form
$$\rho\frac{\mathrm d \mathbf u}{\mathrm d t} = -\nabla P + \mu \nabla^2 \mathbf u + \rho \mathbf g$$where, remember, \(\mathbf g\) is an acceleration and $\rho\mathbf g$ is a force density.
This follows simply from
$$\begin{align} \big(\nabla\cdot\big(\nabla\mathbf u + (\nabla\mathbf u)^T\big) \big)_\alpha &= \partial_{x_\beta}(\partial_{x_\beta} u_\alpha + \partial_{x_\alpha} u_\beta ) \\ &= \partial_{x_\beta} \partial_{x_\beta} u_\alpha + \partial_{x_\alpha} \partial_{x_\beta} u_\beta \\ &= (\nabla^2 \mathbf u)_\alpha + \cancelto{0}{(\nabla(\nabla\cdot \mathbf u))_\alpha} \end{align}$$References
- An Introduction to the Theory of the Boltzmann Equation - Stewart Harris
- Existence of the Chapman-Enskog solution and its relation with first-order dissipative fluid theories - A. L. García-Perciante, et al (link )
- The Lattice Boltzmann Method, Principles and Practice - Timm Krueger, et al (link )
- The Mathematical Theory of Non-uniform Gases - Chapman & Cowling (link )
- The Methods of Chapman-Enskog and Grad and Applications - Gilberto Medeiros Kremer (PDF link )
- Time Derivative Expansion in Knudsen Number for the Chapman-Enskog expansion - @Syrocco (link )
I found all of these sources useful, but probably the clearest and most complete description of the C-E expansion was in Harris. It is certainly the one source I found myself referring to more than any of the others when trying to comprehend the material. That being said there are still some large jumps, and results that are stated but either not proved or at least the proofs are not referenced. I’ve tried to be very explicit here by showing all the steps (as I’ve understood them) and by justifying any ancillary results that are used.