User:Egm6936.f09/When to use incompressible flow model?

Fluid motion at steady state is considered incompressible when Mach number is much less than unity, given in Eq.(24). In fluid flow analysis, including turbulent flow, it is very frequent to take the fluid, either a gas or a liquid, to be incompressible because this assumption simplifies the governing equations. For example, McComb (1992, p.2) states that the fluid motion is considered incompressible throughout the book and refers to the book by Batchelor (1967, p.167). In practice, when the flow speed is less than half of the sound speed, incompressibility condition holds (Bolster et al., 2011).

 = Introduction: What is incompressibility? =

The velocity and the pressure fields, $$\displaystyle \mathbf v $$ and $$\displaystyle p $$ respectively, are obtained from a solution to the momentum balance and mass balance equations. The latter, also known as the continuity equation, is expressed as


 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial \rho} {\partial t} + \rm div (\rho \mathbf v) = 0 $$  (1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

or, in terms of material time derivative,


 * {| style="width:100%" border="0"

$$  \displaystyle \frac {1}{\rho} \frac {D \rho}{D t} + \rm div (\mathbf v) = 0 $$  (2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

where $$\displaystyle \rho $$ is the density of the fluid, which is found from an equation of state of the form $$\displaystyle \rho = \rho(p,T) $$ where $$\displaystyle T $$ is the temperature. A fluid (gas or liquid) is said to be incompressible if its volume (or density) is constant upon an applied external force, e.g., pressure. Then, for an incompressible fluid, the continuity equation, i.e. Eq.(1) simplifies considerably and is expressed as


 * {| style="width:100%" border="0"

$$  \displaystyle \rm div (\mathbf v) = 0 $$     (3)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * 
 * }

Liquids are usually assumed to be incompressible because the effect of the pressure on its density is negligible (see Schlichting 1968, p.10 ). In an isothermal system, even though the density of the gases is pressure dependent, the question is whether this change is important enough to take into consideration. Therefore, the goal in this section is to derive the criteria under which a fluid can be considered incompressible. We will follow three approaches presented by Malvern (1969, p.443), by Landau and Lifshitz 1975, p.245 , and by Batchelor 1967, p.167. Even though the starting points and the steps they take are different, they come to the same conclusion.

= A quick derivation using gas dynamical equation =

Here, we obtain the condition for using incompressible flow model starting with the gas dynamical equation. By assuming that the flow velocity is very small compared to the speed of sound, and that the gradient of the velocity is bounded by some finite number, then the divergence of the velocity field is negligible, i.e., approaching the incompressibility condition; Malvern 1969, p.443

The gas dynamical equation


 * {| style="width:100%" border="0"

$$  \displaystyle c^2 {\rm div \,} \mathbf u - \mathbf u \cdot {\rm grad \,} \mathbf u \cdot \mathbf u = 0 $$     (4)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

or


 * {| style="width:100%" border="0"

$$  \displaystyle {\rm div \,} \mathbf u = \frac{\mathbf u \cdot {\rm grad \,} \mathbf u \cdot \mathbf u}{c^2} $$     (5)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

thus


 * {| style="width:100%" border="0"

$$  \displaystyle \parallel \mathbf u \parallel \ll c {\ \rm and \ } \parallel {\rm grad \,} \mathbf u \parallel < M \Rightarrow {\rm div \,} \mathbf u \approx 0 $$     (6)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

for some finite number $$\displaystyle M$$.

The downside of this approach is that the condition on the velocity field (being less than the speed of sound) had to be assumed ahead. Instead, Landau and Lifshitz 1975 and by Batchelor 1967 took different starting points to arrive at the velocity condition for the validity of using incompressible-flow model.

Appendix: Derivation of the gas dynamical equation
The steady state Euler equation, dotted with the velocity vector is given as


 * {| style="width:100%" border="0"

$$  \displaystyle \left(\rho (\rm grad \, \mathbf v) \cdot \mathbf v \right) \cdot \mathbf v = \left(-\rm grad p\right) \cdot \mathbf v $$ (7)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

Using the definition of the material time derivative, the right hand side of Eq.(7) can be written as


 * {| style="width:100%" border="0"

$$  \displaystyle -\rm grad p \cdot \mathbf v = -\frac{Dp}{Dt} + \frac{\partial p} {\partial t} $$ (8)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

As we are concerned with steady state flow, the second term in Eq.(8) vanishes and the first term can be expanded as in Eq.(28) where $$\displaystyle DS/Dt=0$$ as isentropic (constant entropy, i.e., reversible adiabatic) process is assumed. Hence, Eq.(8) reduces to


 * {| style="width:100%" border="0"

$$  \displaystyle -\rm grad p \cdot \mathbf v = - c^2 \frac {D \rho}{D t} = -c^2 \left(-\rho \, \nabla \cdot \mathbf v \right) $$  (9)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

For the second inequality in Eq.(9), the continuity equation, Eq.(2) is used. Combining Eq.(9) with Eq.(7) yields the gas dynamical equation.

= Derivation of the incompressibility conditions using sound wave equations =

Here, we derive the incompressibility conditions from the sound wave equation (see Landau and Lifshitz 1975, p.245 ). First, we start from the momentum equation, derive the wave equation, and obtain the incompressibility conditions.

A sound wave is an oscillatory motion of small amplitude in a compressible fluid. The pressure and the density may be written as

<span id="(10)">
 * {| style="width:100%" border="0"

$$  \displaystyle p = p_0 + p^ \prime $$  (10)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and

<span id="(11)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho = \rho_0 + \rho^ \prime $$  (11)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

where $$\displaystyle p_0 $$ and $$\displaystyle  \rho_0 $$ represent the undisturbed/equilibrium pressure and density and $$\displaystyle  p^ \prime $$ and $$\displaystyle  \rho^ \prime $$ are the small changes/disturbances to the equilibrium quantities. The continuity equation, given in Eq.(1) becomes

<span id="(12)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial \rho^ \prime } {\partial t} + \rho_0 \rm div \mathbf v = 0 $$  (12)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Navier Stokes equation, which is

<span id="(13)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial \mathbf v} {\partial t} + (\rm grad \, \mathbf v) \cdot \mathbf v = -\frac{1}{\rho} \rm grad p + \nu \nabla^2 \mathbf v $$ (13)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For an inviscid fluid, Eq.(13) becomes Euler’s equation, which is

<span id="(14)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial \mathbf v} {\partial t} + (\rm grad \, \mathbf v) \cdot \mathbf v = -\frac{1}{\rho} \rm grad p $$ (14)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Omitting higher order terms, Eq.(14) reduces to

<span id="(15)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial \mathbf v} {\partial t} + \frac{1}{\rho_0} \rm grad p^ \prime = 0 $$  (15)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For an ideal fluid (deformation of fluid elements is an isentropic process, i.e. adiabatic and reversible), the small changes in pressure and density are related to each other as

<span id="(16)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {p^ \prime}{\rho^ \prime} = \left( \frac {\partial p} {\partial \rho_0} \right)_S $$  (16)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Then, the continuity equation, i.e. Eq.(12) becomes

<span id="(17)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial p^ \prime } {\partial t} + \rho_0 \left( \frac {\partial p} {\partial \rho_0} \right)_S \rm div \mathbf v = 0 $$  (17)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For an ideal fluid, the velocity field may be expressed as a potential $$\displaystyle \mathbf v = \rm grad \phi $$. Then, Eq.(15) becomes

<span id="(18)">
 * {| style="width:100%" border="0"

$$  \displaystyle p^ \prime = -\rho_0 \frac {\partial \phi} {\partial t} $$ (18)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

Continuity, i.e. Eq.(17) becomes <span id="(19)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial^2 \phi } {\partial t^2} - c^2 \rm div (\rm grad \phi) = 0 $$  (19) which is called a wave equation. Here, $$\displaystyle c := \sqrt{(\partial p/\partial \rho_0)_S} $$. Now, let’s consider a wave equation in one dimension as
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(20)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial^2 \phi} {\partial x^2} - \frac {1} {c^2} \frac {\partial^2 \phi} {\partial t^2} = 0 $$  (20) which has a solution of the form <span id="(21)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \phi = f_1(x-ct) + f_2(x+ct) $$  (21) This form of the solution is also valid for $$\displaystyle (p^\prime, \rho^\prime, \mathbf v) $$. The function $$\displaystyle f_1 $$ implies that the value of $$\displaystyle \phi $$ at $$\displaystyle t=0 $$ is also observed after time $$\displaystyle t $$ at a distance $$\displaystyle ct $$ from the previous point. Therefore $$\displaystyle c $$ represents the speed of the sound and the function $$\displaystyle f_1 $$ is called a travelling plane wave. Note that $$\displaystyle f_2 $$ propagates in the opposite direction. For a travelling wave, as $$\displaystyle \phi = f(x-ct) $$, $$\displaystyle v_x = v = f^\prime (x-ct) $$ and $$\displaystyle p^\prime = -\rho_0 \partial \phi/ \partial t = \rho_0 c f^\prime(x-ct) $$. Therefore, <span id="(22)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle v = p^\prime/(\rho_0 c) $$ (22) Substituting Eq.(16) into Eq.(22), with the definition of $$\displaystyle c$$ yields <span id="(23)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{\rho^\prime}{\rho_0} = \frac{v}{c} $$  (23) The fluid may be treated as incompressible if $$\displaystyle \rho^\prime/\rho_0 \ll 1 $$ or
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(24)">
 * {| style="width:100%" border="0"

$$  \displaystyle v/c \ll 1 $$     (24)
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

This last quantity is also known as Mach number. At $$\displaystyle 25^{\circ}{\rm C} $$, the speed of sound in air is $$\displaystyle 346.1 \rm {m/s} $$ and it is $$\displaystyle 1497\rm {m/s} $$  in water. Therefore, for all practical purposes, the liquids are assumed to be incompressible and for gases, one needs to check the Mach number.

For an unsteady flow, a second criteria needs to be satisfied

<span id="(25)">
 * {| style="width:100%" border="0"

$$  \displaystyle \tau \gg l/c $$     (25) where $$\displaystyle \tau $$  and $$\displaystyle l$$  are the time and length scales over which the fluid goes significant changes. The above equation is obtained when the $$\displaystyle \partial \mathbf v/\partial t$$ and $$\displaystyle (1/\rho) \rm grad p$$ in Euler’s equation are comparable. The order of magnitude of each term are $$\displaystyle v/\tau \sim \Delta p/(l \rho)$$. Then, $$\displaystyle \Delta \rho \sim l\rho v /(\tau c^2)$$. When the $$\displaystyle \partial \rho / \partial t$$ and $$\displaystyle \rho \rm{div} \mathbf v$$ terms in the continuity equation are compared, $$\displaystyle \partial \rho / \partial t$$ can be neglected if $$\displaystyle \Delta \rho/\tau \ll \rho v/l$$, which yields Eq.(25). Note that, for steady flow, $$\displaystyle \tau = \infty$$, and the second criteria, Eq.(25) is always fulfilled.
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

= Alternative derivation =

Here, we follow primarily the exposition in Batchelor (1967, p.167). Let the characteristic length and velocity scales be $$\displaystyle L$$ and $$\displaystyle U$$, respectively. Then, the spatial derivative of the velocity field has an order of magnitude of $$\displaystyle U/L$$. The velocity field can be approximated as solenoidal if

<span id="(26)">
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{vmatrix} \nabla \cdot \mathbf v \end{vmatrix} \ll \frac{U}{L} $$  (26)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

i.e., using Eq.(2)

<span id="(27)">
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{vmatrix} \displaystyle \frac {1}{\rho} \frac {D \rho}{D t} \end{vmatrix} \ll \frac{U}{L} $$  (27)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

For a homogeneous fluid, the degrees of freedom $$\displaystyle F$$, given as $$ \displaystyle F=C-P+2$$ is two. Here, $$\displaystyle C $$ and $$\displaystyle P$$ are the number of components and phases, respectively and they are both one for a homogenous fluid. We choose the density $$\displaystyle \rho$$ and the entropy per unit mass $$\displaystyle S$$ as two independent variables. Then,

<span id="(28)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D p}{D t} = c^2 \frac {D \rho}{D t} + \left( \frac {\partial p} {\partial S} \right)_\rho \frac {D S}{Dt}
 * style="width:95%" |
 * style="width:95%" |

$$  (28) where $$\displaystyle c := \sqrt{(\partial p/\partial \rho)_S} $$. Then the incompressibility condition Eq.(27) becomes
 * <p style="text-align:right">
 * }

<span id="(29)">
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{vmatrix} \displaystyle \frac {1}{\rho c^2} \frac {D p}{D t} - \frac {1}{\rho c^2} \left( \frac {\partial p} {\partial S} \right)_\rho \frac {D S}{Dt} \end{vmatrix}
 * style="width:95%" |
 * style="width:95%" |

\ll \frac{U}{L} $$  (29) Using the definitions for the material time derivative, i.e.
 * <p style="text-align:right">
 * }

<span id="(30)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {D p}  {Dt} =  \frac {\partial p}  {\partial t} + {\rm grad} \, p  \cdot \mathbf v $$ (30)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

and the Euler equation with the body force term, which is

<span id="(31)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {\partial \mathbf v} {\partial t} + (\rm grad \, \mathbf v) \cdot \mathbf v = -\frac{1}{\rho} \rm grad p + \rho \mathbf g $$ (31) eq.(29) becomes
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(32)">
 * {| style="width:100%" border="0"

$$  \displaystyle \begin{vmatrix} \displaystyle \frac {1}{\rho c^2} \frac {\partial p}{\partial t} - \frac {1}{2 c^2} \frac {D u^2}{D t} + \frac {1}{c^2} \left( \mathbf u \cdot \mathbf g \right) - \frac {1}{\rho c^2} \left( \frac {\partial p} {\partial S} \right)_\rho \frac {D S}{Dt} \end{vmatrix}
 * style="width:95%" |
 * style="width:95%" |

\ll \frac{U}{L}

$$  (32)
 * <p style="text-align:right">
 * }

The most general way to satisfy the inequality in Eq.(32) is that each term has a magnitude much smaller than $$\displaystyle U/L$$. The second term in Eq.(32) has the same magnitude as $$\displaystyle {\rm grad} \, v^2 \cdot   \mathbf v $$, which is $$\displaystyle U^3/L$$. Then, the first criteria for the incompressibility is

<span id="(33)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac{U^2}{c^2} \ll 1 $$     (33) which is the condition in Eq.(25). The pressure term in Eq.(32) has a magnitude of rate of momentum per area, i.e., $$\displaystyle \rho L U/ \tau $$. Then, the first term yields
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

<span id="(34)">
 * {| style="width:100%" border="0"

$$  \displaystyle \tau^2 \gg l^2/c^2 $$     (34) which is the condition in Eq.(D17). Note that if $$\displaystyle \tau=L/U $$, then the condition in Eq.(34) becomes same as in Eq.(33). If $$\displaystyle \tau=L/c $$, then $$\displaystyle L $$ is the wavelength of a sound wave. The gravitational term in Eq.(32) yields
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * <p style="text-align:right">
 * }

<span id="(35)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {gL}{c^2} \ll 1
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$     (35)
 * <p style="text-align:right">
 * }

The term with the rate of entropy change in Eq.(32) can be decomposed into two parts. The first one

<span id="(36)">
 * {| style="width:100%" border="0"

$$  \displaystyle \frac {1}{\rho c^2} \left( \frac {\partial p} {\partial S} \right)_\rho = - \frac {1}{\rho} \frac {1}{ \left( \partial p / \partial \rho \right)_S} \left( \frac {\partial p }{ \partial \rho}\right)_S \left( \frac {\partial \rho }{ \partial S}\right)_p = \frac{\beta T}{C_p} $$  (36) where (p.26, Batchelor, 1967 ) <span id="(37)">
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \left( \frac {\partial p} {\partial S} \right)_\rho = -\left( \frac {\partial p} {\partial \rho} \right)_S \left( \frac {\partial \rho} {\partial S} \right)_p $$  (37) In Eq.(36), the thermal expansion coefficient is defined as $$\displaystyle \beta = -1/\rho \left(\partial \rho /  \partial T \right)_p $$, which is a positive quatity as the density decreases with temperature. The constant pressure specific heat is defined as $$\displaystyle C_p = T \left(\partial S / \partial T \right)_p $$. The time rate of change of the entropy is given as (Bird et al., 2006, p.372 Eq.(11D.1-3) )
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

<span id="(38)">
 * {| style="width:100%" border="0"

$$  \displaystyle \rho \frac {DS}{Dt} = -\frac {1}{T} \left( \nabla \cdot \mathbf q \right) -\frac {1}{T} \left( \boldsymbol \tau : \nabla \mathbf v \right) $$  (38)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

= References =