$$\newcommand{\D}{\displaystyle} \renewcommand{\eqref}[1]{Eq.~(\ref{#1})}$$

# 3.5 Exercises

## Exercise 5: Stokes first problem for a non-Newtonian fluid

Figure 51: Stokes first problem for a non-Newtonian fluid: schematic representation of the velocity profile.

A non-Newtonian fluid flows in the plane with constant velocity $$U_0$$, paralell with the $$x$$-axis. At time $$t=0$$ we insert a thin plate into the flow, paralell to the $$x$$-axis. See Figure 51 (The plate must be concidered as infinitesimal thin, with infinite extent in th $$x$$-directions.) We are interested in finding the velocity field $$u(y,t)$$ which developes because of the adherence between the fluid and the plate. The equation of motion for the problem may be stated as: \begin{align} \rho \frac{\partial u}{\partial t} = \frac{\partial \tau_{xy}}{\partial y} \tag{3.152} \end{align}

The relation between the shear stress $$\tau_{xy}$$ and the velocitygradient $$\frac{\partial u}{\partial y}$$ is given by: \begin{align} \tau_{xy} = K \cdot \left| \frac{\partial u}{\partial y} \right|^{\frac{1}{2}} \frac{\partial u}{\partial y} \tag{3.153} \end{align}

where $$K$$ is a positive constant.

We introduce the following transformation which reduces the system to an ordinary differential equation: \begin{align} \eta = C \cdot \frac{y}{t^{\frac{2}{5}}}, \quad f(\eta) = \frac{u}{U_0} \tag{3.154} \end{align}

where $$C$$ is a positive constant, $$y \geq 0$$ and $$t \geq 0$$. By inserting equation (3.153) into equation (3.152) and using (3.154), we get: \begin{align} f''(\eta) + 2 \, \eta \sqrt{f'(\eta)} = 0 \tag{3.155} \end{align}

With boundaryconditions: \begin{align} f(0) = 0, f(\delta) = 1 \tag{3.156} \end{align}

where $$\delta$$ is the boundary-layer thickness.

Pen and paper: The following problems should be done using pen and paper:

a) We are first gonna solve equation (3.155) as an initial-value problem. In this problem we thus set $$f'(0)=1.25$$. Calculate $$f$$, and $$f'$$ for $$\eta=0.6$$ using Euler's method. Use $$\Delta \eta = 0.2$$

b) Write Eq (3.155) as a system of equations. Explain how you could solve the problem using shooting technique, when $$f(0)=0$$, $$f(\delta)=1$$, and $$\delta$$ is considered known.

c) During the solution of the system above we find that $$f'(\delta)=0$$. This information makes it possible to determine $$\delta$$ as a part of the solution of the system in b). Explain how.

d) Solve equation (3.155) when $$f_1(0)=f'(0)$$ is given. Show that the boundary condition $$f(\delta)=1$$ give $$f_1(0)= \frac{225}{128}=1.25311 \dots$$

Programming:

Write a python program that solves the problem defined in b. Use the information about $$\delta$$ obtained in c.

Pen and paper:

a) Solutions: $$f(0.6)=0.732$$, $$f'(0.6)=0.988$$.

c) Perform the differentiation $$\frac{d}{d \eta} \left[ \sqrt{f'(\eta)}\right]$$, and integrate equation (3.155) once. Solution: $$\delta=\sqrt{2 \sqrt{f'(0)}}$$

Programming:

You may use the template script below:

# src-ch2/stokes.py;ODEschemes.py @ git@lrhgit/tkt4140/src/src-ch2/ODEschemes.py;

#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()
import numpy as np

from ODEschemes import euler, heun, rk4
from math import sqrt
from scipy.interpolate import splev, splrep

#### set default plot values: ####
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT

#### a: ####
def func(f, eta):

f0 = f[0]
f1 = f[1]
if f1<0:
# f1 should not be negative
f1=0
df0 = ?
df1 = ?

dF = np.array([df0, df1])

return dF

N = 5

Eta = np.linspace(0, 1, N + 1)
d_eta = Eta[1] - Eta[0]

f_0 = ?
df_0 = ?

F_0 = [f_0, df_0]

F = euler(func, F_0, Eta)

print "f(0.6) = {0}, f'(0.6) = {1}".format(F[N*3/5, 0], F[N*3/5, 1])

#### b/c: ####

N = 100 # using more gridpoints

d_eta = Eta[1] - Eta[0]

f_0 = 0

s0, s1 = 1., 1.1
delta = ?

Eta = np.linspace(0, delta, N + 1)
F_0 = [?, ?]

F = euler(func, F_0, Eta)

phi0 = ?

# Solve using shooting technique:
for n in range(4):
F_0 = [?, ?]
delta = ?
Eta = np.linspace(0, delta, N + 1)

F = euler(func, F_0, Eta)

phi1 = ?

s = ?

s1, s0 = s, s1
print "n = {0}, s = {1}, ds = {2}, delta = {3} ".format(n, s0, s1 - s0, delta)

# Transform to time-space variables:

C = 1.
U0 = 1
dt = 0.05
t0, tend = dt, 2.

dy = 0.01
y0, yend = 0, 2.

time = np.arange(dt, tend + dt, dt)

Y = np.arange(y0, yend + dy, dy)

Usolutions = np.zeros((len(time), len(Y)))

f = F[:, 0]*U0
tck = splrep(Eta, f) # create interpolation splines

for n, t in enumerate(time):
Eta_n = C*Y/(t**(2./5))

Eta_n = np.where(Eta_n>delta, delta, Eta_n) # if eta>delta set eta to delta

Usolutions[n, :] = splev(Eta_n, tck)

from Visualization import myAnimation
myAnimation(Y, Usolutions, time)


You also need this file in order to run simulation:

# src-ch2/Visualization.py

import matplotlib.pylab as plt
import numpy as np
from matplotlib import animation

def myAnimation(Y, Usolutions, time):
fig = plt.figure()
ax = plt.axes(xlim=(0, 1.1), ylim=(0, Y[-1]))

lines=[]     # list for plot lines for solvers and analytical solutions
legends=[]   # list for legends for solvers and analytical solutions

line, = ax.plot([], [])
time_text = ax.text(0.1, 0.95, '', transform=ax.transAxes)
dt = time[1]-time[0]
plt.xlabel('u')
plt.ylabel('y')

# initialization function: plot the background of each frame
def init():
time_text.set_text('')
line.set_data([], [])
return lines,

# animation function.  This is called sequentially
def animate(i):
time = i*dt
time_text.set_text('time = %.4f' % time)
U = Usolutions[i, :]
line.set_data(U, Y)
return lines,

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, frames=len(time), init_func=init, blit=False)
#    Writer = animation.writers['ffmpeg']
#    writer = Writer(fps=15, metadata=dict(artist='Me'), bitrate=1800)
#    anim.save('../mov/stokes.mov', writer=writer)

plt.show()


## Exercise 6: The Falkner-Skan equation

Blasius equation is a special case of a more general equation called the Falkner-Skan equation. The Falkner-Skan Transformation (1930) allows to write the original equation in a similarity form given by: $$$$f'''+ff''+\beta\cdot[1-(f')^2]=0. \tag{3.157}$$$$

It can be shown that the velocity $$U(x)$$ is defined as: $$$$U(x)=U_0x^m,\ m= \text{constant}. \tag{3.158}$$$$

(3.157) is sometimes called the wedge flow equation because parameter $$\beta$$ has the geometric meaning as shown in Figure 52:

Figure 52: Schematic representation of the Falkner-Skan equation reference frame and problem setup.

We see that half the wedge angle is $$\beta\cdot \frac{\pi}{2}$$. For $$\beta=0$$, we get a flat plate and (3.157) becomes Blasius equation.

A detailed derivation can be found in Appendix C.3 in the Numeriske Beregninger. Here we reproduce some of the derivations for the sake of completeness.

Connection between $$m$$ and $$\beta$$: $$\beta=\frac{2m}{m+1}$$ and $$m=\frac{\beta}{2-\beta}$$. We see that $$\beta=0$$ results in $$U(x)=U_0$$.

Falkner-Skan transformation is given by: $$$$\eta=\sqrt{ \frac{U}{(2-\beta)\nu x} }\cdot y,\ f(\eta)=\frac{\psi}{\sqrt{(2-\beta)U\nu x}}, \tag{3.159}$$$$

where $$\psi$$ is a potential function such that $$\begin{equation*} u = \frac{\partial \psi}{\partial y}, \,\, v = \frac{\partial \psi}{\partial x} \end{equation*}$$

If we choose the same boundary conditions as in the example (3.3.3 Example: The Blasius equation), we get: $$$$f'''+ff''+\beta\cdot [1-(f')^2]=0, \tag{3.160}$$$$

with boundary conditions: $$$$\tag{3.161} \left.\begin{matrix} &f(0)=f'(0)=0, \\ &f'(\eta_\infty)=1. \end{matrix}\right.\$$$$