Categories: Fluid mechanics, Physics, Thermodynamics.

Boltzmann equation

Consider a collection of particles, each with its own position r\vb{r} and velocity v\vb{v}. We can thus define a probability density function f(r,v,t)f(\vb{r}, \vb{v}, t) describing the expected particle count at (r,v)(\vb{r}, \vb{v}) at time tt. Let the total number of particles NN be conserved, then clearly:

N=∬−∞∞f(r,v,t)drdv\begin{aligned} N = \iint_{-\infty}^\infty f(\vb{r}, \vb{v}, t) \dd{\vb{r}} \dd{\vb{v}} \end{aligned}

At equilibrium, all processes affecting the particles no longer have a net effect, so ff is fixed:

dfdt=0\begin{aligned} \dv{f}{t} = 0 \end{aligned}

If each particle’s momentum only changes due to collisions, then a non-equilibrium state can be described as follows, very generally:

dfdt=( ⁣∂f∂t ⁣) ⁣col\begin{aligned} \dv{f}{t} = \bigg(\! \pdv{f}{t} \!\bigg)_\mathrm{\!col} \end{aligned}

Where the right-hand side simply means “all changes in ff due to collisions”. Applying the chain rule to the left-hand side then yields:

( ⁣∂f∂t ⁣) ⁣col=∂f∂t+(∂f∂xdxdt ⁣+ ⁣∂f∂ydydt ⁣+ ⁣∂f∂zdzdt)+(∂f∂vxdvxdt ⁣+ ⁣∂f∂vydvydt ⁣+ ⁣∂f∂vzdvzdt)=∂f∂t+(vx∂f∂x ⁣+ ⁣vy∂f∂y ⁣+ ⁣vz∂f∂z)+(ax∂f∂vx ⁣+ ⁣ay∂f∂vy ⁣+ ⁣az∂f∂vz)=∂f∂t+v⋅∇f+a⋅∂f∂v\begin{aligned} \bigg(\! \pdv{f}{t} \!\bigg)_\mathrm{\!col} &= \pdv{f}{t} + \bigg( \pdv{f}{x} \dv{x}{t} \!+\! \pdv{f}{y} \dv{y}{t} \!+\! \pdv{f}{z} \dv{z}{t} \bigg) + \bigg( \pdv{f}{v_x} \dv{v_x}{t} \!+\! \pdv{f}{v_y} \dv{v_y}{t} \!+\! \pdv{f}{v_z} \dv{v_z}{t} \bigg) \\ &= \pdv{f}{t} + \bigg( v_x \pdv{f}{x} \!+\! v_y \pdv{f}{y} \!+\! v_z \pdv{f}{z} \bigg) + \bigg( a_x \pdv{f}{v_x} \!+\! a_y \pdv{f}{v_y} \!+\! a_z \pdv{f}{v_z} \bigg) \\ &= \pdv{f}{t} + \vb{v} \cdot \nabla f + \vb{a} \cdot \pdv{f}{\vb{v}} \end{aligned}

Where we have introduced the shorthand ∂f/∂v\ipdv{f}{\vb{v}}. Inserting Newton’s second law F=ma\vb{F} = m \vb{a} leads us to the Boltzmann equation or Boltzmann transport equation (BTE):

∂f∂t+v⋅∇f+Fm⋅∂f∂v=( ⁣∂f∂t ⁣) ⁣col\begin{aligned} \boxed{ \pdv{f}{t} + \vb{v} \cdot \nabla f + \frac{\vb{F}}{m} \cdot \pdv{f}{\vb{v}} = \bigg(\! \pdv{f}{t} \!\bigg)_\mathrm{\!col} } \end{aligned}

But what about the collision term? Expressions for it exist, which are almost exact in many cases, but unfortunately also quite difficult to work with. In addition, ff is a 7-dimensional function, so the BTE is already hard to solve without collisions! We only present the simplest case, known as the Bhatnagar-Gross-Krook approximation: if the equilibrium state f0(r,v)f_0(\vb{r}, \vb{v}) is known, then each collision brings the system closer to f0f_0:

∂f∂t+v⋅∇f+Fm⋅∂f∂v=f0−fτ\begin{aligned} \pdv{f}{t} + \vb{v} \cdot \nabla f + \frac{\vb{F}}{m} \cdot \pdv{f}{\vb{v}} = \frac{f_0 - f}{\tau} \end{aligned}

Where τ\tau is the average collision period. The right-hand side is called the Krook term.

Moment equations

From the definition of ff, we see that integrating over all v\vb{v} yields the particle density nn:

n(r,t)=∫−∞∞f(r,v,t)dv\begin{aligned} n(\vb{r}, t) = \int_{-\infty}^\infty f(\vb{r}, \vb{v}, t) \dd{\vb{v}} \end{aligned}

Consequently, a purely velocity-dependent quantity Q(v)Q(\vb{v}) can be averaged like so:

⟨Q⟩=1n∫−∞∞Q(r,v,t) f(r,v,t)dv\begin{aligned} \Expval{Q} = \frac{1}{n} \int_{-\infty}^\infty Q(\vb{r}, \vb{v}, t) \: f(\vb{r}, \vb{v}, t) \dd{\vb{v}} \end{aligned}

With that in mind, we multiply the collisionless BTE equation by Q(v)Q(\vb{v}) and integrate, assuming that F\vb{F} does not depend on v\vb{v}:

0=∫−∞∞Q(∂f∂t+v⋅∇f+Fm⋅∂f∂v)dv=∫Q∂f∂tdv+∫(v⋅∇f) Qdv+Fm⋅∫Q∂f∂vdv=∂∂t∫Qfdv+∫(∇⋅(vf)−f(∇⋅v))Qdv+Fm⋅∫(∂∂v(Qf)−f∂Q∂v)dv\begin{aligned} 0 &= \int_{-\infty}^\infty Q \bigg( \pdv{f}{t} + \vb{v} \cdot \nabla f + \frac{\vb{F}}{m} \cdot \pdv{f}{\vb{v}} \bigg) \dd{\vb{v}} \\ &= \int Q \pdv{f}{t} \dd{\vb{v}} + \int (\vb{v} \cdot \nabla f) \: Q \dd{\vb{v}} + \frac{\vb{F}}{m} \cdot \int Q \pdv{f}{\vb{v}} \dd{\vb{v}} \\ &= \pdv{}{t}\int Q f \dd{\vb{v}} + \int \Big( \nabla \cdot (\vb{v} f) - f (\nabla \cdot \vb{v}) \Big) Q \dd{\vb{v}} + \frac{\vb{F}}{m} \cdot \int \bigg( \pdv{}{\vb{v}} (Q f) - f \pdv{Q}{\vb{v}} \bigg) \dd{\vb{v}} \end{aligned}

The first integral is simply n⟨Q⟩n \Expval{Q}. In the second integral, note that v\vb{v} is a coordinate and hence not dependent on r\vb{r}, so ∇⋅v=0\nabla \cdot \vb{v} = 0. Since ff is a probability density, f→0f \to 0 for v→±∞\vb{v} \to \pm\infty, so the first term in the third integral vanishes after it is integrated:

0=∂∂t(n⟨Q⟩)+∫∇⋅(vf) Qdv+Fm⋅([Qf]−∞∞−∫f∂Q∂vdv)=∂∂t(n⟨Q⟩)+∇⋅∫Qvfdv−Fm⋅∫f∂Q∂vdv\begin{aligned} 0 &= \pdv{}{t}\big(n \Expval{Q}\big) + \int \nabla \cdot (\vb{v} f) \: Q \dd{\vb{v}} + \frac{\vb{F}}{m} \cdot \bigg( \Big[ Q f \Big]_{-\infty}^\infty - \int f \pdv{Q}{\vb{v}} \dd{\vb{v}} \bigg) \\ &= \pdv{}{t}\big(n \Expval{Q}\big) + \nabla \cdot \int Q \vb{v} f \dd{\vb{v}} - \frac{\vb{F}}{m} \cdot \int f \pdv{Q}{\vb{v}} \dd{\vb{v}} \end{aligned}

We thus arrive at the prototype of the BTE’s so-called moment equations:

0=∂∂t(n⟨Q⟩)+∇⋅(n⟨Qv⟩)−Fm⋅(n⟨∂Q∂v⟩)\begin{aligned} \boxed{ 0 = \pdv{}{t}\big(n \Expval{Q}\big) + \nabla \cdot \big(n \Expval{Q \vb{v}}\big) - \frac{\vb{F}}{m} \cdot \bigg( n \Expval{\pdv{Q}{\vb{v}}} \bigg) } \end{aligned}

If we set Q=mQ = m, then the mass density ρ=n⟨Q⟩\rho = n \Expval{Q}, and we find that the zeroth moment of the BTE describes conservation of mass, where V≡⟨v⟩=∫vfdv\vb{V} \equiv \Expval{\vb{v}} = \int \vb{v} f \dd{\vb{v}} is the fluid velocity:

0=∂ρ∂t+∇⋅(ρV)\begin{aligned} \boxed{ 0 = \pdv{\rho}{t} + \nabla \cdot \big(\rho \vb{V}\big) } \end{aligned}

We insert Q=mQ = m into our prototype, and since mm is constant, the rest is trivial:

0=∂∂t(n⟨m⟩)+∇⋅(n⟨mv⟩)−Fm⋅(n⟨∂m∂v⟩)=∂ρ∂t+∇⋅(ρ⟨v⟩)−0\begin{aligned} 0 &= \pdv{}{t}\big(n \Expval{m}\big) + \nabla \cdot \big(n \Expval{m \vb{v}}\big) - \frac{\vb{F}}{m} \cdot \bigg( n \Expval{\pdv{m}{\vb{v}}} \bigg) \\ &= \pdv{\rho}{t} + \nabla \cdot \big(\rho \Expval{\vb{v}}\big) - 0 \end{aligned}

If we instead choose the momentum Q=mvQ = m \vb{v}, we find that the first moment of the BTE describes conservation of momentum, where P^\hat{P} is the Cauchy stress tensor:

0=∂∂t(ρV)+ρV(∇⋅V)+∇⋅P^−nF\begin{aligned} \boxed{ 0 = \pdv{}{t}\big(\rho \vb{V}\big) + \rho \vb{V} (\nabla \cdot \vb{V}) + \nabla \cdot \hat{P} - n \vb{F} } \end{aligned}

We insert Q=mvQ = m \vb{v} into our prototype and recognize ρ\rho wherever possible:

0=∂∂t(n⟨mv⟩)+∇⋅(n⟨mvv⟩)−Fm⋅(n⟨∂(mv)∂v⟩)=∂∂t(ρ⟨v⟩)+∇⋅(ρ⟨vv⟩)−F⋅(n⟨∂v∂v⟩)\begin{aligned} 0 &= \pdv{}{t}\big(n \Expval{m \vb{v}}\big) + \nabla \cdot \big(n \Expval{m \vb{v} \vb{v}}\big) - \frac{\vb{F}}{m} \cdot \bigg( n \Expval{\pdv{(m \vb{v})}{\vb{v}}} \bigg) \\ &= \pdv{}{t}\big(\rho \Expval{\vb{v}}\big) + \nabla \cdot \big(\rho \Expval{\vb{v} \vb{v}}\big) - \vb{F} \cdot \bigg( n \Expval{\pdv{\vb{v}}{\vb{v}}} \bigg) \end{aligned}

With vv\vb{v} \vb{v} being a dyadic product. To give it a physical interpretation, we split v=V ⁣+ ⁣w\vb{v} = \vb{V} \!+\! \vb{w}, where V\vb{V} is the average velocity vector, and w\vb{w} is the local deviation from V\vb{V}:

⟨vv⟩=⟨(V ⁣+ ⁣w)(V ⁣+ ⁣w)⟩=⟨VV+2Vw+ww⟩=VV+2V⟨w⟩+⟨ww⟩\begin{aligned} \Expval{\vb{v} \vb{v}} &= \Expval{(\vb{V} \!+\! \vb{w}) (\vb{V} \!+\! \vb{w})} = \Expval{\vb{V} \vb{V} + 2 \vb{V} \vb{w} + \vb{w} \vb{w}} = \vb{V} \vb{V} + 2 \vb{V} \Expval{\vb{w}} + \Expval{\vb{w} \vb{w}} \end{aligned}

Since w\vb{w} represents a deviation from the mean, ⟨w⟩=0\Expval{\vb{w}} = 0. We define the pressure tensor:

P^≡ρ⟨ww⟩=ρ⟨(v ⁣− ⁣V)(v ⁣− ⁣V)⟩\begin{aligned} \hat{P} \equiv \rho \Expval{\vb{w} \vb{w}} = \rho \Expval{(\vb{v} \!-\! \vb{V}) (\vb{v} \!-\! \vb{V})} \end{aligned}

This leads to the desired result, where ∇⋅(ρVV)\nabla \cdot (\rho \vb{V}\vb{V}) is the fluid momentum, and ∇⋅P^\nabla \cdot \hat{P} is the viscous/pressure momentum:

0=∂∂t(ρV)+∇⋅(ρVV+P^)−nF\begin{aligned} 0 &= \pdv{}{t}\big(\rho \vb{V}\big) + \nabla \cdot \big(\rho \vb{V} \vb{V} + \hat{P}\big) - n \vb{F} \end{aligned}

Finally, if we choose the kinetic energy Q=m∣v∣2/2Q = m |\vb{v}|^2 / 2, we find that the second moment gives conservation of energy, where UU is the thermal energy density and J\vb{J} is the heat flux:

0=∂∂t(ρ2∣V∣2+U)+∇⋅(ρ2∣V∣2V+V⋅P^+UV+J)−F⋅(nV)\begin{aligned} \boxed{ 0 = \pdv{}{t}\bigg(\frac{\rho}{2} |\vb{V}|^2 + U \bigg) + \nabla \cdot \bigg(\frac{\rho}{2} |\vb{V}|^2 \vb{V} + \vb{V} \cdot \hat{P} + U \vb{V} + \vb{J} \bigg) - \vb{F} \cdot \big( n \vb{V} \big) } \end{aligned}

We insert Q=m∣v∣2/2Q = m |\vb{v}|^2 / 2 into our prototype and recognize ρ\rho wherever possible:

0=∂∂t(n⟨m∣v∣22⟩)+∇⋅(n⟨m∣v∣22v⟩)−Fm⋅(n⟨∂∂vm∣v∣22⟩)=∂∂t(ρ2⟨∣v∣2⟩)+∇⋅(ρ2⟨∣v∣2v⟩)−F2⋅(n⟨∂∣v∣2∂v⟩)\begin{aligned} 0 &= \pdv{}{t}\bigg(n \Expval{\frac{m |\vb{v}|^2}{2}}\bigg) + \nabla \cdot \bigg(n \Expval{\frac{m |\vb{v}|^2}{2} \vb{v}}\bigg) - \frac{\vb{F}}{m} \cdot \bigg( n \Expval{\pdv{}{\vb{v}} \frac{m |\vb{v}|^2}{2}} \bigg) \\ &= \pdv{}{t}\bigg(\frac{\rho}{2} \Expval{|\vb{v}|^2}\bigg) + \nabla \cdot \bigg(\frac{\rho}{2} \Expval{|\vb{v}|^2 \vb{v}}\bigg) - \frac{\vb{F}}{2} \cdot \bigg( n \Expval{\pdv{|\vb{v}|^2}{\vb{v}}} \bigg) \end{aligned}

We handle these terms one by one. Substituting v=V+w\vb{v} = \vb{V} + \vb{w} in the first gives:

⟨∣v∣2⟩=⟨(V ⁣+ ⁣w)⋅(V ⁣+ ⁣w)⟩=⟨∣V∣2+2V⋅w+∣w∣2⟩=∣V∣2+2V⋅⟨w⟩+⟨∣w∣2⟩=∣V∣2+⟨∣w∣2⟩\begin{aligned} \Expval{|\vb{v}|^2} &= \Expval{(\vb{V} \!+\! \vb{w}) \cdot (\vb{V} \!+\! \vb{w})} = \Expval{|\vb{V}|^2 + 2 \vb{V} \cdot \vb{w} + |\vb{w}|^2} \\ &= |\vb{V}|^2 + 2 \vb{V} \cdot \Expval{\vb{w}} + \Expval{|\vb{w}|^2} = |\vb{V}|^2 + \Expval{|\vb{w}|^2} \end{aligned}

And likewise for the second term, where we recognize the stress tensor ⟨ww⟩\Expval{\vb{w} \vb{w}}:

⟨∣v∣2v⟩=⟨(V ⁣+ ⁣w)⋅(V ⁣+ ⁣w)(V ⁣+ ⁣w)⟩=⟨(∣V∣2+2V⋅w+∣w∣2)(V ⁣+ ⁣w)⟩=⟨∣V∣2V+∣V∣2w+2(V⋅w)V+2(V⋅w)w+∣w∣2V+∣w∣2w⟩=∣V∣2V+∣V∣2⟨w⟩+2(V⋅⟨w⟩)V+2⟨(V⋅w)w⟩+⟨∣w∣2⟩V+⟨∣w∣2w⟩=∣V∣2V+0+0+2V⋅⟨ww⟩+⟨∣w∣2⟩V+⟨∣w∣2w⟩\begin{aligned} \Expval{|\vb{v}|^2 \vb{v}} &= \Expval{(\vb{V} \!+\! \vb{w}) \cdot (\vb{V} \!+\! \vb{w}) (\vb{V} \!+\! \vb{w})} = \Expval{(|\vb{V}|^2 + 2 \vb{V} \cdot \vb{w} + |\vb{w}|^2) (\vb{V} \!+\! \vb{w})} \\ &= \Expval{|\vb{V}|^2 \vb{V} + |\vb{V}|^2 \vb{w} + 2 (\vb{V} \cdot \vb{w}) \vb{V} + 2 (\vb{V} \cdot \vb{w}) \vb{w} + |\vb{w}|^2 \vb{V} + |\vb{w}|^2 \vb{w}} \\ &= |\vb{V}|^2 \vb{V} + |\vb{V}|^2 \Expval{\vb{w}} + 2 (\vb{V} \cdot \Expval{\vb{w}}) \vb{V} + 2 \Expval{(\vb{V} \cdot \vb{w}) \vb{w}} + \Expval{|\vb{w}|^2} \vb{V} + \Expval{|\vb{w}|^2 \vb{w}} \\ &= |\vb{V}|^2 \vb{V} + 0 + 0 + 2 \vb{V} \cdot \Expval{\vb{w} \vb{w}} + \Expval{|\vb{w}|^2} \vb{V} + \Expval{|\vb{w}|^2 \vb{w}} \end{aligned}

The third term is fairly obvious, but we calculate it rigorously just to be safe:

∂∣v∣2∂v=∂∂v(vx2+vy2+vz2)=e^x∂vx2∂vx+e^y∂vy2∂vy+e^z∂vz2∂vz=2v\begin{aligned} \pdv{|\vb{v}|^2}{\vb{v}} &= \pdv{}{\vb{v}} \big( v_x^2 + v_y^2 + v_z^2 \big) = \vu{e}_x \pdv{v_x^2}{v_x} + \vu{e}_y \pdv{v_y^2}{v_y} + \vu{e}_z \pdv{v_z^2}{v_z} = 2 \vb{v} \end{aligned}

To clarify the physical interpretation, we define UU, J\vb{J} and P^\hat{P} as follows:

U≡ρ2⟨∣w∣2⟩=ρ2⟨(v ⁣− ⁣V)⋅(v ⁣− ⁣V)⟩J≡ρ2⟨∣w∣2w⟩=ρ2⟨(v ⁣− ⁣V)⋅(v ⁣− ⁣V)(v ⁣− ⁣V)⟩P^≡ρ⟨ww⟩=ρ⟨(v ⁣− ⁣V)(v ⁣− ⁣V)⟩\begin{aligned} U &\equiv \frac{\rho}{2} \Expval{|\vb{w}|^2} = \frac{\rho}{2} \Expval{(\vb{v} \!-\! \vb{V}) \cdot (\vb{v} \!-\! \vb{V})} \\ \vb{J} &\equiv \frac{\rho}{2} \Expval{|\vb{w}|^2 \vb{w}} = \frac{\rho}{2} \Expval{(\vb{v} \!-\! \vb{V}) \cdot (\vb{v} \!-\! \vb{V})(\vb{v} \!-\! \vb{V})} \\ \hat{P} &\equiv \rho \Expval{\vb{w} \vb{w}} = \rho \Expval{(\vb{v} \!-\! \vb{V}) (\vb{v} \!-\! \vb{V})} \end{aligned}

Putting it all together, we arrive at the expected result, namely:

0=∂∂t(ρ2∣V∣2+U)+∇⋅(ρ2∣V∣2V+V⋅P^+UV+J)−F⋅(nV)\begin{aligned} 0 &= \pdv{}{t}\bigg(\frac{\rho}{2} |\vb{V}|^2 + U \bigg) + \nabla \cdot \bigg(\frac{\rho}{2} |\vb{V}|^2 \vb{V} + \vb{V} \cdot \hat{P} + U \vb{V} + \vb{J} \bigg) - \vb{F} \cdot \big( n \vb{V} \big) \end{aligned}

For the sake of clarity, we write out the pressure term, including the outer divergence:

∇⋅(V⋅P^)=(∇⋅P^⊤)⋅V=∇⋅[PxxPxyPxzPyxPyyPyzPzxPzyPzz]⋅[VxVyVz]=[∂Pxx∂x+∂Pxy∂y+∂Pxz∂z∂Pyx∂x+∂Pyy∂y+∂Pyz∂z∂Pzx∂x+∂Pzy∂y+∂Pzz∂z]⊤⋅[VxVyVz]=∑i=13∑j=13∂Pij∂xjVi\begin{aligned} \nabla \cdot (\vb{V} \cdot \hat{P}) &= (\nabla \cdot \hat{P}{}^\top) \cdot \vb{V} = \nabla \cdot \begin{bmatrix} P_{xx} & P_{xy} & P_{xz} \\ P_{yx} & P_{yy} & P_{yz} \\ P_{zx} & P_{zy} & P_{zz} \end{bmatrix} \cdot \begin{bmatrix} V_x \\ V_y \\ V_z \end{bmatrix} \\ &= \begin{bmatrix} \displaystyle \pdv{P_{xx}}{x} + \pdv{P_{xy}}{y} + \pdv{P_{xz}}{z} \\ \displaystyle \pdv{P_{yx}}{x} + \pdv{P_{yy}}{y} + \pdv{P_{yz}}{z} \\ \displaystyle \pdv{P_{zx}}{x} + \pdv{P_{zy}}{y} + \pdv{P_{zz}}{z} \end{bmatrix}^{\top} \cdot \begin{bmatrix} V_x \\ V_y \\ V_z \end{bmatrix} = \sum_{i=1}^{3} \sum_{j=1}^{3} \pdv{P_{ij}}{x_j} V_i \end{aligned}

References

  1. M. Salewski, A.H. Nielsen, Plasma physics: lecture notes, 2021, unpublished.