quad8/continuummechanics/homework/hw8/main.tex

\documentclass[11pt,a4paper,dvipsnames]{article}
\usepackage[utf8]{inputenc}

\input{../preamble.tex}

\title{Homework 8\\Continuum Mechanics}
\author{Adrià Vilanova Martínez}
\date{May 23, 2021}

\showcorrectionstrue % Change "true" to "false" in order to show corrections as
                     % if they weren't corrections (in black instead of red).

\begin{document}

\maketitle

\begin{Problem}
  Two incompressible and immiscible fluids of similar densities $\rho_1$ and $\rho_2$, and of viscosities $\eta_1$ and $\eta_2$ are confined between two horizontal parallel plates. The gap between the plates is $H$. Since $\rho_1 \approx \rho_2$, we can safely neglect gravitational effects.

  Fluid 1, in contact with the lower plate, forms a layer of thickness $h_1$. Fluid 2, in contact with the upper plate, forms a layer on top of fluid 1 of thickness $H - h_1$. The upper plate moves with velocity $\vec{U} = U \hat{i}$.

  \begin{enumerate}[a)]
    \item Determine the steady-state velocity profile in each fluid layer. Sketch both velocity profiles assuming $\eta_2 = 2 \eta_1$.

    \item Compute the tangential force per unit area exerted on the upper plate.

    \item Think about the pressure and whether it is constant or not and why.
  \end{enumerate}
\end{Problem}

\correction{\textsc{Note}: I understood $\eta_2 = 2 \eta_1$ was a restriction throughout the whole problem, and not only for the velocity profiles.}

\textbf{Solution for a):} \\
Because of the temporal and spatial symmetries of the problem and some of the boundary conditions, the velocity field will only depend on $z$, so $\vec{v} \equiv \vec{v}(z)$, and we'll also have that $\vec{v} \cdot \hat{j} = \vec{v} \cdot \hat{k} = 0$.

Imposing Navier-Cauchy's equation in the case of incompressible fluids (for a generic viscosity $\eta$ and a velocity of the previous form) taken into account that we're in a steady state, we obtain:
\[ \begin{cases}
  0 = - \partial_x p + \eta v_x''(z), \\
  0 = - \partial_y p \implies p \not\equiv p(y), \\
  0 = - \partial_z p \implies p \not\equiv p(z).
\end{cases} \]

In the case of parallel flow in the $x$-direction in stationary flows, we know that $\partial_x P$ is constant throughout the volume of the fluid (and in fact equal to 0 due to the fact that the flow is only driven by the plates and is not pressure-driven). Thus, by integrating the first equation we have
\[ v_x(z) = c z + d, \]
which is the expression for the only non-zero component of the velocity field in both fluids (each fluid will have its own constants, which will be denoted by $c_i$ and $d_i$ where $i$ is the index of the fluid which has density $\rho_i$ and viscosity $\eta_i$).

Let's now impose all the boundary conditions: due to the fact that the fluids are real, there is a no slip condition in the interfaces with the plates (the tangential components of the velocity of the plate and fluid will be equal). This means:
\[ \begin{cases}
  v_x(H) = U, \\
  v_x(0) = 0
\end{cases} \implies \begin{cases}
  c_2 H + d_2 = U, \\
  d_1 = 0.
\end{cases} \implies \]
\[ \implies v_x(z) = \left\{\begin{array}{ll}
  c_1 z, & (z < h_1) \\
  c_2 (z - H) + U. & (z > h_1)
\end{array}\right. \]

At the interface between both fluids, since it is planar, both fluids are incompressible and Newtonian, and the velocity field is of the form $\vec{v} = v_x(z) \hat{i}$, we have seen in class that in that interface we have
\[ \eta_1 \partial_z v_x^{(1)} = \eta_2 \partial_z v_x^{(2)} \implies \partial_z v_x^{(1)} = 2 \partial_z v_x^{(2)}, \]
where the superscript indicates whether the expression of the velocity taken is the limit coming from the upper (2) or lower (1) section.

Thus, derivating and using the previous equality we obtain:
\[ c := c_1 = 2 c_2 \implies v_x(z) = \left\{\begin{array}{ll}
  \correction{2} c z, & (z < h_1) \\
  c (z - H) + U. & (z > h_1)
\end{array}\right. \]

\correction{I missed the following boundary condition which lets us determine $c$: we must impose continuity of the velocity at the liquid-liquid interface:
\[ v_x^\text{(1)} (h_1) = v_x^\text{(2)} (h_1) \implies 2 c h_1 = c(h_1 - H) + U \implies c = \frac{1}{h_1} (U - cH). \]}

At the interface between both fluids we can also impose the continuity of the tangential stresses. For that, let's calculate the stress tensor:
\[ \grad \vec{v} = v_x'(z) \hat{k} \hat{i} \implies \TT{e} = \frac{1}{2} v_x'(z) (\hat{k} \hat{i} + \hat{i} \hat{k}). \]
Finally, since $\TT{\sigma} = -p \TT{I} + \TT{\sigma}' = -p \TT{I} + 2 \eta \TT{e}$:
\[ \TT{\sigma} = -p \TT{I} + \eta v_x'(z) (\hat{k} \hat{i} + \hat{i} \hat{k}). \]
Now we can impose the boundary condition $(\TT{\sigma}_1 \cdot \hat{n}) \cdot \hat{n} = (\TT{\sigma}_2 \cdot \hat{n}) \cdot \hat{n} = 0$, which gives us:
\[ p_1 = p_2. \]
This means the pressure is identical in both fluids.

\correction{However, here we've supposed that superficial tension is 0. If that wasn't the case, although pressure would remain constants in each of the fluids, we would have a difference of pressures $\Delta p = 2 \gamma H$ at the interface. As the mean curvature is $H = 0$, we would also get $p_1 = p_2$ in this case, though.}

\begin{figure}[H]
  \centering
  \input{vprofile1.tex}
  \caption{Sketch of the velocity profile when $\eta_2 = 2 \eta_1$.}
\end{figure}

\textbf{Solution for b):} \\
The tangential force per unit area exerted on the upper plate is:
\[ \sigma_{xz}(H) = \eta v_x'(H) = \eta U. \]

\textbf{Solution for c):} \\
As shown in section a), in both fluids the pressure in all points throughout each fluid is constant, and not only that, but both pressures are equal to each other.

\newpage

\begin{Problem}
  A liquid of viscosity $\eta$ flows under the action of gravity through a cylindrical pipe of radius $R$ and length $L$. The pipe is inclined an angle $\alpha$ relative to the horizontal. There is also a pressure difference $\Delta p$ between the two ends of the pipe.

  \begin{enumerate}[a)]
    \item Compute the steady-state velocity profile.
    \item Find the flow rate across the pipe.
    \item Calculate the viscous drag force exerted by the fluid on the pipe.
  \end{enumerate}
\end{Problem}

Let's use cylindrical coordinates, $(r, \theta, z)$, with the symmetry axis of the pipe along $z$.

In general, $\vec{v} = (v_r, v_\theta, v_z)$, but considering laminar flow and symmetry arguments (no end effects), we have $\vec{v} = v_z(r, \theta, z) \hat{e}_z$ and $v_r = v_\theta = 0$.

In addition, the fluid is assumed to be incompressible:
\[ \div \vec{v} = 0 \implies \partial_z v_z = 0 \implies v_z \not\equiv v_z(z). \]

Due to symmetry, the velocity field has no angular dependence either, so that $v_z \not\equiv v_z(\theta)$. Hence:
\[ \vec{v} = v_z(r) \hat{e}_z, \]
where $v_z(r)$ is the radial velocity profile.

The convective acceleration is then identically zero (parallel shear flow): $(\vec{v} \cdot \nabla) \vec{v} = v_z \partial_z \vec{v} = 0$. The flow is also stationary: $\partial_t \vec{v} = 0$. The Navier-Stokes equation reduces to:
\[ 0 = - \grad p + \rho \vec{g} + \eta \nabla^2 \vec{v}. \]
For the three velocity components:
\[ \begin{cases}
  0 = - \partial_r p + \rho g \cos \alpha \cos \theta, \\
  0 = - \frac{1}{r} \partial_\theta p - \rho g \cos \alpha \sin \theta, \\
  0 = - \partial_z p + \rho g \sin \alpha + \eta \nabla^2 v_z.
\end{cases} \]

The first 2 equations simply determine the hydrostatic pressure profile, $p(r, \theta, z)$.

The presure gradient in the $z$ direction is constant, by translational invariance along $z$, and given by $\partial_z p = \frac{\Delta p}{L}$ with $\Delta p := p_L - p_0$. The last equation then reads:
\[ \eta \nabla^2 v_z = -\rho g \sin \alpha + \frac{\Delta p}{L} \implies \eta \frac{1}{r} \partial_r (r \partial_r) v_z(r) = - \rho g \sin \alpha + \frac{\Delta p}{L} \implies \]
\[ \implies \partial_r (r \partial_r) v_z(r) = \frac{1}{\eta} \left( - \rho g \sin \alpha + \frac{\Delta p}{L} \right) r \implies \]
\[ \implies r \partial_{rr} v_z(r) = \frac{1}{\eta} \left( - \rho g \sin \alpha + \frac{\Delta p}{L} \right) r - \partial_{r} v_z(r) \implies \]
\[ \implies \partial_{rr} v_z(r) = \frac{1}{\eta} \left( - \rho g \sin \alpha + \frac{\Delta p}{L} \right) - \frac{1}{r} \partial_{r} v_z(r) \implies \]
\[ \implies \partial_r v_z(r) = \frac{1}{2 \eta} \left( - \rho g \sin \alpha + \frac{\Delta p}{L} \right) r + \frac{A}{r} \implies \]
\[ \implies v_z(r) = \frac{1}{4 \eta} \left( - \rho g \sin \alpha + \frac{\Delta p}{L} \right) r^2 + A \log(r) + B. \]
The constants $A$, $B$ are determined by the boundary conditions:
\[ \begin{cases}
  v_z(0) \not\to \infty \implies A = 0, \\
  v_z(R) = 0 \text( (no-slip)) \implies B = -\frac{1}{4 \eta} \left( - \rho g \sin \alpha + \frac{\Delta p}{L} \right) R^2.
\end{cases} \]

Finally:
\[ v_z(r) = \frac{1}{4 \eta} \left( \frac{\Delta p}{L} - \rho g \sin \alpha \right) (r^2 - R^2). \]
This is a parabolic velocity profile, which results from the superposition of a profile due to $\frac{\Delta p}{L}$ and to $g$ separately, since the reduced Navier-Stokes equation is linear.

\textbf{Solution for b):} \\
Flow rate:
\[ Q = \int_0^{2 \pi} d\theta \int_0^R dr \, r v_z(r) = \frac{2 \pi}{4 \eta} \left( \frac{\Delta p}{L} - \rho g \sin \alpha \right) \int_0^R r (r^2 - R^2) \, dr \implies \]
\[ \implies Q = \frac{\pi}{8 \eta} R^4 \left( \rho g \sin \alpha - \frac{\Delta p}{L} \right). \]

\textbf{Solution for c):} \\
Drag force on the pipe:
\[ \frac{d \vec{F}}{dS} = \TT{\sigma}' \cdot \hat{n}, \quad \TT{\sigma}' = \eta (\grad \vec{v} + (\grad \vec{v})^T). \]

In this case the only non-zero elements of the viscous stress tensor are $\sigma'_{rz} = \sigma'_{zr} = \eta \partial_r v_z$. Hence:
\[ \frac{d \vec{F}}{dS} = \sigma'_{zr} |_{r = R} \hat{e}_z = \eta \left. \frac{\partial v_z}{\partial r} \right|_{r = R} \hat{e}_z = \frac{1}{2} \left( \frac{\Delta p}{L} - \rho g \sin \alpha \right) R \hat{e}_z. \]

Integrating over the surface of the pipe, we get the total drag force:
\[ F_z = \frac{1}{2} \left( \frac{\Delta p}{L} - \rho g \sin \alpha \right) R \int_0^L dz \int_0^{2 \pi} d\theta \, R, \]
so that
\[ F_z = \pi R^2 L \left( \frac{\Delta p}{L} - \rho g \sin \alpha \right); \quad \vec{F} = F_z \hat{e}_z. \]

\begin{Problem}
  Vorticity in Stokes' 1st problem - we have looked at this problem in class. Specifically, we obtained the velocity profile by solving the Navier-Stokes equation, which reduces in this problem to the diffusion equation.

  \begin{enumerate}[a)]
    \item Using the result in class, find the vorticity.

    \item Calculate the vorticity flux through rectangle in the $xy$-plane with dimension along $x$ equal to $L$. (Hint: use Stokes' theorem to find the corresponding circulation.)

    Note that the vorticity-flux is independent of time; it remains unchanged as the Stokes' boundary layer grows.

    \item Using what we called in class the alternative form of the Navier-Stokes equation, show that vorticity obeys a diffusion equation.

    This tells us that vorticity cannot be generated inside the viscous Stokes' layer during its growth, but only redistributed. We thus conclude that vorticity must arise at the plate surface during the (in reality, nearly) instantaneous accelereation of the fluid, and afterwards diffuse away from the plate into the fluid at large without changing the total vorticity-flux.
  \end{enumerate}
\end{Problem}

\textbf{Solution for a):} \\
In class we found that
\[ v_x(y, t) = v_0 \left[ 1 - \erf \left( \frac{y}{2 \sqrt{\nu t}} \right) \right]. \]
Therefore:
\[ \vec{\omega} = \curl \vec{v} = \begin{vmatrix}
  \hat{i} & \partial_x & v_0 \left[ 1 - \erf \left( \frac{y}{2 \sqrt{\nu t}} \right) \right] \\
  \hat{j} & \partial_y & 0 \\
  \hat{k} & \partial_z & 0
\end{vmatrix} = \frac{v_0}{\sqrt{\pi \nu t}} \exp\left[ - \left( \frac{y}{2 \sqrt{\nu t}} \right)^2 \right] \hat{k}. \]

\textbf{Solution for b):}
\[ I = \int_{R} (\curl \vec{v}) dS = \oint_{\partial R} \vec{v} \cdot dl = \left( \int_{AB} + \int_{BC} + \int_{CD} + \int_{DA} \right) \vec{v} \cdot dl, \]
where $A, B, C, D$ are the vertexs of the rectangle starting from the bottom-left corner (lowest $x$ and $y$) in the anti-clockwise direction.
Since $\vec{v} = v_x(y) \hat{i}$:
\[ I = \left( \int_{AB} + \bcancel{\int_{BC}} + \int_{CD} + \bcancel{\int_{DA}} \right) \vec{v} \cdot dl = L (v_x(y_1) - v_x(y_2)). \]

\correction{Considering the bottom side at the origin, and taking the limit of the upper side of the rectangle going to infinity, given that at the infinity limit $v_x = 0$ and at the origin $v_x = v_0$, we get
\[ I = L v_0. \]}

\textbf{Solution for c):} \\
The alternative form of the Navier-Stokes equation for incompressible flows is:
\[ \partial_t \vec{\omega} + \curl (\vec{\omega} \cross \vec{v}) = \frac{\eta}{\rho} \nabla^2 \vec{\omega} \]

We can use the identity for the curl of the vectorial product of 2 vector fields to see that the second term in the LHS vanishes, thus leaving us with the diffusion equation in 3D:
\[ \curl (\vec{\omega} \cross \vec{v}) = \vec{\omega} \cancel{\div \vec{v}} - \vec{v} \cancel{\div \vec{\omega}} + (\vec{v} \cdot \nabla) \vec{\omega} - (\vec{\omega} \cdot \nabla) \vec{v} = \]
\[ = v_x \cancel{\partial_x \vec{\omega}} - \omega_z \cancel{\partial_z \vec{v}} = 0, \]
where we have used the results in sections a) and b), the fact that the flow is incompressible and that $\div \vec{\omega} = \div (\curl \vec{v}) = 0$.

Therefore:
\[ \partial_t \vec{\omega} = \frac{\eta}{\rho} \nabla^2 \vec{\omega}. \]

\begin{Problem}
  Consider the kinetic energy density of an incompressible fluid, $\frac{1}{2} \rho v^2$. In this problem we are going to derive the conservation law for energy applied to a control (fixed) volume $V$.

  \begin{enumerate}[a)]
    \item Use the equation of motion for a fluid to show that:
    \[ \partial_t \left( \frac{1}{2} \rho v^2 \right) = \frac{v^2}{2} \partial_t \rho - \rho v_i v_j \partial_j v_i - v_i \partial_i p + v_i \partial_j \sigma'_{ij} + v_i f_i. \]
  \end{enumerate}
\end{Problem}

\textbf{Solution:}
\[ \partial_t \left( \frac{1}{2} \rho v^2 \right) = \frac{v^2}{2} \partial_t \rho + \rho v_i \partial_t v_i = \]
\[ = \frac{v^2}{2} \partial_t \rho + \rho v_i (- (\vec{v} \cdot \nabla) v_i + \frac{1}{\rho} f_i^* ) = \]
\[ = \frac{v^2}{2} \partial_t \rho + \rho v (- (v_j \partial_j) v_i + \frac{1}{\rho} (f_i + (\div \TT{\sigma}^T)_i) ) = \]
\[ = \frac{v^2}{2} \partial_t \rho - \rho v_i v_j \partial_j v_i - v_i \partial_i p + v_i \partial_j \sigma'_{ij} + v_i f_i. \]

\begin{FreeProblem}
  \begin{enumerate}[a)]
    \setcounter{enumi}{1}
    \item Now use the continuity equation to show that:
    \[ v_j \partial_j \left( \frac{\rho v^2}{2} \right) = v_j \rho v_i \partial_j v_i - \frac{v^2}{2} \frac{\partial \rho}{\partial t}. \]
  \end{enumerate}
\end{FreeProblem}

\Solution
First of all, we have:
\[ v_j \partial_j \left( \frac{\rho v^2}{2} \right) = v_j \rho v_i \partial_j v_i - \frac{v^2}{2} v_j \partial_j \rho. \]

Also, from the continuity equation it follows that
\[ \partial_t \rho = - \div (\rho \vec{v}) = - \rho \cancel{\div \vec{v}} - \vec{v} \cdot \grad \rho = - v_j \partial_j \rho, \]
which completes our proof.

\begin{FreeProblem}
  \begin{enumerate}[a)]
    \setcounter{enumi}{2}
    \item Use this result, that $\partial_j (v_i \sigma'_{ij}) = v_i \partial_j \sigma'_{ij} + \sigma'_{ij} \partial_j v_i$ and that $\TT{\sigma}'$ is symmetric, to obtain:
    \[ \partial_t \left( \frac{1}{2} \rho v^2 \right) = - (\vec{v} \cdot \grad) \left( \frac{1}{2} \rho v^2 + p \right) + \div (\TT{\sigma}' \cdot \vec{v}) - \TT{\sigma}' : \grad \vec{v} + \vec{v} \cdot \vec{f}. \]
  \end{enumerate}
\end{FreeProblem}

\textbf{Solution:}
\[ \partial_t \left( \rho \frac{v^2}{2} \right) = \cancel{\frac{v^2}{2} \partial_t \rho} - v_j \partial_j \left( \frac{\rho v^2}{2} \right) - \cancel{\frac{v^2}{2} \partial_t \rho} - v_i \partial_i p + \partial_j (v_i \sigma'_{ij}) - \sigma'_{ij} \partial_j v_i + v_i f_i = \]
\[ \notate[X]{{}={}}{1}{\scriptstyle v_i \sigma'_{ij} = v_i \sigma'_{ji} = \sigma'_{ji} v_i} - (\vec{v} \cdot \grad) \left( \frac{\rho v^2}{2} + p \right) + \div (\TT{\sigma}' \cdot \vec{v}) - \TT{\sigma}' : \grad \vec{v} + \vec{v} \cdot \vec{f}. \]

\begin{FreeProblem}
  \begin{enumerate}[a)]
    \setcounter{enumi}{3}
    \item Now use the identity $\div (\alpha \vec{A}) = \alpha \div \vec{A} + \vec{A} \cdot \grad \alpha$, with $\alpha$ a scalar and $\vec{A}$ a vector, to finally obtain:
    \[ \partial_t \left( \frac{1}{2} \rho v^2 \right) = - \div[ \vec{v} \left( \frac{\rho v^2}{2} + p \right) - \TT{\sigma}' \cdot \vec{v} ] - \TT{\sigma}' : \grad \vec{v} + \vec{v} \cdot \vec{f}. \]
  \end{enumerate}
\end{FreeProblem}

\Solution
As per the identity above:
\[ \div[ \vec{v} \left( \frac{\rho v^2}{2} + p \right) ] = \cancel{\left( \frac{\rho v^2}{2} + p \right) \underbrace{\div \vec{v}}_{= 0}} + (\vec{v} \cdot \grad) \left( \frac{\rho v^2}{2} + p \right) = (\vec{v} \cdot \grad) \left( \frac{\rho v^2}{2} + p \right). \]

\begin{FreeProblem}
  \begin{enumerate}[a)]
    \setcounter{enumi}{4}
    \item By integrating this equation to control (fixed) volume $V$, and using Gauss theorem, we arrive at the conservation law for the energy in volume $V$:
    \[ \begin{split}
      \frac{d}{dt} \int_V \frac{1}{2} \rho v^2 \, dV' = & - \oint_A \frac{1}{2} \rho v^2 \vec{v} \cdot \hat{n} \, dS - \oint_A p \vec{v} \cdot \hat{n} \, dS \\
      & \oint_A \vec{v} \cdot \TT{\sigma}' \cdot \hat{n} \, dS + \int_V \vec{v} \cdot \vec{f} \, dV' - \int_V \TT{\sigma}' : \grad \vec{v} \, dV'.
    \end{split} \]
  \end{enumerate}
\end{FreeProblem}

\Solution
Given that $V$ is fixed:
\[ \int_V \partial_t \left( \frac{1}{2} \rho v^2 \right) \, dV' = \frac{d}{dt} \int_V \left( \frac{1}{2} \rho v^2 \right) \, dV'. \]

Then, using Gauss' theorem:
\[ \int_V \div[ \vec{v} \left( \frac{\rho v^2}{2} + p \right) - \TT{\sigma}' \cdot \vec{v} ] \, dV' \lnotate[X]{{}={}}{1}{\scriptstyle \TT{\sigma}' \cdot \vec{v} = \vec{v} \cdot \TT{\sigma}' \text{ because } \TT{\sigma}' \text{ symm.}} \oint_A \left[ \frac{\rho v^2}{2} \vec{v} + p \vec{v} - \vec{v} \cdot \TT{\sigma}' \right] \cdot \hat{n}. \]

\begin{FreeProblem}
  \begin{enumerate}[a)]
    \setcounter{enumi}{5}
    \item Let's work with the last energy dissipation term. Show that for a Newtonian (and incompressible) fluid:
    \[ \TT{\sigma}' : \grad \vec{v} = \TT{\sigma}' : \TT{e} = 2 \eta \TT{e}^2. \]
  \end{enumerate}
\end{FreeProblem}

\textbf{Solution:} \\
\[ \TT{\sigma}' : \grad \vec{v} = \sigma'_{ij} \partial_k v_n \delta_{jk} \delta_{in} = \sigma'_{ij} \partial_j v_i = \frac{1}{2} (2 \sigma'_{ij} \partial_j v_i) = \]
\[ = \frac{1}{2} \sigma'_{ij} \partial_j v_i + \frac{1}{2} \sigma'_{ji} \partial_i v_j \notate[X]{{}={}}{1}{\scriptstyle \TT{\sigma}' \text{ symmetric}} \frac{1}{2} \sigma'_{ij} \partial_j v_i + \frac{1}{2} \sigma'_{ij} \partial_i v_j = \sigma'_{ij} \underbrace{\frac{1}{2} (\partial_i v_j + \partial_j v_i)}_{= e_{ji}} = \]
\[ = \sigma'_{ij} e_{ji} = \TT{\sigma}' : \TT{e}. \]

Then, for an incmopressible newtonian fluid we have $\TT{\sigma}' = 2 \eta \TT{e}$, and thus:
\[ \TT{\sigma}' : \TT{e} = 2 \eta \TT{e} : \TT{e} = 2 \eta e_{ij} e_{kn} \delta_{jk} \delta_{in} = 2 \eta e_{ij} e_{ji} \notate[X]{{}={}}{1}{\scriptstyle \TT{e} \text{ symmetric}} 2 \eta e_{ij} e_{ij} = 2 \eta \sum_{i, j} e_{ij}^2. \]

Therefore, we obtain:
\[ \TT{\sigma}' : \grad \vec{v} = \TT{\sigma}' : \TT{e} = 2 \eta \TT{e}^2. \]

\end{document}