- September 28, Numerical Solutions to the Heat Equation
- September 30, Coding Heat Equation and Markov Chain Solvers
- October 3, Introduction to Fluid Dynamics
- October 5, More Navier-Stokes
- October 7, Still More Navier-Stokes!
- October 10, Nondimensionalization
- October 12, Reynolds Number from Navier-Stokes
- October 14, Review
- October 19, Final Project Selection
- October 24
- October 26
- October 28
- October 31
- November 2
- Homework and Projects

In [4]:

```
#PYTHON GENERIC IMPORTS
import numpy as np
import matplotlib.pyplot as plt
import random
import math
import ipywidgets
from ipywidgets import widgets
from IPython.display import display
```

- Dirichlet: $u(x,t)$ is given on the boundary
- Neumann: $\frac{du}{dt}$ is given on the boundary
- Robin: Linear combination of Dirichlet and Neumann (Ex: $\alpha u+\beta \frac{du}{dt}=c \,\text{on}\,\partial \Omega$)
- Periodic: $u(x,t)$ is periodic in the spacial domain

Note: $\Delta$ is the Laplacian and $\Delta = \nabla^2$ and $\nabla = \begin{pmatrix}\frac{\partial}{\partial x}\\ \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z}\end{pmatrix}$.

$$var = 2Dt$$

Normal Probability Distribution: $$f(x)=\frac{1}{2\sigma^2\pi}e^{\frac{-(x-\mu)^2}{2\sigma^2}}$$ Given a Heat Equation starting from $\delta_0$, $$u(x,t)=\frac{1}{\sqrt{4Dt\pi}}e^{\frac{-(x-0)^2}{4Dt}}$$

Moving to a numerical solution,
$$u(x,t)\approx \vec{u}(t) = \begin{pmatrix}u_1(t)\\u_2(t)\\\vdots\\u_n(t)\end{pmatrix}$$
$$\frac{d\vec{u}}{dt} = \frac{D}{h^2}\text{A}\vec{u}$$
With derivative matrix A:
$$A = \begin{pmatrix}
&&&\vec{b_0}\\
1 & -2 & 1 & 0 & 0 & \dots & 0\\
0 & 1 & -2 & 1 & 0 & \dots & 0\\
\vdots\\
0 & 0 & \dots & 0 & 1 & -2 & 1\\
&&&\vec{b_1}
\end{pmatrix}$$
The row vectors $\vec{b_0}$ and $\vec{b_1}$ control the boundary conditions. For Dirchlet, assume they are zero vectors.
For Periodic, $\vec{b_0} = \begin{pmatrix} -2 & 1 & 0 & \dots & 1\end{pmatrix}$ and $\vec{b_1} = \begin{pmatrix} 1 & 0 & \dots & 1 & -2\end{pmatrix}$

So, what we're doing is solving
$$\frac{\partial u}{\partial t} = \Delta u$$
$$\frac{u(x, t+h)-u(x,t)}{h} = Au(x,t)$$
$$u(x, t+h)= hAu(x,t) + u(x,t)$$
Which means we can calculate the spacial spread at time $t+h$ using the spacial second derivative matrix and the spacial spread at time $t$.

In [6]:

```
from enum import Enum
#size = 11
#boundary = np.concatenate([np.mat([-2,1]), np.zeros([1, size - 3])])
#print (boundary)
class BoundaryConditions(Enum):
periodic = 1
dirchlet = 0
def getLaplacianMatrix(boundary, size): #may not be stable -- be careful, not sure I understand laplacian transform matrix
A = np.zeros([size,size])
A += -2*np.identity(size)
A += np.diagonal(np.ones([size-1,1]), 1)
A += np.diagonal(np.ones([size,1]), -1)
if(boundary == BoundaryConditions.dirchlet):
A[0] = np.zeros([size, 1])
A[size-1] = np.zeros([size, 1])
elif(boundary == BoundaryConditions.periodic):
tmp = A[0]
A[0] = A[size-1]
A[size-1] = t
return A
def run1DHeatSim(initx, dt, steps, boundary):
r = [initx]
A = dt*getLaplacianMatrix(boundary, np.size(initx))
for i in range(0, steps):
r.append(A*r[i]+r[i])
return r
def run(size, dt, steps, boundary):
initx = np.zeros([size,1])
initx[round(size/2)] = 1
solution = run1DHeatSim(initx, dt, steps, boundary)
print(solution[steps-1])
#run(5, 1, 20, BoundaryConditions.dirchlet)
```

$p(\text{left}) = 0.45$

$p(\text{right}) = 0.55$

Walk for 80 steps

In [96]:

```
#RANDOM WALK in ONE DIMENSION
seeds = [874083802,29305632,94086608,-60502738,32466780,-8312131,72793714,35637054,53822097,-10185695,-94233408]
boundmin = -100
boundmax = 100
def walk1D(bias=0.55):
pos = 0
history = [pos]
t = [0]
steps = 100
boundmin = -10
boundmax = 10
random.seed(seeds[2])
for i in range(0,steps):
t.append(i)
pos = step(pos, bias, 1)
history.append(pos)
plt.plot(t, history, "#84B12F")
plt.xlabel("time")
plt.ylabel("position")
plt.show()
def roll(bias):
return (random.random() < bias)
def step(ipos, bias, size):
if(ipos == boundmin):
if(roll(bias)):
return boundmax - size
else:
return ipos - size
if(ipos == boundmax):
if(roll(bias)):
return ipos + size
else:
return boundmin + size
if(roll(bias)):
return ipos + size
return ipos - 1
def update(x):
walk1D(x)
widgets.interact(update, x=0.55)
```

Out[96]:

In [41]:

```
#FOR FUN, RANDOM WALK in TWO DIMENSIONS
def walk2D(xpos=0, ypos=0, steps=10000, seed=4, boundmin=-100, boundmax=100):
xhistory = [xpos]
yhistory = [ypos]
random.seed(seeds[seed])
for i in range(0,steps):
stp = random.random()*2
ypos = step(ypos, 0.5, stp)
xpos = step(xpos, 0.5, stp)
xhistory.append(xpos)
yhistory.append(ypos)
plt.plot(yhistory, xhistory, "#84B12F")
plt.xlabel("x position")
plt.ylabel("y position")
plt.show()
_steps = 1000000
_seed = 4
_boundmin = -100
_boundmax = 100
def func(btn):
walk2D(steps=_steps, seed=_seed, boundmin=_boundmin, boundmax=_boundmax)
def fsteps(x):
_steps = x
def fseed(x):
_seed = x
def fmin(x):
_boundmin = x
def fmax(x):
_boundmax = x
widgets.interact(fsteps, x=1000000)
widgets.interact(fseed, x=4)
widgets.interact(fmin, x=-100)
widgets.interact(fmax, x=100)
btn = widgets.Button(description="Run")
btn.on_click(func)
display(btn)
```

In physics, the **Navier–Stokes equations**, named after Claude-Louis Navier and George Gabriel Stokes, describe the motion of viscous fluid substances. These balance equations arise from applying Newton's second law to fluid motion, together with the assumption that the stress in the fluid is the sum of a diffusing viscous term (proportional to the gradient of velocity) and a pressure term—hence describing *viscous flow*. The main difference between them and the simpler Euler equations for inviscid flow is that Navier–Stokes equations also in the Froude limit (no external field) are not conservation equations, but rather a dissipative system, in the sense that they cannot be put into the quasilinear homogeneous form:

From our assumption of conservation of mass: $$\frac{\partial m}{\partial t} = \frac{d}{dt}\int_{\chi(\omega, t)}\rho(\vec{x},t) = 0$$

And from conservation of momentum: $$p = \int_{\chi(\omega,t)}\rho\vec{u}d\vec{x}$$ $$\frac{dp}{dt} = \frac{d}{dt}\int_{\chi(\omega,t)}\rho\vec{u}d\vec{x} = $$ Which by the transport theorem is $$ = \int_{\chi(\omega,t)} \frac{\partial}{\partial t}(\rho\vec{u}) + \nabla\cdot(\vec{u}\otimes\vec{u}\rho)d\vec{x}$$

$$\int_{\chi(\omega,t)}\rho\left(\frac{\partial \vec{u}}{\partial t} +\vec{u}\cdot\nabla\vec{u}\right)d\vec{x} = \int_{\chi(\omega,t)} g(\omega, \vec{x}) d\vec{x}$$Where $g(\omega, \vec{x})$ is the function of all forces acting on the fluid.

(note: $\text{stress}=\frac{\text{force}}{\text{area}}$)

Using the divergence theorem,
$$\int_{\chi}\rho\frac{D\vec{u}}{Dt}d\vec{x}=\int_{\chi}\nabla\cdot\sigma\vec{x}$$
The stress tensor can be built by
$$\sigma = -PI+\tau$$
With $-P$ as a *scalar* for pressure and $\tau$ as a matrix of the viscous forces.

So $$\nabla\cdot\sigma = \nabla\cdot(-PI) + \nabla\cdot\tau$$ $$\nabla\cdot(-PI) = \begin{pmatrix}\frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z}\end{pmatrix} \cdot \begin{pmatrix}-P&0&0\\0&-P&0\\0&0&-P\end{pmatrix} = \begin{pmatrix}-\frac{\partial P}{\partial x}\\-\frac{\partial P}{\partial y}\\-\frac{-\partial P}{\partial z}\end{pmatrix}$$

For Newtonian Fluids, we can write the viscous forces as $$\tau=\mu\left(\nabla\vec{u}+(\nabla\vec{u})^T\right)$$

Extending the viscous forces, $$\tau=\mu\left(\nabla\vec{u}+(\nabla\vec{u})^T\right) = $$ $$ = \mu\begin{pmatrix} \frac{\partial \vec{u_1}}{\partial x} & \frac{\partial \vec{u_1}}{\partial y} & \frac{\partial \vec{u_1}}{\partial z}\\ \frac{\partial \vec{u_2}}{\partial x} & \frac{\partial \vec{u_2}}{\partial y} & \frac{\partial \vec{u_2}}{\partial z}\\ \frac{\partial \vec{u_3}}{\partial x} & \frac{\partial \vec{u_3}}{\partial y} & \frac{\partial \vec{u_3}}{\partial z} \end{pmatrix} +\mu \begin{pmatrix} \frac{\partial \vec{u_1}}{\partial x} & \frac{\partial \vec{u_2}}{\partial x} & \frac{\partial \vec{u_3}}{\partial x}\\ \frac{\partial \vec{u_1}}{\partial y} & \frac{\partial \vec{u_2}}{\partial y} & \frac{\partial \vec{u_3}}{\partial y}\\ \frac{\partial \vec{u_1}}{\partial z} & \frac{\partial \vec{u_2}}{\partial z} & \frac{\partial \vec{u_3}}{\partial z} \end{pmatrix} =$$ $$ = \begin{pmatrix} \mu\left(\frac{\partial \vec{u_1}}{\partial x}+\frac{\partial \vec{u_1}}{\partial x} \right) & \mu\left(\frac{\partial \vec{u_1}}{\partial y}+\frac{\partial \vec{u_2}}{\partial x} \right) & \mu\left(\frac{\partial \vec{u_1}}{\partial z}+\frac{\partial \vec{u_3}}{\partial x} \right) \\ \mu\left(\frac{\partial \vec{u_2}}{\partial x}+\frac{\partial \vec{u_1}}{\partial y} \right) & \mu\left( \frac{\partial \vec{u_2}}{\partial y}+\frac{\partial \vec{u_2}}{\partial y} \right) & \mu\left(\frac{\partial \vec{u_2}}{\partial z}+\frac{\partial \vec{u_3}}{\partial y} \right) \\ \mu\left(\frac{\partial \vec{u_1}}{\partial z}+\frac{\partial \vec{u_3}}{\partial x} \right) & \mu\left(\frac{\partial \vec{u_3}}{\partial y}+\frac{\partial \vec{u_2}}{\partial z} \right) & \mu\left(\frac{\partial \vec{u_3}}{\partial z}+\frac{\partial \vec{u_3}}{\partial z} \right) \end{pmatrix}$$

Writing this matrix as

$$\mu \begin{pmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \end{pmatrix}$$We can interpret $\tau_{nm}$ as the stress acting on the $n$ face from the $m$ direction.

Finally, taking the divergence of $\tau$ $$(\nabla\cdot\tau)= $$

So, finally, $$\rho\left(\frac{\partial\vec{u}}{\partial t} + \vec{u}\cdot\nabla\vec{u}\right)=-\nabla P + \mu\Delta\vec{u}$$

Nondimensionalization can also recover characteristic properties of a system. For example, if a system has an intrinsic resonance frequency, length, or time constant, nondimensionalization can recover these values. The technique is especially useful for systems that can be described by differential equations. One important use is in the analysis of control systems. One of the simplest characteristic units is the doubling time of a system experiencing exponential growth, or conversely the half-life of a system experiencing exponential decay; a more natural pair of characteristic units is mean age/mean lifetime, which correspond to base e rather than base 2.

Many illustrative examples of nondimensionalization originate from simplifying differential equations. This is because a large body of physical problems can be formulated in terms of differential equations.

- List of dynamical systems and differential equations topics
- List of partial differential equation topics
- Differential equations of mathematical physics

Although nondimensionalization is well adapted for these problems, it is not restricted to them. An example of a non-differential-equation application is dimensional analysis; another example is normalization in statistics.

Measuring devices are practical examples of nondimensionalization occurring in everyday life. Measuring devices are calibrated relative to some known unit. Subsequent measurements are made relative to this standard. Then, the absolute value of the measurement is recovered by scaling with respect to the standard.

The Reynolds number is defined as the ratio of inertial forces to viscous forces and consequently quantifies the relative importance of these two types of forces for given flow conditions.[5] Reynolds numbers frequently arise when performing scaling of fluid dynamics problems, and as such can be used to determine dynamic similitude between two different cases of fluid flow. They are also used to characterize different flow regimes within a similar fluid, such as laminar or turbulent flow:

laminar flow occurs at low Reynolds numbers, where viscous forces are dominant, and is characterized by smooth, constant fluid motion; turbulent flow occurs at high Reynolds numbers and is dominated by inertial forces, which tend to produce chaotic eddies, vortices and other flow instabilities. In practice, matching the Reynolds number is not on its own sufficient to guarantee similitude. Fluid flow is generally chaotic, and very small changes to shape and surface roughness can result in very different flows. Nevertheless, Reynolds numbers are a very important guide and are widely used.

Reynolds number interpretation has been extended into the area of arbitrary complex systems as well: financial flows,[6] nonlinear networks,[7] etc. In the latter case an artificial viscosity is reduced to nonlinear mechanism of energy distribution in complex network media. Reynolds number then represents a basic control parameter which expresses a balance between injected and dissipated energy flows for open boundary system. It has been shown [7] that Reynolds critical regime separates two types of phase space motion: accelerator (attractor) and decelerator. High Reynolds number leads to a chaotic regime transition only in frame of strange attractor model.

The Reynolds number can be defined for several different situations where a fluid is in relative motion to a surface.[n 1] These definitions generally include the fluid properties of density and viscosity, plus a velocity and a characteristic length or characteristic dimension. This dimension is a matter of convention – for example radius and diameter are equally valid to describe spheres or circles, but one is chosen by convention. For aircraft or ships, the length or width can be used. For flow in a pipe or a sphere moving in a fluid the internal diameter is generally used today. Other shapes such as rectangular pipes or non-spherical objects have an equivalent diameter defined. For fluids of variable density such as compressible gases or fluids of variable viscosity such as non-Newtonian fluids, special rules apply. The velocity may also be a matter of convention in some circumstances, notably stirred vessels. The Reynolds number is defined as:

t is in seconds, set $T\tau=t$ such that $\tau$ is dimensionless. $$f=f(t(\tau))$$ $$\frac{df}{dt} = \frac{d\tau}{dt}\frac{df}{d\tau}$$ So, $$\frac{d}{dt}=\frac{1}{T}\frac{d}{d\tau}$$ And, $$\frac{d^2}{dt^2}=\left(\frac{1}{T}\frac{d}{d\tau}\right)^2=\frac{1}{T^2}\frac{d^2}{d\tau^2}$$ With these, we can find the system $$m\frac{1}{T^2}\frac{d^2}{d\tau^2}(L\chi) +b\frac{1}{T}\frac{d}{d\tau}(L\chi)+k(L\chi)=0$$ Dividing by $\frac{mL}{T^2}$, $$\frac{d^2\chi}{d\tau}+\frac{bT}{m}\frac{d\chi}{d\tau}+\frac{k}{m}T^2\chi=0$$ Now we can set $T=\sqrt{\frac{m}{k}}$ so that $$\frac{d^2\chi}{d\tau^2}+\frac{b}{\sqrt{mk}}\frac{d\chi}{d\tau}+\chi=0$$ Finally, just to make the equation look nice, we can define $\zeta=\frac{b}{2\sqrt{mk}}$ so, $$\frac{d^2\chi}{d\tau^2}+2\zeta\frac{d\chi}{d\tau}+\chi=0$$ This equation can be solved as $$\left(\alpha^2+2\zeta\alpha+1\right)e^{\alpha\tau}=0$$ And, $$\alpha=-\zeta\pm\sqrt{\zeta^2-1}$$ So, $$\chi(\tau)=c_1e^{\left(-\zeta+\sqrt{\zeta^2-1}\right)\tau}+c_2e^{\left(-\zeta-\sqrt{\zeta^2-1}\right)\tau}$$

Review.

- Gillespie Algorithm
- Stochastic Differential Equations
- Integral Spread Equations (related to Einstein's derivation of the heat equation)
- Particle Image Velocimetry
- Delayed Differential Equations
- Graph Theory and Dynamic Graphs

notes and discussion

Also, Paper Rough Draft and Group Notes

Problem 4, ManyWalk2D

In [86]:

```
#function many_walk2D(N, M);
#Simulates many random walkers that have an equally likely chance of moving
#left, right, up or down.
#N is the number of time steps
#M is the number of walkers
N = 20
M = 3
x = np.zeros([N,1]) #initialize the vector to hold the x-positions
y = np.zeros([N,1]) #initialize the vector to hold the y-positions
x[0] = 0 #set initial x-position of walker
y[0] = 0 #set initial y-position of walker
#set parameters here
stepsize = 1 #stepsize of walker
bins = (N/80) #size of bins
binmax = (N/2) #maximum bin
binnum = round(binmax/bins) #number of bins
random.seed(seeds[0])
rad = np.zeros([M,1]) #final distance of walker
nn = np.arange(0, bins, binmax)
size1 = np.size(nn)
for j in range(1, M): #loop through all walkers
for i in range(2, N): #loop through all of the time steps
r = random.random() #pick a random number to determine up/down or right/left
if r > 0.5: #go right/left
r1 = random.random() #pick a random number to determine right or left
if r1 > 0.5: #go right
y[i] = y[i-1]
x[i] = x[i-1] + stepsize
else: #go left
y[i] = y[i-1]
x[i] = x[i-1] - stepsize
else: #go up/down
r2 = random.random() #pick a random number to go up or down
if r2 > 0.5: #go up
x[i] = x[i-1]
y[i] = y[i-1] + stepsize;
else: #go down
x[i] = x[i-1]
y[i] = y[i-1] - stepsize;
rad[j] = math.sqrt(x[N-1]**2+y[N-1]**2) #This holds the final distance from the starting point of each walker
#use these to set the axes
xmax = max(abs(x))+1;
ymax = max(abs(y))+1;
#figure(1)
plt.plot(x, y, 'r'); #plot the path
plt.axis([-N/2, N/2, -N/2, N/2]);
plt.title('Two-dimensional ' + str(N) + '-step random walk')
plt.show()
#zoom on
#calculate theoretical distribution
x1 = np.zeros([100, 1])
pr = np.zeros([100, 1])
dx = nn[size1-1]/100
for k in range(0, 99):
x1[k]= (k-1)*dx;
pr[k]= ( 2*x1[k] / ((stepsize^2)*N))* math.exp(-1*x1[k]*x1[k]/((stepsize^2)*N) )
pr = pr * (bins*M); #rescale theoretical distribution
#figure(2)
#plot histogram and theoretical distribution
plt.hist(rad, nn)
plt.plot(x1, pr, 'r')
plt.xlabel('Distance')
plt.ylabel('Number of walkers')
plt.show()
```

Problem 5, (chain.m translated to python *by hand*)

In [95]:

```
# Thanks for making me do this conversion by hand...
# chain.m simulates a Markov chain on {1, 2, ..., n} given an initial
# distribution and a transition matrix.
# Program to simulate a Markov chain
# Below is a sample
k = 3 #number of states
v = np.mat([1/3, 1/3, 1/3]) # initial distribution
M = np.mat([[.5, 0.5, 0], [0.5, 0.5, 0], [0, 0.5, 0.5]]) # Markov transition matrix
n = 80 # number of time steps to take
x = v
t = [range(0, n)] # time indices
#print(M)
#print(v)
#Calculate v for n times steps.
for i in range (0,n):
x = np.concatenate((x,(x[i]*M)), axis=0)
#print (x)
xt = np.transpose(x)
#print(np.transpose(xt[0]))
#plot the results
plt.plot(t, xt[0], 'ro')
plt.plot(t, xt[1], 'go')
plt.plot(t, xt[2], 'bo')
plt.axis([0,n,0,1]);
plt.xlabel('time step')
plt.ylabel('probability')
#plt.legend('Island A', 'Island B', 'Island C')
plt.show()
```

Problem 6; Chain Final

In [110]:

```
# Program to simulate a Markov chain
# Below is a sample
k = 3 #number of states
v = np.mat([0.25, 0.25, 0.25, 0.25]) # initial distribution
M = np.mat([[0.00, 0.33, 0.34, 0.33],
[0.34, 0.00, 0.33, 0.33],
[0.33, 0.33, 0.00, 0.34],
[0.33, 0.34, 0.33, 0.00]]) # Markov transition matrix
n = 80 # number of time steps to take
x = v
t = [range(0, n)] # time indices
#print(M)
#print(v)
#Calculate v for n times steps.
for i in range (0,n):
x = np.concatenate((x,(x[i]*M)), axis=0)
#print (x)
xt = np.transpose(x)
#print(np.transpose(xt[0]))
#plot the results
plt.plot(t, xt[0], 'r.', label="Island A")
plt.plot(t, xt[1], 'g.', label="Island B")
plt.plot(t, xt[2], 'b.', label="Island C")
plt.plot(t, xt[3], 'y.', label="Island D")
plt.axis([0,n,0,1]);
plt.xlabel('time step')
plt.ylabel('proportion of population')
plt.legend(['Island A', 'Island B', 'Island C', 'Island D'])
plt.title("Question 6 part B; No one stays")
plt.show()
```

In [42]:

```
#oh joy another matlab file to transcribe. One of the two of us is really stubborn.
#(It might be the one writing pithy comments no living human but himself will ever read)
def run(numx, periodic):
numx = 101 #number of grid points in x
numt = 2000 #number of time steps to be iterated
dx = 1/(numx)
dt = 0.00005
pi = math.pi
x = np.arange(0,1,dx) #vector of x values, to be used for plotting
periodic = False
C = np.zeros([numx,numt]) #initialize everything to zero
#specify initial conditions
t = [0] #t=0
for i in range(0,numx-1):
C[i] = math.cos(2*pi*x[i]-pi)+1
#iterate difference equations
if(periodic):
for j in range(0, numt-1):
t.append(t[j] + dt) #update time
for i in range(0, numx-1): #loop through the middle spatial steps
C[i][j+1] = C[i][j] + (dt/dx**2)*(C[i+1][j] - 2*C[i][j] + C[i-1][j])
tmp = C[0][j+1]
C[0][j+1] = C[numx-1][j+1]
C[numx-1][j+1] = tmp
else:
for j in range(0, numt-1):
t.append(t[j] + dt) #update time
for i in range(1, numx-2): #loop through the middle spatial steps
C[i][j+1] = C[i][j] + (dt/dx**2)*(C[i+1][j] - 2*C[i][j] + C[i-1][j])
C[0][j+1] = 0 #Dirichlet boundary conditions
C[numx-1][j+1] = 0 #Dirichlet boundary conditions
D = np.transpose(C)
plt.plot(x,D[1], '-')
plt.plot(x,D[9], 'b-')
plt.plot(x,D[99], 'r--')
plt.plot(x,D[999], 'y-.')
plt.plot(x,D[1999], 'g:')
plt.xlabel('x')
plt.ylabel('u(x,t)')
plt.show()
numx = 0
periodic = False
```

In this problem you will consider two different fluid models and find the fluid flow profile within a rectangular channel, as in the figure below

a)

In the first part, we consider the inviscid, incompressible Navier-Stokes equations. As the name suggests, the inviscid models assumes the fluid does not undergo any shearing-typeforces. This model is also known as the Euler Equations. The governing equations for this situation are given by the following equations

Find the velocity of the fluid through the channel.

b)

Now in this part of the problem, we will consider the viscous (Newtonian), incompressible Navier-Stokes equations to model the situation in Figure 1. Again, as the name suggests here, we will reintroduce viscous forces back into the model. The governing equations for this situation are given by the following equations

$$u_x(\vec{x},t)
=-\left(\frac{\partial P(\vec{x},t)}{\partial x}\frac{1}{\rho} + \frac{\partial u_x(\vec{x},t)}{\partial t} + \frac{\partial u_x(\vec{x},t)}{\partial y}\cdot u_y(\vec{x},t) + \frac{\partial u_z(\vec{x},t)}{\partial z}\cdot u_z(\vec{x},t)\right) \left(\frac{\partial u_x(\vec{x},t)}{\partial x}\right)^{-1}$$

Now we add the viscous forces to the left hand side
$$
\rho
\left(
\begin{pmatrix} \frac{\partial u_x(\vec{x},t)}{\partial t} \\ \frac{\partial u_y(\vec{x},t)}{\partial t} \\ \frac{\partial u_z(\vec{x},t)}{\partial t} \end{pmatrix}
+
\begin{pmatrix}
\frac{\partial u_x(\vec{x},t)}{\partial x} & \frac{\partial u_x(\vec{x},t)}{\partial y} & \frac{\partial u_x(\vec{x},t)}{\partial z}\\
\frac{\partial u_y(\vec{x},t)}{\partial x} & \frac{\partial u_y(\vec{x},t)}{\partial y} & \frac{\partial u_y(\vec{x},t)}{\partial z}\\
\frac{\partial u_z(\vec{x},t)}{\partial x} & \frac{\partial u_z(\vec{x},t)}{\partial y} & \frac{\partial u_z(\vec{x},t)}{\partial z}
\end{pmatrix}
\begin{pmatrix}u_x(\vec{x},t)\\u_y(\vec{x},t)\\u_z(\vec{x},t)\end{pmatrix}
\right)=
-\begin{pmatrix} \frac{\partial P(\vec{x},t)}{\partial x} \\ \frac{\partial P(\vec{x},t)}{\partial y} \\ \frac{\partial P(\vec{x},t)}{\partial z}\end{pmatrix} +
\begin{pmatrix}
\frac{\partial^2u_x(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial z^2} \\
\frac{\partial^2u_y(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_y(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_y(\vec{x},t)}{\partial z^2} \\
\frac{\partial^2u_z(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_z(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_z(\vec{x},t)}{\partial z^2}
\end{pmatrix}
$$

$$u_x(\vec{x},t)
=-\left( \left(\frac{\partial P(\vec{x},t)}{\partial x} + \frac{\partial^2u_x(\vec{x},t)}{\partial x^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial y^2} + \frac{\partial^2u_x(\vec{x},t)}{\partial z^2} \right) \frac{1}{\rho} + \frac{\partial u_x(\vec{x},t)}{\partial t} + \frac{\partial u_x(\vec{x},t)}{\partial y}\cdot u_y(\vec{x},t) + \frac{\partial u_z(\vec{x},t)}{\partial z}\cdot u_z(\vec{x},t)\right) \left(\frac{\partial u_x(\vec{x},t)}{\partial x}\right)^{-1}$$

In this problem you are asked to calculate the $Re$ for a few biological systems. Recall that the Reynolds Number, $Re$, is

$$Re=\frac{\rho LV}{\mu}$$a)

First we consider a human embryonic heart around 3 weeks post fertilization. At this stage of development, the heart is a valveless tube, that pumps blood using peristalsis. The diameter of the heart tube is roughly $0.5\cdot10^{-4}m$ and the peak velocity of the blood is approximately $1\frac{mm}{s}$. Although physical properties of embryonic blood are not explicitly known, one tends to approximate the density as $1025\frac{kg}{m^3}$ and kinematic viscosity of $2\cdot10^{-3}Pa\cdot s$. What is the $Re$ for this cardiac system?

b)

Once the heart has fully matured to a multi-chambered valvular pumping system, it has a different $Re$. However, rather than compute the $Re$ of the adult heart, you will calculate the maximum velocity of the blood running through the aorta, e.g., the main artery of the body that is responsible for supplying oxygenated blood throughout the entire circulatory system. The adult human ascending aorta is roughly $2.3 cm$ in diameter on average. The properties of adult human blood are better documented with a density of approximately $1060 \frac{kg}{m^3}$ and kinematic viscosity of roughly $3\cdot10^{-3}Pa\cdot s$. If the $Re$ of the adult aorta is 10 000, what is the maximum velocity of the blood running through it?

c)

In the above two cases, you have calculated the either the $Re$ or parameters involved for blood flowing through different vessels. Having learned what the $Re$ was from part (a) and the maximum velocity of blood flow in part (b), what can you say about the human heart as it develops?

a) $$Re=\frac{\rho LV}{\mu} = \frac{\left(1025\frac{kg}{m^3}\right)\left(0.5\cdot10^{-4}m\right)\left(1\cdot10^{-3}\frac{m}{s}\right)}{2\cdot10^{-3}\frac{kg}{s\cdot m}}=0.025625$$

b) $$Re=\frac{\rho LV}{\mu} = \frac{\left(1060\frac{kg}{m^3}\right)\left(2.3\cdot10^{-3}m\right)\left(V\right)}{3\cdot10^{-3}\frac{kg}{s\cdot m}}=10000$$ $$V=12.3052\frac{m}{s}$$

c) The vorticity of the embryonic heart is far less than that off a fully developed one and much of this comes from the massive increase in peak velocity. The embryonic heart seems to deal with a much simpler, less turbulent fluid system.

For this problem, find a non-dimensional number in fluid dynamics (other than the Reynolds Number), and write out what parameters it depends on. You do not have to derive it. Briefly discuss what the non-dimensional number describes, and find two fluids applications it is used in (hint: *cough, cough* - Google “dimensionless fluids numbers”).

One interesting case is the Stanton Number, which governs forced convection. That is, it is used to characterize heat flow from a boundary to a fluid. It is constructed from three other nondimensional numbers: the Nusselt Number, the Prandtl Number, and the Reynolds Number.

The Nusselt Number is the ratio of convective to conductive heat transfer from a boundary to the surface of a fluid. Essentially it seems to describe how temperature differences propagate from the surface of a fluid to that fluid's body. The number is defined as $$Nu=\frac{hL}{k}=$$ $$=\frac{\text{heat transfer coefficient}\left( \frac{kg}{s^3\cdot K} \right)\text{characteristic length}\left( m \right)}{\text{thermal conductivity}\left( \frac{kg\cdot m}{s^3\cdot K} \right)}$$

The Prandtl Number is a ratio of viscous diffusion to thermal diffusion. Thus, it seems to define how much heat diffusion in a material is due to standard conduction versus how much is due to convection currents within the material. It is simply $$Pr=\frac{\nu}{\alpha}=\frac{\text{momentum diffusivity}\left( \frac{m^2}{s} \right)}{\text{thermal diffusivity}\left( \frac{m^2}{s} \right)}$$

The last number, the Reynolds Number, has already been covered exhaustively enough that I will avoid examining it in the same way.

From these three, we can construct the Stanton Number as $$St=\frac{Nu}{Pr\cdot Re}$$ What does this tell us? We seem to be dealing with a ratio between the convection of heat into a fluid and the movement of heat within that fluid (as the Prandtl Number describes the movement of heat by conduction and by current and the Reynolds number describes those currents). Taken together, this number seems to wholly capture the process by which temperature differences propagate from a boundary to within the fluid until the system reaches equilibrium.

I can see many possibly applications for this number. One case would be the pot of warming (though not boiling) water; the hot bottom and warm sides along with the cool top would create convection currents while the heat also conducts itself through the water statically. I could even imagine a more complete model that took into account the diffusion of the pot (via the heat equation) and modeled the air as a separate, large system.

The Stanton Number might also help with situations in which materials must be heated or cooled very precisely (telescope refractor mirrors come to mind; as do silicone wafers for microchip production). How best should heating elements be placed to ensure even heating? How should air currents be directed to best control cooling? How does a change in humidity or air composition affect outcomes? These are questions the Stanton Number could help answer (although exact numerical modelling is probably required, especially for those first two).

Here you are asked to run a fluid-structure interaction simulation using the immersed boundary method, and

Take a screenshot of the simulation

Describe what is happening within the simulation

This is the simulation at t=60:

From what I can gather, this example models a 2 dimensional rubber band in a fluid which is pulled by springs at the Lagrangian points of the system. I am not entirely sure *what* Lagrangian points are in the context of fluids, but in gravitational systems they are points of stability between two bodies, so perhaps they are still points in the system? In either case, this fluid seems to be under periodic and exponentially decaying forcing which creates waves of pressure; at first these are simple with only four main waves ossilating back and forth from t=0 to around t=40. This doubles to eight main waves after t=40 and twenty-four distinct waves can be seen at t=60.