Fluid Approximation
The fluid approximation is valid when the system size $L$ is much larger than the mean free path $l_\text{mfp}$.
For ISM (interstellar medium), the typical mean free path is
l mfp ∼ 1 0 15 cm ( n 1 cm − 3 ) − 1 ( d 1 A ˚ ) − 2 l_\text{mfp}\sim 10^{15}\text{cm}\left(\frac{n}{1\text{ cm}^{-3}}\right)^{-1}\left(\frac{d}{1\ \AA}\right)^{-2} l mfp ∼ 1 0 15 cm ( 1 cm − 3 n ) − 1 ( 1 A ˚ d ) − 2 Inside each fluid particle, we can define thermal dynamical quantities, $T,\rho,P,S$, as well as a velocity field $\vec v(\vec x,t)$.
For any quantity of the fluid element at $(\vec r_0,t_0)$, say $F(\vec r_0,t_0)$, we have
Δ F = F ( r ⃗ 0 + v ⃗ ( r ⃗ 0 ) Δ t , t 0 + Δ t ) − F ( r ⃗ 0 , t 0 ) = ( ∂ F ∂ t + v x ∂ F ∂ x + v y ∂ F ∂ y + v z ∂ F ∂ z ) Δ t + O ( Δ t 2 ) \begin{align*}
\Delta F&=F\left(\vec r_0+\vec v(\vec r_0)\Delta t,t_0+\Delta t\right)-F(\vec r_0,t_0)\\
&=\left(\frac{\partial F}{\partial t}+v_x\frac{\partial F}{\partial x}+v_y\frac{\partial F}{\partial y}+v_z\frac{\partial F}{\partial z}\right)\Delta t+\mathcal O(\Delta t^2)
\end{align*} Δ F = F ( r 0 + v ( r 0 ) Δ t , t 0 + Δ t ) − F ( r 0 , t 0 ) = ( ∂ t ∂ F + v x ∂ x ∂ F + v y ∂ y ∂ F + v z ∂ z ∂ F ) Δ t + O ( Δ t 2 ) In fact, we can write the derivative of $F$
lim Δ t → 0 Δ F Δ t \lim_{\Delta t\to 0}\frac{\Delta F}{\Delta t} Δ t → 0 lim Δ t Δ F in two different ways
Lagrange's derivative
D F D t \frac{\text DF}{\text Dt} D t D F which follows the fluid element everywhere.
Eulerian derivative
( ∂ ∂ t + v ⃗ ⋅ ∇ ) F \left(\frac{\partial}{\partial t}+\vec v\cdot \nabla\right) F ( ∂ t ∂ + v ⋅ ∇ ) F which is more important for numerical calculations, as it considers only fixed fluid cells.
Equation of Continuity - Mass Conservation
In Eulerian picture
For a small, fixed volume $\delta V_0$, in a time span of $\Delta t$, the change in mass inside is
δ m = Δ t ⋅ d d t ( ∭ δ V 0 ρ d V 0 ) \delta m=\Delta t\cdot\frac{\text d}{\text dt}\left(\iiint_{\delta V_0}\rho\text dV_0\right) δ m = Δ t ⋅ d t d ( ∭ δ V 0 ρ d V 0 ) Consider the mass flow at the surface
j ⃗ = Δ t ∯ δ S ( v ⃗ ⋅ n ⃗ ) ρ d S = Δ t ∭ δ V 0 ∇ ⋅ ( ρ v ⃗ ) d V 0 \vec j=\Delta t\oiint_{\delta S}\left(\vec v\cdot\vec n\right)\rho\text dS=\Delta t\iiint_{\delta V_0}\nabla\cdot\left(\rho\vec v\right)\text dV_0 j = Δ t ∬ δ S ( v ⋅ n ) ρ d S = Δ t ∭ δ V 0 ∇ ⋅ ( ρ v ) d V 0 where we have applied Gauss' divergence theorem. Here $\vec n$ is the normal vector of the surface.
If there are no source/sink of mass within the fluid element, the mass flow should exactly account for the change in mass
d d t ( ∭ δ V 0 ρ d V 0 ) + ∭ δ V 0 ∇ ⋅ ( ρ v ⃗ ) d V 0 = 0 \frac{\text d}{\text dt}\left(\iiint_{\delta V_0}\rho\text dV_0\right)+\iiint_{\delta V_0}\nabla\cdot\left(\rho\vec v\right)\text dV_0=0 d t d ( ∭ δ V 0 ρ d V 0 ) + ∭ δ V 0 ∇ ⋅ ( ρ v ) d V 0 = 0 Since we have fixed $\Delta V$, we can write this formula so that for any $\Delta V_0$
∭ δ V 0 ( ∂ ρ ∂ t + ∇ ⋅ ( ρ v ⃗ ) ) d V 0 = 0 \iiint_{\delta V_0}\left(\frac{\partial \rho}{\partial t}+\nabla\cdot\left(\rho\vec v\right)\right)\text dV_0=0 ∭ δ V 0 ( ∂ t ∂ ρ + ∇ ⋅ ( ρ v ) ) d V 0 = 0 The arbitrariness of the choice of $\Delta V$ directly leads to the famous equation of continuity
∂ ρ ∂ t + ∇ ⋅ ( ρ v ⃗ ) = 0 \frac{\partial \rho}{\partial t}+\nabla\cdot\left(\rho\vec v\right)=0 ∂ t ∂ ρ + ∇ ⋅ ( ρ v ) = 0 In Lagrange's picture
Let us consider a volume element $\delta V$ that moves with fluid particles.
If no mass creation / annihilation
0 = D ( δ m ) D t = D D t ( ρ δ V ) = D ρ D t δ V + ρ D D t ( δ V ) 0=\frac{\text D(\delta m)}{\text Dt}=\frac{\text D}{\text Dt}\left(\rho\delta V\right)=\frac{\text D\rho}{\text Dt}\delta V+\rho\frac{\text D}{\text Dt}\left(\delta V\right) 0 = D t D ( δ m ) = D t D ( ρ δ V ) = D t D ρ δ V + ρ D t D ( δ V ) ⇒ D ρ D t = − ρ δ V D D t ( δ V ) = − ρ ∑ i 1 δ x i D ( δ x i ) D t = − ρ ( ∇ ⋅ v ⃗ ) \Rightarrow
\frac{\text D\rho}{\text Dt}=-\frac{\rho}{\delta V}\frac{\text D}{\text Dt}\left(\delta V\right)=-\rho\sum_{i}\frac{1}{\delta x_i}\frac{\text D\left(\delta x_i\right)}{\text Dt}=-\rho\left(\nabla\cdot \vec v\right) ⇒ D t D ρ = − δ V ρ D t D ( δ V ) = − ρ i ∑ δ x i 1 D t D ( δ x i ) = − ρ ( ∇ ⋅ v ) Recall that
D D t = ∂ ∂ t + v ⃗ ⋅ ∇ \frac{\text D}{\text Dt}=\frac{\partial}{\partial t}+\vec v\cdot \nabla D t D = ∂ t ∂ + v ⋅ ∇ Then
( ∂ ∂ t + v ⃗ ⋅ ∇ + ∇ ⋅ v ⃗ ) ρ = 0 ⇒ ∂ ρ ∂ t + ∇ ⋅ ( ρ v ⃗ ) = 0 \left(\frac{\partial}{\partial t}+\vec v\cdot \nabla+\nabla\cdot\vec v\right)\rho=0\Rightarrow\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\vec v\right)=0 ( ∂ t ∂ + v ⋅ ∇ + ∇ ⋅ v ) ρ = 0 ⇒ ∂ t ∂ ρ + ∇ ⋅ ( ρ v ) = 0 Two pictures give exactly the same equation!
Now we further explore the relation between two pictures. In the Lagrange's picture,
δ V = d x ⋅ d y ⋅ d z \delta V=\text dx\cdot\text dy\cdot\text dz δ V = d x ⋅ d y ⋅ d z while in the Eulerian picture
δ V 0 = d ξ x ⋅ d ξ y ⋅ d ξ z \delta V_0=\text d\xi_x\cdot\text d\xi_y\cdot\text d\xi_z δ V 0 = d ξ x ⋅ d ξ y ⋅ d ξ z $\delta V$ and $\delta V_0$ are related as follows,
δ V = ∣ ∂ ( x , y , z ) ∂ ( ξ x , ξ y , ξ z ) ∣ δ V 0 ≡ J δ V 0 \delta V=\left|\frac{\partial\left(x,y,z\right)}{\partial\left(\xi_x,\xi_y,\xi_z\right)}\right|\delta V_0\equiv J\delta V_0 δ V = ∂ ( ξ x , ξ y , ξ z ) ∂ ( x , y , z ) δ V 0 ≡ J δ V 0 The ratio $J$ is called the expansion of the fluid. Euler's expansion formula claims that
D J D t = ( ∇ ⋅ v ⃗ ) J \frac{\text DJ}{\text Dt}=\left(\nabla \cdot \vec v\right)J D t D J = ( ∇ ⋅ v ) J Then if we reconsider the conservation of mass in Lagrange's picture,
D D t ∭ δ V ρ d V = 0 \frac{\text D}{\text Dt}\iiint_{\delta V}\rho\text{d}V=0 D t D ∭ δ V ρ d V = 0 ⟺ 0 = D D t ∭ δ V ρ J d V 0 = ∭ δ V [ D ρ D t + ρ ( ∇ ⋅ v ⃗ ) ] J d V 0 \iff 0=\frac{\text D}{\text Dt}\iiint_{\delta V}\rho J\text{d}V_0=\iiint_{\delta V}\left[\frac{\text D \rho}{\text Dt}+\rho\left(\nabla\cdot \vec v\right)\right]J\text dV_0 ⟺ 0 = D t D ∭ δ V ρ J d V 0 = ∭ δ V [ D t D ρ + ρ ( ∇ ⋅ v ) ] J d V 0 which also gives the continuity formula
D ρ D t + ρ ( ∇ ⋅ u ) = 0 \frac{\text D \rho}{\text Dt}+\rho\left(\nabla\cdot u\right)=0 D t D ρ + ρ ( ∇ ⋅ u ) = 0 Before moving to the next section, we can similarly derive several useful identities of the Lagrange's derivative.
For a function $F$,
D D t ∭ V F d V = D D t ∭ V F J d V 0 = ∭ V [ D F D t + ( ∇ ⋅ v ⃗ ) F ] J d V 0 = ∭ V [ D F D t + ( ∇ ⋅ v ⃗ ) F ] d V \frac{\text D }{\text Dt}\iiint_VF\text dV=\frac{\text D }{\text Dt}\iiint_VFJ\text dV_0=\iiint_V\left[\frac{\text D F}{\text Dt}+\left(\nabla \cdot \vec v\right)F\right]J\text dV_0=\iiint_V\left[\frac{\text D F}{\text Dt}+\left(\nabla \cdot \vec v\right)F\right]\text dV D t D ∭ V F d V = D t D ∭ V F J d V 0 = ∭ V [ D t D F + ( ∇ ⋅ v ) F ] J d V 0 = ∭ V [ D t D F + ( ∇ ⋅ v ) F ] d V ⟺ D D t ∭ V F d V = ∭ V [ ∂ F ∂ t + ∇ ⋅ ( F v ⃗ ) ] d V \iff \frac{\text D }{\text Dt}\iiint_VF\text dV=\iiint_V\left[\frac{\partial F}{\partial t}+\nabla\cdot\left(F \vec v\right)\right]\text dV ⟺ D t D ∭ V F d V = ∭ V [ ∂ t ∂ F + ∇ ⋅ ( F v ) ] d V This is known as the Reynolds transport theorem . Simply let $F=\alpha\rho$, we have
D D t ∭ V α ρ d V = ∭ V [ D D t ( α ρ ) + ( ∇ ⋅ v ⃗ ) α ρ ] d V = ∭ V { α [ D ρ D t + ( ∇ ⋅ v ⃗ ) ρ ] + ρ D α D t } d V = ∭ V ρ D α D t d V \begin{align*}
\frac{\text D }{\text Dt}\iiint_V\alpha\rho\text dV&=\iiint_V\left[\frac{\text D }{\text Dt}\left(\alpha\rho\right)+\left(\nabla \cdot \vec v\right)\alpha\rho\right]\text dV\\
&=\iiint_V\left\{\alpha\left[\frac{\text D \rho}{\text Dt}+\left(\nabla \cdot \vec v\right)\rho\right]+\rho\frac{\text D\alpha}{\text Dt}\right\}\text dV\\
&=\iiint_V\rho\frac{\text D\alpha}{\text Dt}\text dV
\end{align*} D t D ∭ V α ρ d V = ∭ V [ D t D ( α ρ ) + ( ∇ ⋅ v ) α ρ ] d V = ∭ V { α [ D t D ρ + ( ∇ ⋅ v ) ρ ] + ρ D t D α } d V = ∭ V ρ D t D α d V This identity is extremely useful in the following chapters.
Euler Equation - Momentum Conservation
In the Lagrange's picture, the acceleration onto a mass element is given by
D D t ∭ V ( t ) ρ v ⃗ d V = ∭ V ( t ) ρ D v ⃗ D t d V \frac{\text D}{\text Dt}\iiint_{V(t)}\rho\vec v\text dV=\iiint_{V(t)}\rho\frac{\text D\vec v}{\text Dt}\text dV D t D ∭ V ( t ) ρ v d V = ∭ V ( t ) ρ D t D v d V With the Reynolds transport theorem, the EoS thus gives
∭ V ( t ) ρ D v ⃗ D t d V = ∭ V ( t ) f ⃗ d V + ∯ S ( t ) t ⃗ d S \iiint_{V(t)}\rho\frac{\text D\vec v}{\text Dt}\text dV=\iiint_{V(t)}\vec f\text dV+\oiint_{S(t)}\vec t\text dS ∭ V ( t ) ρ D t D v d V = ∭ V ( t ) f d V + ∬ S ( t ) t d S First term - body force (e.g gravity)
Second term - surface force (e.g. pressure)
Consider the $i$-th component of the RHS
∭ V ( t ) f i d V + ∯ S ( t ) t i d S = ∭ V ( t ) f i d V + ∯ S ( t ) T i j n j d S = ∭ V ( t ) ( f i + ∂ j T i j ) d V \iiint_{V(t)}f^i\text dV+\oiint_{S(t)}t^i\text dS=\iiint_{V(t)}f^i\text dV+\oiint_{S(t)}T^{ij}n_j\text dS=\iiint_{V(t)}\left(f^i+\partial_jT^{ij}\right)\text dV ∭ V ( t ) f i d V + ∬ S ( t ) t i d S = ∭ V ( t ) f i d V + ∬ S ( t ) T ij n j d S = ∭ V ( t ) ( f i + ∂ j T ij ) d V where $\mathcal T$ is the stress tensor , which satisfies
t i = T i j n j t^i=T^{ij}n_j t i = T ij n j An important surface force is pressure.
The force caused by pressure is
F ⃗ p = − ∬ S p n ⃗ d S = − ∭ V ∇ p d V \vec F_p=-\iint_S p\vec n\text dS=-\iiint_V\nabla p\text dV F p = − ∬ S p n d S = − ∭ V ∇ p d V Pressure and Stress Tensor
If the $i$-th component
− ∂ i p = ∂ j T i j -\partial_i p=\partial_j T^{ij} − ∂ i p = ∂ j T ij Thus $p$ is strongly correlated with the diagonal components in $\mathcal T$.
Since the choice of $V$ is arbitrary, the EoM is expressed as
ρ D v ⃗ D t = − ∇ p + f ⃗ \rho\frac{\text D\vec v}{\text Dt}=-\nabla p+\vec f ρ D t D v = − ∇ p + f ⟺ ∂ v ⃗ ∂ t + ( v ⃗ ⋅ ∇ ) v ⃗ = − 1 ρ ∇ p + 1 ρ f ⃗ \iff \frac{\partial \vec v}{\partial t}+\left(\vec v\cdot \nabla\right)\vec v=-\frac1\rho\nabla p+\frac1\rho\vec f ⟺ ∂ t ∂ v + ( v ⋅ ∇ ) v = − ρ 1 ∇ p + ρ 1 f which is the famous Euler equation for ideal fluid.
Using the equation of continuity and the Euler equation, we can calculate the time derivative of mass flux
∂ ∂ t ( ρ v i ) = v i ∂ ρ ∂ t + ρ ∂ v i ∂ t = − v i ∂ j ( ρ v j ) − ρ v j ∂ j v i − ∂ j δ i j p = − ∂ j ( ρ v i v j + δ i j p ) ≡ − ∂ j Π i j \begin{align*}
\frac{\partial}{\partial t}\left(\rho v^i\right)&=v^i\frac{\partial\rho}{\partial t}+\rho\frac{\partial v^i}{\partial t}\\
&=-v^i\partial_j\left(\rho v^j\right)-\rho v^j\partial_jv^i-\partial_j\delta^{ij}p\\
&=-\partial_j\left(\rho v^iv^j+\delta^{ij}p\right)\\
&\equiv -\partial_j\Pi^{ij}
\end{align*} ∂ t ∂ ( ρ v i ) = v i ∂ t ∂ ρ + ρ ∂ t ∂ v i = − v i ∂ j ( ρ v j ) − ρ v j ∂ j v i − ∂ j δ ij p = − ∂ j ( ρ v i v j + δ ij p ) ≡ − ∂ j Π ij So in the Eulerian picture, the momentum conservation is
∂ ∂ t ( ρ v i ) + ∂ j Π i j = 0 \frac{\partial}{\partial t}\left(\rho v^i\right)+\partial_j\Pi^{ij}=0 ∂ t ∂ ( ρ v i ) + ∂ j Π ij = 0 Integrate the above equation over a fixed volume $V$
d d t ∭ V ρ v i d V = − ∭ ∂ i p d V − ∭ ∂ j ( ρ v i v j ) d V = − ∯ S p n i d S − ∯ S ρ v i v j n j d S \begin{align*}
\frac{\text d}{\text dt}\iiint_V\rho v^i\text dV&=-\iiint\partial^i p\text dV-\iiint\partial_j\left(\rho v^iv^j\right)\text dV\\
&=-\oiint_Spn^i\text dS-\oiint_S\rho v^iv^jn_j\text dS
\end{align*} d t d ∭ V ρ v i d V = − ∭ ∂ i p d V − ∭ ∂ j ( ρ v i v j ) d V = − ∬ S p n i d S − ∬ S ρ v i v j n j d S First term - thermal pressure
Second term - ram pressure
Equation of State: $f(\rho,p,T)=0$
Two independent quantities are necessary to characterize thermal states.
Ideal gas
p = ρ k B T μ m p p=\frac{\rho k_BT}{\mu m_p} p = μ m p ρ k B T Adiabatic gas
p = K ρ γ , γ = c p c V p=K\rho^\gamma,\quad \gamma=\frac{c_p}{c_V} p = K ρ γ , γ = c V c p The adiabatic index $\gamma=5/3$ for mono-atomic gas.
Ideal + adiabatic is an example of polytropic EoS
If $p=f(\rho)$, the EoS is known as barotropic .
For inviscid fluid
D s D t = ∂ s ∂ t + ( v ⃗ ⋅ ∇ ) s = 0 \frac{\text Ds}{\text Dt}=\frac{\partial s}{\partial t}+\left(\vec v\cdot \nabla\right)s=0 D t D s = ∂ t ∂ s + ( v ⋅ ∇ ) s = 0 that is, entropy doesn't change along each fluid stream line.
Furthermore, if $s=const$ everywhere at $r=r_0$, and $s=const$ is preserved, the. motion is known as isentropic motion .
Bernoulli's theorem
In general, we would like to derive the energy conservation with Euler equation with certain assumptions
Barotropic EoS: $\rho=f(p)$, thus the specific enthalpy is simply
h ( p ) ≡ ∫ d p ρ ( p ) h(p)\equiv\int\frac{\text dp}{\rho(p)} h ( p ) ≡ ∫ ρ ( p ) d p Potential force: $\vec f=-\rho\nabla \phi$
Then we can rewrite Euler equation
∂ v ⃗ ∂ t + ( v ⃗ ⋅ ∇ ) v ⃗ = − 1 ρ ∇ p − ∇ ϕ \frac{\partial \vec v}{\partial t}+\left(\vec v\cdot \nabla\right)\vec v=-\frac1\rho\nabla p-\nabla\phi ∂ t ∂ v + ( v ⋅ ∇ ) v = − ρ 1 ∇ p − ∇ ϕ ⇒ ∂ v ⃗ ∂ t + ∇ ( 1 2 v 2 + h + ϕ ) = v ⃗ × ω ⃗ \Rightarrow \frac{\partial \vec v}{\partial t}+\nabla\left(\frac12v^2+h+\phi\right)=\vec v\times\vec\omega ⇒ ∂ t ∂ v + ∇ ( 2 1 v 2 + h + ϕ ) = v × ω where the vorticity is defined as
ω ⃗ = ∇ × v ⃗ \vec\omega=\nabla\times\vec v ω = ∇ × v and we have applied the identity
v ⃗ × ( ∇ × v ⃗ ) = ( ∇ ⋅ v ⃗ ) v ⃗ − ( v ⃗ ⋅ ∇ ) v ⃗ = 1 2 ∇ v 2 − ( v ⃗ ⋅ ∇ ) v ⃗ \vec v\times\left(\nabla\times\vec v\right)=\left(\nabla\cdot\vec v\right) \vec v-\left(\vec v\cdot\nabla\right) \vec v=\frac12\nabla v^2-\left(\vec v\cdot\nabla\right) \vec v v × ( ∇ × v ) = ( ∇ ⋅ v ) v − ( v ⋅ ∇ ) v = 2 1 ∇ v 2 − ( v ⋅ ∇ ) v For steady flow
∂ u ∂ t = 0 \frac{\partial u}{\partial t}=0 ∂ t ∂ u = 0 thus
v ⃗ ⋅ ∇ ( 1 2 v 2 + h + ϕ ) = v ⃗ ⋅ ( v ⃗ × ω ⃗ ) = 0 \vec v\cdot\nabla\left(\frac12v^2+h+\phi\right)=\vec v\cdot\left(\vec v\times\vec\omega\right)=0 v ⋅ ∇ ( 2 1 v 2 + h + ϕ ) = v ⋅ ( v × ω ) = 0 which means $\frac12v^2+h+\phi$ is conserved along each streamline.
No vorticity
∇ × v ⃗ = 0 ⇒ v ⃗ = − ∇ ψ \nabla\times \vec v=0\Rightarrow \vec v=-\nabla\psi ∇ × v = 0 ⇒ v = − ∇ ψ thus
∂ ψ ∂ t + 1 2 v 2 + h + ϕ = F ( t ) \frac{\partial\psi}{\partial t}+\frac12v^2+h+\phi=F(t) ∂ t ∂ ψ + 2 1 v 2 + h + ϕ = F ( t ) and $\nabla F(t)=0$.
Steady & no vorticity
∂ ψ ∂ t = 0 \frac{\partial \psi}{\partial t}=0 ∂ t ∂ ψ = 0 Then $\frac12v^2+h+\phi$ remains a constant everywhere.
Energy Conservation
For fluid dynamics, the energy change in a unit time is given by
D D t ∭ V ( t ) ρ ( e + 1 2 v 2 ) d V = ∭ V ( t ) f ⃗ ⋅ v ⃗ d V + ∯ S ( t ) t ⃗ ⋅ v ⃗ d S \frac{\text D}{\text D t}\iiint_{V(t)}\rho\left(e+\frac12 v^2\right)\text dV=\iiint_{V(t)}\vec f\cdot\vec v\text dV+\oiint_{S(t)}\vec t\cdot\vec v\text dS D t D ∭ V ( t ) ρ ( e + 2 1 v 2 ) d V = ∭ V ( t ) f ⋅ v d V + ∬ S ( t ) t ⋅ v d S where $e$ is the specific internal energy.
If $T^{ij}$ is simply given by $-p\delta^{ij}$, the second term
∯ S ( t ) t ⃗ ⋅ v ⃗ d S = − ∭ V ( t ) ∇ ⋅ ( p v ⃗ ) d V \oiint_{S(t)}\vec t\cdot\vec v\text dS=-\iiint_{V(t)}\nabla\cdot\left(p\vec v\right)\text dV ∬ S ( t ) t ⋅ v d S = − ∭ V ( t ) ∇ ⋅ ( p v ) d V Again with Reynolds transport theorem,
∭ V ( t ) [ ρ D D t ( e + 1 2 v 2 ) + ∇ ⋅ ( p v ⃗ ) ] d V = ∭ V ( t ) f ⃗ ⋅ v ⃗ d V \iiint_{V(t)}\left[\rho\frac{\text D}{\text Dt}\left(e+\frac12 v^2\right)+\nabla\cdot\left(p\vec v\right)\right]\text dV=\iiint_{V(t)}\vec f\cdot \vec v\text dV ∭ V ( t ) [ ρ D t D ( e + 2 1 v 2 ) + ∇ ⋅ ( p v ) ] d V = ∭ V ( t ) f ⋅ v d V ⇒ ρ D D t ( e + 1 2 v 2 ) + ∇ ⋅ ( p v ⃗ ) = f ⃗ ⋅ v ⃗ \Rightarrow \rho\frac{\text D}{\text Dt}\left(e+\frac12 v^2\right)+\nabla\cdot\left(p\vec v\right)=\vec f\cdot \vec v ⇒ ρ D t D ( e + 2 1 v 2 ) + ∇ ⋅ ( p v ) = f ⋅ v ⟺ D D t [ ρ ( e + 1 2 v 2 ) ] − ( e + 1 2 v 2 ) D ρ D t + ∇ ⋅ ( p v ⃗ ) = f ⃗ ⋅ v ⃗ \iff \frac{\text D}{\text Dt}\left[\rho\left( e+\frac12 v^2\right)\right]-\left(e+\frac12 v^2\right)\frac{\text D\rho}{\text Dt}+\nabla\cdot\left(p\vec v\right)=\vec f\cdot \vec v ⟺ D t D [ ρ ( e + 2 1 v 2 ) ] − ( e + 2 1 v 2 ) D t D ρ + ∇ ⋅ ( p v ) = f ⋅ v ⟺ ∂ ∂ t [ ρ ( e + 1 2 v 2 ) ] + ( ∇ ⋅ v ⃗ + v ⃗ ⋅ ∇ ) [ ρ ( e + 1 2 v 2 ) ] + ∇ ⋅ ( p v ⃗ ) = f ⃗ ⋅ v ⃗ \iff \frac{\partial}{\partial t}\left[\rho\left( e+\frac12 v^2\right)\right]+\left(\nabla\cdot \vec v+\vec v\cdot \nabla\right)\left[\rho\left( e+\frac12 v^2\right)\right]+\nabla\cdot\left(p\vec v\right)=\vec f\cdot \vec v ⟺ ∂ t ∂ [ ρ ( e + 2 1 v 2 ) ] + ( ∇ ⋅ v + v ⋅ ∇ ) [ ρ ( e + 2 1 v 2 ) ] + ∇ ⋅ ( p v ) = f ⋅ v ⟺ ∂ ∂ t [ ρ ( e + 1 2 v 2 ) ] + ∇ ⋅ [ ρ v ⃗ ( e + 1 2 v 2 + p ρ ) ] = f ⃗ ⋅ v ⃗ \iff \frac{\partial}{\partial t}\left[\rho\left( e+\frac12 v^2\right)\right]+\nabla\cdot\left[\rho\vec v\left( e+\frac12 v^2+\frac{p}{\rho}\right)\right]=\vec f\cdot \vec v ⟺ ∂ t ∂ [ ρ ( e + 2 1 v 2 ) ] + ∇ ⋅ [ ρ v ( e + 2 1 v 2 + ρ p ) ] = f ⋅ v which is known as the energy equation .
We can integrate and expand the second term so that
∭ V ∇ ⋅ [ ρ v ⃗ ( e + 1 2 v 2 + p ρ ) ] d V = ∯ ρ v n ⃗ ( e + 1 2 v 2 ) d S + ∯ p v n ⃗ d S \iiint_V \nabla\cdot\left[\rho\vec v\left( e+\frac12 v^2+\frac{p}{\rho}\right)\right]\text dV=\oiint\rho v_{\vec n}\left( e+\frac12 v^2\right)\text dS+\oiint pv_{\vec n}\text dS ∭ V ∇ ⋅ [ ρ v ( e + 2 1 v 2 + ρ p ) ] d V = ∬ ρ v n ( e + 2 1 v 2 ) d S + ∬ p v n d S These two terms correspond to the energy flux and the work , respectively.
Here,
h = e + p ρ h=e+\frac p\rho h = e + ρ p is the specific enthalpy.
If $\vec f$ is given by a static potential as
f ⃗ = − ρ ∇ ϕ , ∂ ϕ ∂ t = 0 \vec f=-\rho\nabla\phi,\quad \frac{\partial \phi}{\partial t}=0 f = − ρ ∇ ϕ , ∂ t ∂ ϕ = 0 we can further write the energy equation as
∂ ∂ t [ ρ ( e + 1 2 v 2 + ϕ ) ] + ∇ ⋅ [ ρ v ⃗ ( 1 2 v 2 + h + ϕ ) ] = 0 \frac{\partial}{\partial t}\left[\rho\left( e+\frac12 v^2+\phi\right)\right]+\nabla\cdot\left[\rho\vec v\left(\frac12 v^2+h+\phi\right)\right]=0 ∂ t ∂ [ ρ ( e + 2 1 v 2 + ϕ ) ] + ∇ ⋅ [ ρ v ( 2 1 v 2 + h + ϕ ) ] = 0 Kinetic Energy
We derive the equation of the kinetic energy from Euler equation,
v ⃗ ⋅ [ ∂ v ⃗ ∂ t + ( v ⃗ ⋅ ∇ ) v ⃗ + 1 ρ ∇ p − 1 ρ f ⃗ ] = 0 \vec v\cdot\left[\frac{\partial \vec v}{\partial t}+\left(\vec v\cdot \nabla\right)\vec v+\frac1\rho\nabla p-\frac1\rho\vec f\right]=0 v ⋅ [ ∂ t ∂ v + ( v ⋅ ∇ ) v + ρ 1 ∇ p − ρ 1 f ] = 0 ⟺ ρ ( ∂ ∂ t + v ⃗ ⋅ ∇ ) ( 1 2 v 2 ) + ( v ⃗ ⋅ ∇ ) p = f ⃗ ⋅ v ⃗ \iff \rho\left(\frac{\partial}{\partial t}+\vec v\cdot\nabla\right)\left(\frac12v^2\right)+\left(\vec v\cdot\nabla\right) p=\vec f\cdot\vec v ⟺ ρ ( ∂ t ∂ + v ⋅ ∇ ) ( 2 1 v 2 ) + ( v ⋅ ∇ ) p = f ⋅ v ⟺ ρ D D t ( 1 2 v 2 ) + ( v ⃗ ⋅ ∇ ) p = f ⃗ ⋅ v ⃗ \iff \rho\frac{\text D}{\text Dt}\left(\frac12v^2\right)+\left(\vec v\cdot\nabla\right) p=\vec f\cdot\vec v ⟺ ρ D t D ( 2 1 v 2 ) + ( v ⋅ ∇ ) p = f ⋅ v Thermal Energy (Total - Kinetic)
ρ D D t ( e + 1 2 v 2 ) − ρ D D t ( 1 2 v 2 ) + ∇ ⋅ ( p v ⃗ ) − ( v ⃗ ⋅ ∇ ) p = 0 \rho\frac{\text D}{\text Dt}\left(e+\frac12 v^2\right)-\rho\frac{\text D}{\text Dt}\left(\frac12v^2\right)+\nabla\cdot\left(p\vec v\right)-\left(\vec v\cdot\nabla\right) p=0 ρ D t D ( e + 2 1 v 2 ) − ρ D t D ( 2 1 v 2 ) + ∇ ⋅ ( p v ) − ( v ⋅ ∇ ) p = 0 ⟺ ρ D e D t + ( ∇ ⋅ v ⃗ ) p = 0 \iff \rho\frac{\text De}{\text Dt}+\left(\nabla\cdot\vec v\right) p=0 ⟺ ρ D t D e + ( ∇ ⋅ v ) p = 0 Since
D D t ( 1 ρ ) = 1 ρ ∇ ⋅ v ⃗ \frac{\text D}{\text Dt}\left(\frac1\rho\right)=\frac1\rho\nabla\cdot\vec v D t D ( ρ 1 ) = ρ 1 ∇ ⋅ v We have
D e D t + p D D t ( 1 ρ ) = 0 \frac{\text De}{\text Dt}+p\frac{\text D}{\text Dt}\left(\frac1\rho\right)=0 D t D e + p D t D ( ρ 1 ) = 0 Recall that the first law of thermaldynamics gives
D e D t = − p D D t ( 1 ρ ) + T D s D t \frac{\text De}{\text Dt}=-p\frac{\text D}{\text Dt}\left(\frac1\rho\right)+T\frac{\text Ds}{\text Dt} D t D e = − p D t D ( ρ 1 ) + T D t D s thus
D s D t = 0 \frac{\text Ds}{\text Dt}=0 D t D s = 0 which suggests adiabatic motion of the fluid. In this way, as long as we assume
T i j = − p δ i j T^{ij}=-p\delta^{ij} T ij = − p δ ij we obtain ideal fluid .
Viscosity & Navier-Stokes Equation
In general, the stress tensor is not diagonal,
T i j = − p δ i j + σ i j T^{ij}=-p\delta^{ij}+\sigma^{ij} T ij = − p δ ij + σ ij where $\sigma$ is the viscous tensor .
The $i$-th component of the surface force is thus
F S i = ∯ T i j n j d S = − ∯ p n i d S + ∯ σ i j n j d S F_S^i=\oiint T^{ij}n_j\text dS=-\oiint pn^i\text dS+\oiint\sigma^{ij}n_j\text dS F S i = ∬ T ij n j d S = − ∬ p n i d S + ∬ σ ij n j d S and the force vector is
F ⃗ S = ∭ V ( ∇ ⋅ σ − ∇ p ) d V \vec F_S=\iiint_V\left(\nabla\cdot \sigma-\nabla p\right)\text dV F S = ∭ V ( ∇ ⋅ σ − ∇ p ) d V We claim that $\mathcal T$ is symmetric, that is, $T^{ij}=T^{ji}$. Consider a box with $\delta V=\delta x\delta y\delta z$. In the $z$-direction, the net torque caused leads to the change in angular momentum with respect to the center of $\delta S_{x,y}$
d L z d t = ( r ⃗ × F ⃗ ) z \frac{\text dL_z}{\text dt}=\left(\vec r\times \vec F\right)_z d t d L z = ( r × F ) z ρ δ V ⋅ δ x 2 + δ y 2 6 ω ˙ z ∼ ( T x , y − T y , x ) δ V \rho\delta V\cdot\frac{\delta x^2+\delta y^2}{6}\dot\omega_z\sim\left(T_{x,y}-T_{y,x}\right)\delta V ρ δ V ⋅ 6 δ x 2 + δ y 2 ω ˙ z ∼ ( T x , y − T y , x ) δ V for $\delta x\to0,\delta y\to0$, we have $T{x,y}=T {y,x}$, otherwise the $\dot \omega$ goes to infinite.
Now that we know $\mathcal T$ is symmetric, it has 6 independent components.
Viscous Tensor for Fluid
In a flow with $v_x$ that has a gradient of $\text dv_x/\text dy$, the stress to the $x$-axis is well approximated as
σ x y = η d v x d y \sigma_{xy}=\eta\frac{\text d v_x}{\text dy} σ x y = η d y d v x where $\eta$ is the (dynamic) viscosity determined by microscopic processes of particles. Fluid that satisfies such approximation is known as Nowton fluid .
In astrophysics, a more frequently applied viscosity is the kinematic viscosity $\nu$, defined as
ν = η ρ \nu=\frac{\eta}{\rho} ν = ρ η $\nu$ in a microscopic perspective
Consider a flow with coherent bulk velocity $v_x(y)$ along the $x$-axis, we try to count how much momentum is exchanged at the surface of $y=y_0$ per unit area per unit time, that is, the shear stress.
For simplicity, all we need to consider is a thin layer within $[y-l\text{mfp},y+l \text{mfp}]$, since we assume mometum is fully exchanged in the first collision. In this way, statistically, particles outside this layer cannot transport momentum to the surface $y=y_0$.
The number flux of particles hitting $y=y_0$ from one side can be estimated as
N = 1 3 n c ˉ × 1 2 N=\frac13 n\bar c\times\frac12 N = 3 1 n c ˉ × 2 1 where $n$ is the number density of particles and $\bar c$ is the average velocity dispersion ($\sim c_s$). We introduce a factor of $1/3$ because particles move in all three dimensions and $1/2$ further distinguish motion upwards and downwards. As a result, the shear stress is given by
σ x y = 1 6 n c ˉ m [ v x ( y 0 + l mfp ) − v x ( y 0 − l mfp ) ] = 1 3 ρ c ˉ l mfp ∂ v x ( y 0 ) ∂ x \sigma_{xy}=\frac16n\bar cm\left[v_x(y_0+l_\text{mfp})-v_x(y_0-l_\text{mfp})\right]=\frac13\rho\bar cl_\text{mfp}\frac{\partial v_x(y_0)}{\partial x} σ x y = 6 1 n c ˉ m [ v x ( y 0 + l mfp ) − v x ( y 0 − l mfp ) ] = 3 1 ρ c ˉ l mfp ∂ x ∂ v x ( y 0 ) Thus
ν = 1 3 c ˉ l mfp \nu=\frac13\bar cl_\text{mfp} ν = 3 1 c ˉ l mfp Unfortunately, momentum transition through collisions of particles is not the main cause of viscosity in astrophysics, and this formula is useless... In fact we can estimate the viscous timescale , in which viscosity drives materials to fall into the accretor.
t vis = r 2 ν ∼ r 2 c ˉ l mfp t_\text{vis}=\frac{r^2}{\nu}\sim\frac{r^2}{\bar cl_\text{mfp}} t vis = ν r 2 ∼ c ˉ l mfp r 2 For the supermassive black hole (SMBH), this assumption gives a $t\text{vis}$ if 10 Gyr, which is orders of magtitude longer than the typical AGN timescale ($\sim$ Myr). At present, it is widely agreed that the Magneto-Rotational Instability (MRI) causes the viscosity in accretion disks. In this case, the $l \text{mfp}$ is given by the thickness of the disk $H\sim0.01 r$, and the viscous timescale is quite satistying.
Construction of Stress Tensor
$\sigma{ij}$ is symmetric since $T {ij}$ is symmetric, but the original definition.
σ i j = η ∂ v i ∂ x j \sigma_{ij}=\eta\frac{\partial v_i}{\partial x^j} σ ij = η ∂ x j ∂ v i is not. In this way, we expand it into symmetric and anti-symmetric parts, as we can do to any second-order tensor
σ i j = η 2 ( ∂ v i ∂ x j + ∂ v j ∂ x i ) + η 2 ( ∂ v i ∂ x j − ∂ v j ∂ x i ) \sigma_{ij}=\frac{\eta}{2}\left(\frac{\partial v_i}{\partial x^j}+\frac{\partial v_j}{\partial x^i}\right)+\frac{\eta}{2}\left(\frac{\partial v_i}{\partial x^j}-\frac{\partial v_j}{\partial x^i}\right) σ ij = 2 η ( ∂ x j ∂ v i + ∂ x i ∂ v j ) + 2 η ( ∂ x j ∂ v i − ∂ x i ∂ v j ) The second (anti-symmetric) term corresponds to the rotation. It has no viscosity, so we will consider a linear relation between $\sigma{ij}$ and $E {ij}$, which is defined as
E i j ≡ 1 2 ( ∂ v i ∂ x j + ∂ v j ∂ x i ) E_{ij}\equiv\frac{1}{2}\left(\frac{\partial v_i}{\partial x^j}+\frac{\partial v_j}{\partial x^i}\right) E ij ≡ 2 1 ( ∂ x j ∂ v i + ∂ x i ∂ v j ) Fluid doesn't have a preferred direction.
Intuitionally, we write down the (multi-)linear relation as
σ i j = α i j k l E k l \sigma_{ij}=\alpha_{ijkl}E^{kl} σ ij = α ijk l E k l We don't have any information of $\alpha$, but since fluid is isotropic , we would like to believe the tensor itself is isotropic, which means it has the same components in all rotated coordinate systems. There are only three basic 4-order isotropic tensors, all of which are the combination of the famous Kronecker delta $\delta_{ij}$.
A i j k l = δ i j δ k j , δ i k δ j l , δ i l δ j k A_{ijkl}=\delta_{ij}\delta_{kj},\ \delta_{ik}\delta_{jl},\ \delta_{il}\delta_{jk} A ijk l = δ ij δ kj , δ ik δ j l , δ i l δ jk So we assume $\alpha_{ijkl}$ is simply their linear combination
σ i j = ( B δ i j δ k j + C δ i k δ j l + D δ i l δ j k ) E k l = B δ i j E k k + C E i j + D E j i = B ( ∇ ⋅ v ⃗ ) δ i j + ( C + D ) E i j \begin{align*}
\sigma_{ij}&=\left(B\delta_{ij}\delta_{kj}+C\delta_{ik}\delta_{jl}+D\delta_{il}\delta_{jk}\right)E^{kl}\\
&=B\delta_{ij}E^k_k+CE_{ij}+DE_{ji}\\
&=B\left(\nabla\cdot \vec v\right)\delta_{ij}+(C+D)E_{ij}
\end{align*} σ ij = ( B δ ij δ kj + C δ ik δ j l + D δ i l δ jk ) E k l = B δ ij E k k + C E ij + D E ji = B ( ∇ ⋅ v ) δ ij + ( C + D ) E ij Decompose $\sigma_{ij}$ into trace part and traceless part
tr ( σ ) = δ i j σ i j = ( 3 B + C + D ) ∇ ⋅ v ⃗ \text{tr}\left(\sigma\right)=\delta^{ij}\sigma_{ij}=(3B+C+D)\nabla\cdot \vec v tr ( σ ) = δ ij σ ij = ( 3 B + C + D ) ∇ ⋅ v Thus
σ i j tr = 1 3 ( 3 B + C + D ) ( ∇ ⋅ v ⃗ ) δ i j \sigma_{ij}^{\text{tr}}=\frac{1}3(3B+C+D)\left(\nabla\cdot \vec v\right)\delta_{ij} σ ij tr = 3 1 ( 3 B + C + D ) ( ∇ ⋅ v ) δ ij σ i j trf = − C + D 3 ( ∇ ⋅ v ⃗ ) δ i j + C + D 2 ( ∂ v i ∂ x j + ∂ v j ∂ x i ) \sigma_{ij}^\text{trf}=-\frac{C+D}3\left(\nabla\cdot \vec v\right)\delta_{ij}+\frac{C+D}2\left(\frac{\partial v_i}{\partial x^j}+\frac{\partial v_j}{\partial x^i}\right) σ ij trf = − 3 C + D ( ∇ ⋅ v ) δ ij + 2 C + D ( ∂ x j ∂ v i + ∂ x i ∂ v j ) Redefine $\eta$ and $\zeta$ as
η ≡ 1 2 ( C + D ) , ζ ≡ 1 3 ( 3 B + C + D ) \eta\equiv \frac12(C+D),\quad\zeta\equiv\frac13(3B+C+D) η ≡ 2 1 ( C + D ) , ζ ≡ 3 1 ( 3 B + C + D ) Finally, with only two free parameters, we may obtain
σ i j = η [ ∂ v i ∂ x j + ∂ v j ∂ x i − 2 3 ( ∇ ⋅ v ⃗ ) δ i j ] + ζ ( ∇ ⋅ v ⃗ ) δ i j \sigma_{ij}=\eta\left[\frac{\partial v_i}{\partial x^j}+\frac{\partial v_j}{\partial x^i}-\frac23\left(\nabla\cdot \vec v\right)\delta_{ij}\right]+\zeta\left(\nabla\cdot \vec v\right)\delta_{ij} σ ij = η [ ∂ x j ∂ v i + ∂ x i ∂ v j − 3 2 ( ∇ ⋅ v ) δ ij ] + ζ ( ∇ ⋅ v ) δ ij First term - traceless, pure shear viscosity
Second term - bulk viscosity
When Stokes first derived this formula, he believed the bulk viscosity vanished. In fact this term is usually close to 0.
Navier-Stokes Equation
Right now,
T i j = − p δ i j + σ i j T^{ij}=-p\delta^{ij}+\sigma^{ij} T ij = − p δ ij + σ ij and we can rewrite the momentum tensor by adding the viscous tensor
Π i j = ρ v i v j + δ i j p − σ i j \Pi^{ij}=\rho v^iv^j+\delta^{ij}p-\sigma^{ij} Π ij = ρ v i v j + δ ij p − σ ij From this correction, the conservation of momentum (Euler equation) gives
D v ⃗ D t = 1 ρ f ⃗ + 1 ρ ∇ ⋅ T = − 1 ρ ∇ p + 1 ρ f ⃗ + 1 ρ ∇ ⋅ σ \frac{\text D\vec v}{\text Dt}=\frac1\rho\vec f+\frac1\rho\nabla\cdot\mathcal T=-\frac1\rho\nabla p+\frac1\rho\vec f+\frac1\rho\nabla\cdot\sigma D t D v = ρ 1 f + ρ 1 ∇ ⋅ T = − ρ 1 ∇ p + ρ 1 f + ρ 1 ∇ ⋅ σ Since
( ∇ ⋅ σ ) i = η ∂ j [ ∂ j v i + ∂ i v j − 2 3 ∂ k v k δ i j ] + ζ ∂ j ∂ k v k δ i j = η [ ∇ 2 v i + 1 3 ∂ i ( ∇ ⋅ v ⃗ ) ] + ζ ∂ i ( ∇ ⋅ v ⃗ ) \begin{align*}
(\nabla\cdot\sigma)_i&=\eta\partial^j\left[\partial_jv_i+\partial_iv_j-\frac23\partial_kv^k\delta_{ij}\right]+\zeta\partial^j\partial_kv^k\delta_{ij}\\
&=\eta\left[\nabla^2v_i+\frac13\partial_i(\nabla\cdot\vec v)\right]+\zeta\partial_i(\nabla\cdot\vec v)
\end{align*} ( ∇ ⋅ σ ) i = η ∂ j [ ∂ j v i + ∂ i v j − 3 2 ∂ k v k δ ij ] + ζ ∂ j ∂ k v k δ ij = η [ ∇ 2 v i + 3 1 ∂ i ( ∇ ⋅ v ) ] + ζ ∂ i ( ∇ ⋅ v ) ⟺ ∇ ⋅ σ = η ∇ 2 v ⃗ + ( ζ + 1 3 η ) ∇ ( ∇ ⋅ v ⃗ ) \iff \nabla\cdot\sigma=\eta\nabla^2\vec v+\left(\zeta+\frac13\eta\right)\nabla\left(\nabla\cdot\vec v\right) ⟺ ∇ ⋅ σ = η ∇ 2 v + ( ζ + 3 1 η ) ∇ ( ∇ ⋅ v ) Finally we derive the Navier-Stokes equation
∂ v ⃗ ∂ t + ( v ⃗ ⋅ ∇ ) v ⃗ = − 1 ρ ∇ p + 1 ρ f ⃗ + 1 ρ η ∇ 2 v ⃗ + 1 ρ ( ζ + 1 3 η ) ∇ ( ∇ ⋅ v ⃗ ) \frac{\partial \vec v}{\partial t}+\left(\vec v\cdot \nabla\right)\vec v=-\frac1\rho\nabla p+\frac1\rho\vec f+\frac1\rho\eta\nabla^2\vec v+\frac1\rho\left(\zeta+\frac13\eta\right)\nabla\left(\nabla\cdot\vec v\right) ∂ t ∂ v + ( v ⋅ ∇ ) v = − ρ 1 ∇ p + ρ 1 f + ρ 1 η ∇ 2 v + ρ 1 ( ζ + 3 1 η ) ∇ ( ∇ ⋅ v ) In the end, let us consider the energy conservation in the viscous situation
D D t ∭ V ( t ) ρ ( e + 1 2 v 2 ) d V = ∭ V ( t ) f ⃗ ⋅ v ⃗ d V + ∯ S ( t ) t ⃗ ⋅ v ⃗ d S \frac{\text D}{\text D t}\iiint_{V(t)}\rho\left(e+\frac12 v^2\right)\text dV=\iiint_{V(t)}\vec f\cdot\vec v\text dV+\oiint_{S(t)}\vec t\cdot\vec v\text dS D t D ∭ V ( t ) ρ ( e + 2 1 v 2 ) d V = ∭ V ( t ) f ⋅ v d V + ∬ S ( t ) t ⋅ v d S ⟺ ∭ V ( t ) ρ D D t ( e + 1 2 v 2 ) d V = ∭ V ( t ) f ⃗ ⋅ v ⃗ d V + ∭ V ( t ) ∇ ⋅ ( T ⋅ v ⃗ ) d V \iff \iiint_{V(t)}\rho\frac{\text D}{\text D t}\left(e+\frac12 v^2\right)\text dV=\iiint_{V(t)}\vec f\cdot\vec v\text dV+\iiint_{V(t)}\nabla\cdot\left(\mathcal T\cdot\vec v\right)\text dV ⟺ ∭ V ( t ) ρ D t D ( e + 2 1 v 2 ) d V = ∭ V ( t ) f ⋅ v d V + ∭ V ( t ) ∇ ⋅ ( T ⋅ v ) d V ⟺ ρ D D t ( e + 1 2 v 2 ) = f ⃗ ⋅ v ⃗ + ∇ ⋅ ( T ⋅ v ⃗ ) \iff\rho\frac{\text D}{\text D t}\left(e+\frac12 v^2\right)=\vec f\cdot\vec v+\nabla\cdot\left(\mathcal T\cdot\vec v\right) ⟺ ρ D t D ( e + 2 1 v 2 ) = f ⋅ v + ∇ ⋅ ( T ⋅ v ) By multiplying N-S equation with $\rho\vec v$, we have
ρ D D t ( 1 2 v 2 ) = f ⃗ ⋅ v ⃗ + v ⃗ ⋅ ( ∇ ⋅ T ) \rho\frac{\text D}{\text Dt}\left(\frac12v^2\right)=\vec f\cdot\vec v+\vec v\cdot\left(\nabla\cdot\mathcal T\right) ρ D t D ( 2 1 v 2 ) = f ⋅ v + v ⋅ ( ∇ ⋅ T ) ⇒ ρ D e D t = ( T ⋅ ∇ ) ⋅ v ⃗ \Rightarrow \rho\frac{\text De}{\text D t}=\left(\mathcal T\cdot\nabla\right)\cdot \vec v ⇒ ρ D t D e = ( T ⋅ ∇ ) ⋅ v Now we expand the RHS
( T ⋅ ∇ ) ⋅ v ⃗ = ( ∂ i v j ) T i j = − ∂ i v i p + ( ∂ i v j ) σ i j = − p ( ∇ ⋅ v ⃗ ) + η ( ∂ i v j ) [ 2 E i j − 2 3 ( ∇ ⋅ v ⃗ ) δ i j ] + ζ ( ∂ i v i ) ( ∇ ⋅ v ⃗ ) = − p ( ∇ ⋅ v ⃗ ) + η ( ∂ i v j ) ( ∂ j v i + ∂ i v j ) + ( ζ − 2 3 η ) ( ∇ ⋅ v ⃗ ) 2 \begin{align*}
\left(\mathcal T\cdot\nabla\right)\cdot \vec v&=(\partial_iv_j)T^{ij}\\
&=-\partial_iv^ip+(\partial_iv_j)\sigma^{ij}\\
&=-p(\nabla\cdot\vec v)+\eta(\partial_iv_j)\left[2E^{ij}-\frac23\left(\nabla\cdot \vec v\right)\delta_{ij}\right]+\zeta(\partial_iv^i)\left(\nabla\cdot \vec v\right)\\
&=-p(\nabla\cdot\vec v)+\eta(\partial_iv_j)(\partial^jv^i+\partial^iv^j)+\left(\zeta-\frac23\eta\right)\left(\nabla\cdot \vec v\right)^2
\end{align*} ( T ⋅ ∇ ) ⋅ v = ( ∂ i v j ) T ij = − ∂ i v i p + ( ∂ i v j ) σ ij = − p ( ∇ ⋅ v ) + η ( ∂ i v j ) [ 2 E ij − 3 2 ( ∇ ⋅ v ) δ ij ] + ζ ( ∂ i v i ) ( ∇ ⋅ v ) = − p ( ∇ ⋅ v ) + η ( ∂ i v j ) ( ∂ j v i + ∂ i v j ) + ( ζ − 3 2 η ) ( ∇ ⋅ v ) 2 where
( ∂ i v j ) ( ∂ j v i + ∂ i v j ) = 1 2 ( ∂ i v j + ∂ j v i ) ( ∂ j v i + ∂ i v j ) − 1 2 ( ∂ i v j − ∂ j v i ) ( ∂ j v i + ∂ i v j ) = 2 E i j E i j \begin{align*}
(\partial_iv_j)(\partial^jv^i+\partial^iv^j)&=\frac12(\partial_iv_j+\partial_jv_i)(\partial^jv^i+\partial^iv^j)\\
&-\frac12(\partial_iv_j-\partial_jv_i)(\partial^jv^i+\partial^iv^j)\\
&=2E_{ij}E^{ij}
\end{align*} ( ∂ i v j ) ( ∂ j v i + ∂ i v j ) = 2 1 ( ∂ i v j + ∂ j v i ) ( ∂ j v i + ∂ i v j ) − 2 1 ( ∂ i v j − ∂ j v i ) ( ∂ j v i + ∂ i v j ) = 2 E ij E ij as the second term obviously vanishes. Finally,
( T ⋅ ∇ ) ⋅ v ⃗ = − p ( ∇ ⋅ v ⃗ ) + 2 η E i j E i j + ( ζ − 2 3 η ) ( ∇ ⋅ v ⃗ ) 2 \left(\mathcal T\cdot\nabla\right)\cdot \vec v=-p(\nabla\cdot\vec v)+2\eta E_{ij}E^{ij}+\left(\zeta-\frac23\eta\right)\left(\nabla\cdot \vec v\right)^2 ( T ⋅ ∇ ) ⋅ v = − p ( ∇ ⋅ v ) + 2 η E ij E ij + ( ζ − 3 2 η ) ( ∇ ⋅ v ) 2 Define
Φ vis = 2 η E i j E i j + ( ζ − 2 3 η ) ( ∇ ⋅ v ⃗ ) 2 \Phi_\text{vis}=2\eta E_{ij}E^{ij}+\left(\zeta-\frac23\eta\right)\left(\nabla\cdot \vec v\right)^2 Φ vis = 2 η E ij E ij + ( ζ − 3 2 η ) ( ∇ ⋅ v ) 2 we have
ρ [ D e D t + p D D t ( 1 ρ ) ] = Φ vis = ρ T D s D t \rho\left[\frac{\text De}{\text D t}+p\frac{\text D}{\text Dt}\left(\frac1\rho\right)\right]=\Phi_\text{vis}=\rho T\frac{\text Ds}{\text Dt} ρ [ D t D e + p D t D ( ρ 1 ) ] = Φ vis = ρT D t D s Thus,
D S D t = D D t ∭ V ρ s d V = ∭ V ρ D s D t d V = ∭ V Φ vis T d V \frac{\text DS}{\text Dt}=\frac{\text D}{\text Dt}\iiint_{V}\rho s\text dV=\iiint_{V}\rho \frac{\text Ds}{\text Dt}\text dV=\iiint_{V}\frac{\Phi_\text{vis}}{T}\text dV D t D S = D t D ∭ V ρ s d V = ∭ V ρ D t D s d V = ∭ V T Φ vis d V The second law of thermaldynamics tells us
D S D t ≥ 0 \frac{\text DS}{\text Dt}\ge0 D t D S ≥ 0 In fact, we can prove
Φ vis ≥ 0 \Phi_\text{vis}\ge0 Φ vis ≥ 0 The Reynolds Number
The Reynolds number is defined as
Re = U L ν \text{Re}=\frac{UL}{\nu} Re = ν UL where $U$ is the flow speed, $L$ is a characteristic linear dimension, and $\nu$ is the kinematic viscosity. It is naturally connected with the Navier-Stokes equation when we try to derive a dimensionless form of the equation.
By assuming incompressible fluid ($\nabla\cdot \vec v=0$) and neglecting the external force, the Navier-Stokes equation goes like
∂ v ⃗ ∂ t + ( v ⃗ ⋅ ∇ ) v ⃗ = − 1 ρ ∇ p + ν ∇ 2 v ⃗ \frac{\partial \vec v}{\partial t}+\left(\vec v\cdot \nabla\right)\vec v=-\frac1\rho\nabla p+\nu\nabla^2\vec v ∂ t ∂ v + ( v ⋅ ∇ ) v = − ρ 1 ∇ p + ν ∇ 2 v Define several dimensionless quantities,
r ′ = r L , v ⃗ ′ = v ⃗ U , t ′ = t L / u , ρ ′ = ρ ρ 0 , p ′ = p ρ 0 v 2 r'=\frac rL,\quad \vec v'=\frac{\vec v}{U},\quad t'=\frac{t}{L/u},\quad \rho'=\frac{\rho}{\rho_0},\quad p'=\frac{p}{\rho_0v^2} r ′ = L r , v ′ = U v , t ′ = L / u t , ρ ′ = ρ 0 ρ , p ′ = ρ 0 v 2 p We have
∂ v ⃗ ′ ∂ t ′ + ( v ⃗ ′ ⋅ ∇ ′ ) v ⃗ ′ = − 1 ρ ′ ∇ ′ p ′ + ν U L ∇ ′ 2 v ⃗ ′ ≡ − 1 ρ ′ ∇ ′ p ′ + 1 Re ∇ ′ 2 v ⃗ ′ \frac{\partial \vec v'}{\partial t'}+\left(\vec v'\cdot \nabla'\right)\vec v'=-\frac1\rho'\nabla' p'+\frac\nu{UL}\nabla'^2\vec v'\equiv -\frac1\rho'\nabla' p'+\frac1{\text{Re}}\nabla'^2\vec v' ∂ t ′ ∂ v ′ + ( v ′ ⋅ ∇ ′ ) v ′ = − ρ 1 ′ ∇ ′ p ′ + UL ν ∇ ′2 v ′ ≡ − ρ 1 ′ ∇ ′ p ′ + Re 1 ∇ ′2 v ′ The Reynolds number is simply there.
If $\text{Re}\to\infty$, the equation goes back to the Euler equation (ideal fluid).
Small $\text{Re}$ corresponds to laminar flow , a fluid flowing in parallel layers with no disruption between the layers, while large $\text{Re}$ introduces annoying turbulence , a fluid motion characterized by chaotic changes in pressure and flow velocity. In astrophysics, $\text{Re}$ is usually large.
Notes
Dimensionless Navier-Stokes equation does not depend on $\rho_0$ - scale-free . Once a hydrodynamical simulation is done with certain density, it is done with any other density. This beautiful scale-free law breaks when
$P-\rho$ relation is not a simple power-law, and has certain breaks.
the simulation has to take opacity/cooling/heating into account, none of which is scale-free.
If we introduce gravity as the external force
f ⃗ ext = − G M r 2 r ⃗ ^ \vec f_\text{ext}=-\frac{GM}{r^2}\hat{\vec r} f ext = − r 2 GM r ^ the corresponding term in the dimensionless Naiver-Stokes equation is
− 1 L U 2 G M r ′ 2 -\frac{1}{LU^2}\frac{GM}{r'^2} − L U 2 1 r ′2 GM It seems to have dependence on $M$, however, if we set
U = c , L = r g ≡ G M c 2 U=c,\quad L=r_\text{g}\equiv\frac{GM}{c^2} U = c , L = r g ≡ c 2 GM where $r_\text g$ is the Schiwarzchild radius. This term is simply
Around a black hole, Navier-Stokes equation is again scale-free, this time for the mass of the accretor.
