Governing Equations and Numerical Foundations

The continuity equation expresses conservation of mass: mass is neither created nor destroyed.

Derivation

Take a fixed small control volume (a cube) of sides $dx,\ dy,\ dz$ with density $\rho$ and velocity components $u,v,w$.

        y
        ^   +---------+
        |  /         /|
        | +---------+ |   mass in  : rho*u dy dz
        | |         | |   at x face
        | |         | +   mass out : [rho*u +
        | |         |/      d(rho*u)/dx dx] dy dz
        | +---------+
        +-------------> x

Net mass efflux in the x-direction:

$$\left[\rho u+\frac{\partial(\rho u)}{\partial x}dx\right]dy\,dz-\rho u\,dy\,dz=\frac{\partial(\rho u)}{\partial x}dx\,dy\,dz$$

Similarly for y and z. Total net outflow of mass per unit time:

$$\left[\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}\right]dx\,dy\,dz$$

Rate of decrease of mass inside the element is $-\dfrac{\partial\rho}{\partial t}dx\,dy\,dz$. Equating the two and dividing by the volume:

$$\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}=0$$

or in vector form

$$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\vec V)=0$$

Special cases

• Steady flow ($\partial\rho/\partial t=0$):

$$\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}=0$$

• Incompressible flow ($\rho$ constant):

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=\nabla\cdot\vec V=0$$

This means the velocity field is divergence-free, and volume is conserved for a fluid element.

Given: $u=x^2y$, incompressible 2D flow, $v=0$ at $y=0$.

Continuity

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$ $$\frac{\partial u}{\partial x}=2xy\ \Rightarrow\ \frac{\partial v}{\partial y}=-2xy$$

Integrate with respect to $y$:

$$v=-xy^2+f(x)$$

The condition $v=0$ at $y=0$ for all $x$ gives $f(x)=0$. Hence

$$v=-xy^2$$

Vorticity

$$\zeta_z=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}=(-y^2)-(x^2)=-(x^2+y^2)$$

The vorticity is not zero, so the flow is rotational. Its magnitude increases with distance from the origin, and the minus sign means clockwise rotation.

Answer: $v=-xy^2$; $\zeta_z=-(x^2+y^2)$ s$^{-1}$ (rotational).

Integral form

For a fixed control volume (CV) with control surface (CS):

$$\frac{\partial}{\partial t}\int_{CV}\rho\,dV+\int_{CS}\rho\,(\vec V\cdot\hat n)\,dA=0$$

The first term is the rate of change of mass stored in the CV; the second is the net mass outflow through the surface. For steady flow the first term is zero, so total mass inflow equals total outflow. For one inlet and one outlet:

$$\rho_1A_1V_1=\rho_2A_2V_2=\dot m$$

The integral form is the basis of the finite volume method, since it is applied to each cell. The differential form follows by applying the divergence theorem to a very small volume.

Numerical

Mass flow rate:

$$\dot m=\rho_1A_1V_1=1.2\times0.020\times20=0.48\ \text{kg/s}$$

Outlet velocity:

$$V_2=\frac{\dot m}{\rho_2A_2}=\frac{0.48}{1.5\times0.010}=32\ \text{m/s}$$

Check: the velocity increases because the area is halved (factor 2) while the density rises by 1.25: $V_2=V_1\times\dfrac{\rho_1A_1}{\rho_2A_2}=20\times\dfrac{0.024}{0.015}=32$ m/s.

Answer: $\dot m=0.48$ kg/s; $V_2=32$ m/s.

The Navier-Stokes (N-S) equations express Newton's second law ($\text{mass}\times\text{acceleration}=\text{force}$) for a fluid element.

Equation (vector form)

$$\rho\left(\frac{\partial\vec V}{\partial t}+(\vec V\cdot\nabla)\vec V\right)=-\nabla p+\mu\nabla^2\vec V+\rho\vec g$$

In Cartesian form, x-direction:

$$\rho\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}\right)=-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}+\frac{\partial^2u}{\partial z^2}\right)+\rho g_x$$

with similar equations for $y$ and $z$. It is used together with $\nabla\cdot\vec V=0$.

Meaning of terms

The left side is mass times acceleration per unit volume; the right side is the sum of forces.

Assumptions

• Continuum, Newtonian fluid (stress proportional to strain rate).
• Constant density and viscosity (incompressible, isothermal).
• Stokes' hypothesis for bulk viscosity is not needed when $\nabla\cdot\vec V=0$.

Why difficult to solve

• The convective term $u\,\partial u/\partial x$ is non-linear, so superposition fails and no general analytical solution exists.
• Four unknowns ($u,v,w,p$) are coupled with four equations, and pressure has no equation of its own (needs pressure-velocity coupling).
• At high $Re$ the solution is turbulent, with a wide range of scales, so very fine meshes are needed (DNS) or modelling (RANS/LES).
• Analytical solutions exist only for simple cases such as Couette and Poiseuille flows, which are used to test CFD codes.

Derivation

Assumptions: steady, incompressible, fully developed ($\partial u/\partial x=0$), $v=w=0$, no body force effect along x, infinitely wide plates ($\partial/\partial z=0$).

   y=h  ======================  fixed plate
         ->    ->     ->
         ->   -->    --->   u(y)
         ->    ->     ->
   y=0  ======================  fixed plate

The x-momentum equation reduces to

$$0=-\frac{dp}{dx}+\mu\frac{d^2u}{dy^2}\ \Rightarrow\ \frac{d^2u}{dy^2}=\frac1\mu\frac{dp}{dx}$$

(the y-momentum equation shows that $p$ depends only on $x$). Integrate twice:

$$u=\frac{1}{2\mu}\frac{dp}{dx}y^2+C_1y+C_2$$

No-slip: $u=0$ at $y=0$ gives $C_2=0$; $u=0$ at $y=h$ gives $C_1=-\dfrac{1}{2\mu}\dfrac{dp}{dx}h$. So

$$u=-\frac{1}{2\mu}\frac{dp}{dx}\,y\,(h-y)$$

The profile is parabolic. Maximum velocity at $y=h/2$:

$$u_{max}=-\frac{h^2}{8\mu}\frac{dp}{dx}$$

Flow rate per unit width $q=\int_0^hu\,dy=-\dfrac{h^3}{12\mu}\dfrac{dp}{dx}$, so mean velocity $\bar u=q/h=\tfrac23u_{max}$. Wall shear stress $\tau_w=\mu\,\dfrac{du}{dy}\Big|_{0}=-\dfrac h2\dfrac{dp}{dx}$ in magnitude.

Numerical

Data: $\mu=0.08$ Pa·s, $\rho=870$ kg/m³, $h=0.006$ m, $G=-dp/dx=1500$ Pa/m.

$$u_{max}=\frac{G h^2}{8\mu}=\frac{1500\times(0.006)^2}{8\times0.08}=0.0844\ \text{m/s}$$ $$q=\frac{Gh^3}{12\mu}=\frac{1500\times(0.006)^3}{12\times0.08}=3.375\times10^{-4}\ \text{m}^2/\text{s per m width}$$ $$\bar u=\frac qh=\frac{3.375\times10^{-4}}{0.006}=0.05625\ \text{m/s}\quad(=\tfrac23\times0.0844)$$ $$\tau_w=\frac{Gh}{2}=\frac{1500\times0.006}{2}=4.5\ \text{Pa}$$ $$Re=\frac{\rho\bar uh}{\mu}=\frac{870\times0.05625\times0.006}{0.08}=3.67$$

$Re\ll1$, so the flow is laminar (creeping), and the assumption is valid.

Answer: $u_{max}=0.0844$ m/s; $q=3.375\times10^{-4}$ m²/s per m; $\bar u=0.0563$ m/s; $\tau_w=4.5$ Pa; $Re=3.67$.

The energy equation comes from the first law of thermodynamics applied to a fluid element: the rate of change of energy equals heat added plus work done on the element.

Enthalpy form

With specific enthalpy $h=e+p/\rho$, the transport equation is

$$\rho\frac{Dh}{Dt}=\frac{Dp}{Dt}+\nabla\cdot(k\nabla T)+\Phi+S_h$$

that is

$$\rho\left(\frac{\partial h}{\partial t}+\vec V\cdot\nabla h\right)=\frac{\partial p}{\partial t}+\vec V\cdot\nabla p+\nabla\cdot(k\nabla T)+\Phi+S_h$$

For an ideal gas or incompressible liquid $dh=c_p\,dT$, so the equation can be written in terms of temperature.

Simplification

For incompressible flow with constant $\rho,\ c_p,\ k$ and negligible viscous dissipation and pressure work (low speed):

$$\rho c_p\left(\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}\right)=k\nabla^2T+\dot q$$

or

$$\frac{\partial T}{\partial t}+\vec V\cdot\nabla T=\alpha\nabla^2T+\frac{\dot q}{\rho c_p},\qquad \alpha=\frac{k}{\rho c_p}$$

This is the convection-diffusion equation for temperature.

Special cases

• No flow ($\vec V=0$): $\partial T/\partial t=\alpha\nabla^2T$, the transient heat conduction equation.
• Steady, no flow, no source: $\nabla^2T=0$, the Laplace equation.
• Steady 1D conduction with source: $k\,d^2T/dx^2+\dot q=0$.

The ratio of convection to conduction is the Peclet number $Pe=\rho c_pVL/k=Re\,Pr$.

(a) Euler equations

For an inviscid fluid ($\mu=0$) the Navier-Stokes equation loses the viscous term:

$$\rho\left(\frac{\partial\vec V}{\partial t}+(\vec V\cdot\nabla)\vec V\right)=-\nabla p+\rho\vec g$$

In components (x-direction): $\dfrac{\partial u}{\partial t}+u\dfrac{\partial u}{\partial x}+v\dfrac{\partial u}{\partial y}+w\dfrac{\partial u}{\partial z}=-\dfrac1\rho\dfrac{\partial p}{\partial x}+g_x$.

Bernoulli from Euler: take steady flow along a streamline of length $ds$. The momentum equation along $s$ becomes

$$V\frac{dV}{ds}=-\frac1\rho\frac{dp}{ds}-g\frac{dz}{ds}$$

Multiply by $ds$ and integrate with constant $\rho$:

$$\int\frac{dp}{\rho}+\int V\,dV+\int g\,dz=\text{const}\ \Rightarrow\ \frac p\rho+\frac{V^2}{2}+gz=\text{constant}$$

Equivalently $p+\tfrac12\rho V^2+\rho gz=$ const along a streamline. It is valid for steady, incompressible, inviscid flow with no heat or work addition. If the flow is also irrotational, the constant is the same for the whole field.

In CFD, Euler solvers are used for high-speed external flows where viscous effects outside the boundary layer are small, as they are cheaper than full N-S solvers.

(b) Pitot tube

At the stagnation port the velocity is zero, so

$$p_0-p=\tfrac12\rho V^2$$ $$V=\sqrt{\frac{2(p_0-p)}{\rho}}=\sqrt{\frac{2\times450}{1.2}}=27.4\ \text{m/s}$$

(The Mach number is about 0.08, so incompressible treatment is valid.)

Answer: $V=27.4$ m/s.

Stokes flow

Stokes (creeping) flow is flow at very low Reynolds number ($Re\ll1$, in practice $Re<1$) in which inertia is negligible compared with viscous forces. The Navier-Stokes equation reduces to the linear form

$$\nabla p=\mu\nabla^2\vec V,\qquad \nabla\cdot\vec V=0$$

Conditions: small particle or length scale, high viscosity, or very low velocity (e.g. micro-fluidics, sedimentation, lubrication films). The solution is linear and time-reversible.

For a sphere of diameter $d$ moving at velocity $V$ in an infinite fluid, Stokes' drag is

$$F_D=3\pi\mu dV,\qquad C_D=\frac{24}{Re}$$

Terminal velocity

At terminal velocity: weight $-$ buoyancy $=$ drag

$$\frac{\pi d^3}{6}(\rho_p-\rho_f)g=3\pi\mu dV_t\ \Rightarrow\ V_t=\frac{(\rho_p-\rho_f)\,g\,d^2}{18\mu}$$

Substituting ($d=1\times10^{-4}$ m):

$$V_t=\frac{(2650-998)\times9.81\times(1\times10^{-4})^2}{18\times1.002\times10^{-3}}=\frac{1652\times9.81\times10^{-8}}{0.018036}=8.99\times10^{-3}\ \text{m/s}$$

Check of validity

$$Re=\frac{\rho_fV_td}{\mu}=\frac{998\times8.99\times10^{-3}\times10^{-4}}{1.002\times10^{-3}}=0.895<1$$

The particle Reynolds number is below 1, so Stokes' law is acceptable (slightly at the upper limit).

The drag force is $F_D=3\pi\mu dV_t=8.5\times10^{-9}$ N, equal to the net weight.

Answer: $V_t=8.99$ mm/s (about 9 mm/s), $Re=0.89$, so Stokes' law is valid.