2.13 Fluid Mechanics

Mass and Momentum Balance

In a spatial (Eulerian) frame, the momentum balance equation for a continuum is

\[ \begin{equation} \rho\mathbf{a}=\divg\boldsymbol{\sigma}+\rho\mathbf{b}\,,\label{eq:momentum-balance} \end{equation} \]

where \(\rho\) is the density, \(\boldsymbol{\sigma}\) is the Cauchy stress, \(\mathbf{b}\) is the body force per mass, and \(\mathbf{a}\) is the acceleration, given by the material time derivative of the velocity \(\mathbf{v}\) in the spatial frame,

\[ \begin{equation} \mathbf{a}=\dot{\mathbf{v}}=\frac{\partial\mathbf{v}}{\partial t}+\mathbf{L}\cdot\mathbf{v}\,,\label{eq:acceleration} \end{equation} \]

where \(\mathbf{L}=\grad\mathbf{v}\) is the spatial velocity gradient. The mass balance equation is

\[ \begin{equation} \dot{\rho}+\rho\divg\mathbf{v}=0\,,\label{eq:mass-balance} \end{equation} \]

where the material time derivative of the density in the spatial frame is

\[ \begin{equation} \dot{\rho}=\frac{\partial\rho}{\partial t}+\grad\rho\cdot\mathbf{v}\,.\label{eq:density-material-derivative} \end{equation} \]

Let \(\mathbf{F}\) denote the deformation gradient (the gradient of the motion with respect to the material coordinate). The material time derivative of \(\mathbf{F}\) is related to \(\mathbf{L}\) via

\[ \begin{equation} \dot{\mathbf{F}}=\mathbf{L}\cdot\mathbf{F}\,.\label{eq:F-dot} \end{equation} \]

Let \(J=\det\mathbf{F}\) denote the Jacobian of the motion (the volume ratio, or ratio of current to referential volume, \(J>0\)); then, the dilatation (relative change in volume between current and reference configurations) is given by \(e=J-1\). Using the chain rule, \(J\)'s material time derivative is \(\dot{J}=J\mathbf{F}^{-T}:\dot{\mathbf{F}}\) which, when combined with eq.\eqref{eq:F-dot}, produces a kinematic constraint between \(\dot{J}\) and \(\divg\mathbf{v}\),

\[ \begin{equation} \dot{J}=J\,\divg\mathbf{v}\,.\label{eq:divv-kinematic-relation} \end{equation} \]

Substituting this relation into the mass balance, eq.\eqref{eq:mass-balance}, produces \(\dot{\overline{\rho J}}=0\), which may be integrated directly to yield

\[ \begin{equation} \rho=\rho_{r}/J\,,\label{eq:mass-balance-integrated} \end{equation} \]

where \(\rho_{r}\) is the density in the reference configuration (when \(J=1\)). Since \(\rho_{r}\) is obtained by integrating the above material time derivative of \(\rho J\), it is an intrinsic material property that must be invariant in time and space.

The Cauchy stress is given by

\[ \begin{equation} \boldsymbol{\sigma}=-p\mathbf{I}+\boldsymbol{\tau}\,,\label{eq:stress} \end{equation} \]

where \(\mathbf{I}\) is the identity tensor, \(\boldsymbol{\tau}\) is the viscous stress, \(p\) is the pressure arising from the elastic response,

\[ \begin{equation} p=-\frac{d\Psi_{r}\left(J\right)}{dJ}\,,\label{eq:elastic-pressure} \end{equation} \]

and \(\Psi_{r}\) is the free energy density of the fluid (free energy per volume of the continuum in the reference configuration). The axiom of entropy inequality dictates that \(\Psi_{r}\) cannot be a function of the rate of deformation \(\mathbf{D}=\left(\mathbf{L}+\mathbf{L}^{T}\right)/2\). In contrast, the viscous stress \(\boldsymbol{\tau}\) is generally a function of \(J\) and \(\mathbf{D}\).

Boundary conditions may be derived by satisfying mass and momentum balance across a moving interface \(\Gamma\). Let \(\Gamma\) divide the material domain \(V\) into subdomains \(V_{+}\) and \(V_{-}\) and let the outward normal to \(V_{+}\) on \(\Gamma\) be denoted by \(\mathbf{n}\). The jump condition across \(\Gamma\) derived from the axiom of mass balance requires that

\[ \begin{equation} \left[\left[\rho\mathbf{u}_{\Gamma}\right]\right]\cdot\mathbf{n}=0\,,\label{eq:mass-jump} \end{equation} \]

where \(\mathbf{u}_{\Gamma}\equiv\mathbf{v}-\mathbf{v}_{\Gamma}\) on \(\Gamma\) and \(\mathbf{v}_{\Gamma}\) is the velocity of the interface \(\Gamma\). Thus, \(\mathbf{u}_{\Gamma}\) represents the velocity of the fluid relative to \(\Gamma\). The double bracket notation denotes \(\left[\left[f\right]\right]=f_{+}-f_{-}\), where \(f_{+}\) and \(f_{-}\) represent the value of \(f\) on \(\Gamma\) in \(V_{+}\) and \(V_{-}\), respectively. This jump condition implies that the mass flux normal to \(\Gamma\) must be continuous. In particular, if \(V_{+}\) is a fluid domain and \(V_{-}\) is a solid domain, and \(\Gamma\) denotes the solid boundary (e.g., a wall), we use \(\rho_{+}=\rho\), \(\mathbf{v}_{+}=\mathbf{v}\) for the fluid, and \(\mathbf{v}_{-}=\mathbf{v}_{\Gamma}\) for the solid, such that eq.\eqref{eq:mass-jump} reduces to \(\rho\left(\mathbf{v}-\mathbf{v}_{\Gamma}\right)\cdot\mathbf{n}=0\). The jump condition derived from the axiom of linear momentum balance similarly requires that

\[ \begin{equation} \left[\left[\boldsymbol{\sigma}-\rho\mathbf{u}_{\Gamma}\otimes\mathbf{u}_{\Gamma}\right]\right]\cdot\mathbf{n}=\mathbf{0}\,.\label{eq:momentum-jump} \end{equation} \]

This condition implies that the jump in the traction \(\boldsymbol{\sigma}\cdot\mathbf{n}\) across \(\Gamma\) must be balanced by the jump in momentum flux normal to \(\Gamma\). In addition to jump conditions dictated by axioms of conservation, viscous fluids require the satisfaction of the no-slip condition,

\[ \begin{equation} \left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)\cdot\left[\left[\mathbf{u}_{\Gamma}\right]\right]=\mathbf{0}\,,\label{eq:no-slip-condition} \end{equation} \]

which implies that the velocity component tangential to \(\Gamma\) is continuous across that interface.

In our finite element treatment we use \(\mathbf{v}\) and \(J\) as nodal variables, implying that our formulation automatically enforces continuity of these variables across element boundaries, thus \(\left[\left[\mathbf{v}\right]\right]=\left[\left[\mathbf{u}_{\Gamma}\right]\right]=\mathbf{0}\) and \(\left[\left[J\right]\right]=0\). Based on Eqs.\eqref{eq:mass-balance-integrated} and \eqref{eq:elastic-pressure}, it follows that the density and elastic pressure are continuous across element boundaries in this formulation, \(\left[\left[\rho\right]\right]=0\) and \(\left[\left[p\right]\right]=0\). Thus, the mass jump in eq.\eqref{eq:mass-jump} is automatically satisfied, and the momentum jump in eq.\eqref{eq:momentum-jump} reduces to \(\left[\left[\boldsymbol{\sigma}\right]\right]\cdot\mathbf{n}=\mathbf{0}\), requiring continuity of the traction, or more specifically according to \eqref{eq:stress}, the continuity of the viscous traction \(\boldsymbol{\tau}\cdot\mathbf{n}\), since \(p\) is automatically continuous.

Wall Shear Stress

Consider a no-slip impermeable wall (a wall on which the fluid velocity is equal to the wall velocity, \(\mathbf{v}=\mathbf{v}_{\Gamma}\)). Let \(\left\{ \mathbf{n},\mathbf{s},\mathbf{t}\right\}\) be an orthonormal basis where \(\mathbf{n}\) is the outward normal to the fluid on the wall, \(\mathbf{s}\) is a unit tangent to the wall along the local direction of the flow (in the immediate wall vicinity), and \(\mathbf{t}\) is a unit tangent to the wall in the direction orthogonal to \(\mathbf{n}\) and \(\mathbf{s}\). The components of the velocity gradient \(\mathbf{L}=\grad\mathbf{v}\) on the wall, when expressed in matrix form in the basis \(\left\{ \mathbf{n},\mathbf{s},\mathbf{t}\right\}\), are zero in the directions of the wall tangent (since there is no variation in any of the velocity components along those directions). Moreover, by construction, the tangential velocity component \(v_{t}\) along \(\mathbf{t}\) is zero, not only on the wall but also in its vicinity (thus, \(\partial v_{t}/\partial x_{n}=0\)). Combining these factors produces

\[ \begin{equation} \left[\mathbf{L}\right]=\left[\begin{array}{ccc} \frac{\partial v_{n}}{\partial x_{n}} & \frac{\partial v_{n}}{\partial x_{s}} & \frac{\partial v_{n}}{\partial x_{t}}\\ \frac{\partial v_{s}}{\partial x_{n}} & \frac{\partial v_{s}}{\partial x_{s}} & \frac{\partial v_{s}}{\partial x_{t}}\\ \frac{\partial v_{t}}{\partial x_{n}} & \frac{\partial v_{t}}{\partial x_{s}} & \frac{\partial v_{t}}{\partial x_{t}} \end{array}\right]=\left[\begin{array}{ccc} \frac{\partial v_{n}}{\partial x_{n}} & 0 & 0\\ \frac{\partial v_{s}}{\partial x_{n}} & 0 & 0\\ 0 & 0 & 0 \end{array}\right]\label{eq:wss-L} \end{equation} \]

where \(dx_{n}\), \(dx_{s}\) and \(dx_{t}\) represent line elements along each of the coordinate directions. The resulting rate of deformation tensor is

\[ \begin{equation} \left[\mathbf{D}\right]=\frac{1}{2}\left(\left[\mathbf{L}\right]+\left[\mathbf{L}\right]^{T}\right)=\left[\begin{array}{ccc} \frac{\partial v_{n}}{\partial x_{n}} & \frac{1}{2}\frac{\partial v_{s}}{\partial x_{n}} & 0\\ \frac{1}{2}\frac{\partial v_{s}}{\partial x_{n}} & 0 & 0\\ 0 & 0 & 0 \end{array}\right]\label{eq:wss=D} \end{equation} \]

Therefore the fluid stress evaluated from Eq.\eqref{eq:stress} will typically have the form

\[ \begin{equation} \left[\boldsymbol{\sigma}\right]=-p\left[\mathbf{I}\right]+\left(\kappa-\frac{2}{3}\mu\right)\left(\tr\mathbf{D}\right)\left[\mathbf{I}\right]+2\mu\left[\mathbf{D}\right]=\left[\begin{array}{ccc} -p+\left(\kappa+\frac{4}{3}\eta\right)\frac{\partial v_{n}}{\partial x_{n}} & \mu\frac{\partial v_{s}}{\partial x_{n}} & 0\\ \mu\frac{\partial v_{s}}{\partial x_{n}} & -p & 0\\ 0 & 0 & -p \end{array}\right]\label{eq:wss-stress} \end{equation} \]

where \(\kappa\) is the bulk viscosity and \(\mu\) is the shear viscosity (which may be a function of the rate of deformation in the case of a non-Newtonian fluid), see Eq.(5.16-1) in Section Viscous Fluids. From this result we conclude that the wall shear stress is

\[ \begin{equation} \sigma_{ns}=\sigma_{sn}=\mu\frac{\partial v_{s}}{\partial x_{n}}\label{eq:wss-shear-stress} \end{equation} \]

In practice it is inconvenient to evaluate this wall shear stress in a local coordinate system attached to the wall and directed along the flow in the wall vicinity. As it turns out, the principal normal stresses evaluated from \(\left[\boldsymbol{\sigma}\right]\) in Eq. \eqref{eq:wss-stress} are

\[ \begin{equation} \begin{aligned}\sigma_{1} & =-p & \sigma_{2} & =-p+\frac{1}{2}\left(\sigma_{nn}-\sqrt{\sigma_{nn}^{2}+4\sigma_{sn}^{2}}\right) & \sigma_{3} & =-p+\frac{1}{2}\left(\sigma_{nn}+\sqrt{\sigma_{nn}^{2}+4\sigma_{sn}^{2}}\right)\end{aligned} \label{eq:wss-principal-normal} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\sigma_{nn} & =-p+\left(\kappa+\frac{4}{3}\mu\right)\frac{\partial v_{n}}{\partial x_{n}} & \sigma_{sn} & =\mu\frac{\partial v_{s}}{\partial x_{n}}\end{aligned} \label{eq:wss-normal-shear-tractions} \end{equation} \]

Thus, the maximum shear stress is

\[ \begin{equation} \begin{aligned}\sigma_{s,\text{max}} & =\max\left(\frac{\left|\sigma_{1}-\sigma_{2}\right|}{2},\frac{\left|\sigma_{2}-\sigma_{3}\right|}{2},\frac{\left|\sigma_{3}-\sigma_{1}\right|}{2}\right)\\ & =\max\left(\frac{1}{4}\left|\sigma_{nn}-\sqrt{\sigma_{nn}^{2}+4\sigma_{sn}^{2}}\right|,\frac{1}{2}\left|\sqrt{\sigma_{nn}^{2}+4\sigma_{sn}^{2}}\right|,\frac{1}{4}\left|\sigma_{nn}+\sqrt{\sigma_{nn}^{2}+4\sigma_{sn}^{2}}\right|\right) \end{aligned} \label{eq:wss-max-shear-stress-general} \end{equation} \]

When the flow is incompressible it follows that \(\tr\mathbf{D}=\frac{\partial v_{n}}{\partial x_{n}}=0\), in which case \(\boldsymbol{\sigma}_{nn}=-p\) according to Eq.\eqref{eq:wss-normal-shear-tractions}. Thus,

\[ \begin{equation} \sigma_{s,\text{max}}=\max\left(\frac{\left|\sigma_{ns}\right|}{2},\left|\sigma_{ns}\right|,\frac{\left|\sigma_{ns}\right|}{2}\right)=\left|\sigma_{sn}\right|\label{eq:wss-max-shear-incompressible} \end{equation} \]

Therefore, when the flow is incompressible the wall shear stress is equal to the maximum fluid shear stress. This result makes it convenient to use the maximum fluid shear stress on the wall as the value of the wall shear stress. Conversely, if the flow is compressible (which is generally the case in FEBio) but as long as

\[ \begin{equation} \left|\sigma_{nn}\right|\ll2\left|\sigma_{sn}\right|\label{eq:wss-negligible-normal} \end{equation} \]

the conclusion remains the same. Therefore, in most applications one can use the maximum fluid shear stress as a method for estimating the wall shear stress. This result remains valid for problems where the wall is not stationary, such as fluid-structure interaction problems. In practice one can always examine the principal normal stresses in the fluid, along the wall, from which one can deduce if Eq.\eqref{eq:wss-negligible-normal} is satisfied satisfactorily.

However, if the wall is porous or if slippage is allowed, then the wall velocity components are not necessarily zero or negligible. Thus, the maximum shear stress is not always equal to the wall shear stress.

Energy Balance

The energy balance for a continuum may be written in integral form over a control volume \(V\) as

\[ \begin{equation} \begin{aligned}\frac{d}{dt}\int_{V}\rho\left(\varepsilon+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\,dV & =-\int_{S}\rho\left(\varepsilon+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\left(\mathbf{v}\cdot\mathbf{n}\right)\,dS+\int_{S}\mathbf{t}\cdot\mathbf{v}\,dS+\int_{V}\rho\mathbf{b}\cdot\mathbf{v}\,dV\\ & -\int_{S}\mathbf{q}\cdot\mathbf{n}\,dS+\int_{V}\rho r\,dV \end{aligned} \,,\label{eq:energy-balance-integral} \end{equation} \]

where \(S\) is the control surface bounding \(V\), \(\varepsilon\) is the specific internal energy, \(\mathbf{q}\) is the heat flux across \(S\), and \(r\) is the heat supply per mass to the material in \(V\) resulting from other sources. Bringing the time derivative inside the integral on the left-hand-side, and using the divergence theorem, this integral statement of the energy balance may be written as

\[ \begin{equation} \begin{aligned} & \int_{V}\left[\rho\left(\dot{\varepsilon}+\mathbf{v}\cdot\mathbf{a}\right)+\rho\left(\varepsilon+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\left(\divg\mathbf{v}-\frac{\dot{J}}{J}\right)\right]\,dV\\ & =\int_{V}\left[\boldsymbol{\sigma}:\mathbf{D}-\divg\mathbf{q}+\rho r+\mathbf{v}\cdot\left(\divg\boldsymbol{\sigma}+\rho\mathbf{b}\right)\right]\,dV \end{aligned} \,.\label{eq:energy-integral-redux} \end{equation} \]

This statement must be valid for arbitrary control volumes and arbitrary processes, from which we conventionally derive the differential form of the axioms of mass, momentum and energy balance.

For the specialized conditions of a viscous fluid at constant temperature assumed in our treatment, the only state variables for the functions of state \(\varepsilon\), \(\boldsymbol{\sigma}\) and \(\mathbf{q}\) are \(J\) and \(\mathbf{D}\) (i.e., the temperature is not a state variable since it is assumed constant). Under these conditions the entropy inequality shows that the specific entropy \(\eta\) and the heat flux \(\mathbf{q}\) must be zero, and the Cauchy stress \(\boldsymbol{\sigma}\) must have the form of eq.\eqref{eq:stress} where \(p\) is given by eq.\eqref{eq:elastic-pressure} as a function of \(J\) only, leaving the residual dissipation statement \(\boldsymbol{\tau}:\mathbf{D}\ge0\) as a constraint that must be satisfied by constitutive relations for \(\boldsymbol{\tau}\). (For a Newtonian fluid, this constraint is satisfied when the viscosities \(\mu\) and \(\kappa\) are positive.) From these thermodynamic restrictions we conclude that \(\varepsilon=\psi\), where \(\psi\) is the specific (Helmholtz) free energy, with \(\Psi_{r}=\rho_{r}\psi\).

For the conditions adopted here (isothermal viscous fluid), the axiom of energy balance reduces to \(\rho\dot{\psi}=\boldsymbol{\sigma}:\mathbf{D}+\rho r\); since \(\psi\) is only a function of \(J\), this expression may be further simplified using Eqs.\eqref{eq:divv-kinematic-relation}-\eqref{eq:elastic-pressure} to produce \(\boldsymbol{\tau}:\mathbf{D}+\rho r=0\). In other words, isothermal conditions may be maintained only if heat dissipated by the viscous stress is emitted in the form of a heat supply density \(\rho r=-\boldsymbol{\tau}:\mathbf{D}\) (heat leaving the system). Now, the integral form of the energy balance in eq.\eqref{eq:energy-integral-redux} simplifies to

\[ \begin{equation} \int_{V}\left[\mathbf{v}\cdot\left(\divg\boldsymbol{\sigma}+\rho\left(\mathbf{b}-\mathbf{a}\right)\right)+\rho\left(\psi+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\left(\frac{\dot{J}}{J}-\divg\mathbf{v}\right)\right]\,dV=0\,.\label{eq:energy-isothermal-viscous} \end{equation} \]

A comparison of this statement with the statement of virtual work, presented below in eq.(3.5-1), establishes a clear correspondence between the virtual velocity \(\delta\mathbf{v}\) and \(\mathbf{v}\), and between the virtual energy density \(\delta J\) and \(\rho\left(\psi+\frac{1}{2}\mathbf{v}\cdot\mathbf{v}\right)\), with the latter representing the sum of the internal (free) and kinetic energy densities.