from IPython.display import display, Image
from IPython.display import display_pretty, display_html, display_jpeg, display_png, display_json, display_latex, display_svg
from IPython.display import Latex
from IPython.display import HTML
Image(filename='logo3.jpg', width=100)
For $U$ and $V$ we should use Dirichlet boundary conditions, which are the simplest ones fixing boundary point value for our region. Of course in case of North and South boundaries for $U$ and West and East for $V$ we shall average two neibghouring points in the following manner:
In case of implementation we would like to define initial conditions before entering time loop as set of 4x2 values representing 4 directions and two vector fields in following way:
uN = np.array(x*0)
vN = np.array(avg(x)*0)
uS = np.array(x*0)
vS = np.array(avg(x)*0)
uW = np.array(avg(y)*0+1.0)
vW = np.array(y*0)
uE = np.array(avg(y)*0+1.0)
vE = np.array(y*0)
It would mean we have constant West -> East flow
Those arrays we are going to concatenate with U and V interior arrays in first step of each iteration according to formula 1.2 as follows:
Ue = op.cncat((uW, U, uE))
Ue = op.cncat((2 * op.vec2col(uS) - Ue[:, 0], Ue, 2 * op.vec2col(uN)-Ue[:, -1]), dim='h')
Ve = op.cncat((vS.T, V, vN.T), dim='h')
Ve = op.cncat((2 * op.vec2row(vW) - Ve[1, :], Ve, op.vec2row(vE) - Ve[:, -1]))
We would like to have constant stream through boundaries so Neumann boudary conditions seems like a perfect match. Namely we would like to have zero flow, which means following equation holds:
For example for y-axis we have: $$\frac {P_{i,j+1} - P_{i,j}}{h} = 0 \quad Thus \quad P_{i,j+1} = P_{i,j}$$
In case of Laplacians there is one decisive paremeter for inserting boudary conditions, namely: $$ -L_0 = \frac{1}{h^2}\begin{bmatrix} a_{11} & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & a_{11} \end{bmatrix} \quad(1.5) $$
$a_{11}$ | Boundary condition |
---|---|
1 | Neumann |
2 | Dirichlet |
3 | Dirichlet between points |
If we act with this operator on matrix of coordinates $X$ you should see why:
Usually boudary conditions have a form: $$L_u X = Constant$$ $$$$ $$\frac {a_{11}X_{11} + X_{21}}{h^2} = C$$ $$$$ $$a_{11} = -1 => C = \frac {X_{21} - X_{11}}{h} = \frac {dX}{dx}$$
Same story happens for $a_{11} = -2$
?
Which means pressure Laplacian is equal to:
Lp = sparse.kron(sparse.eye(ny), Lap(nx, hx, boudary_cond='Neumann')) + \
sparse.kron(Lap(ny, hy, boudary_cond='Neumann'), sparse.eye(nx))
We are going to solve two types of linear equations:
For each velocity component and Poisson equation for pressure:
Here Python flexibility plays important role, to solve those linear expression we might use library associeted with scipy, usefull for solving sparse matrix equations:
from scipy.sparse.linalg import spsolve
As I mentioned before we want to obtain optimum speed of calculation, so we are going to use vector form of $U$ and $V$. It means our implementation of equations should look like this:
def generate_laplacians(ubx='Neumann', uby="Neumann", vbx="Neumann", vby="Neumann", pbx="Neumann", pby="Neumann",
qbx="Dirichlet_between", qby="Dirichlet_between"):
"""Returns laplacians for vector field (U, V) and scalar fields P and stream function Q
Arguments are just names of boundary conditions
"""
Lu = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny), Lap(nx-1, hx, boudary_cond=ubx)) +
sparse.kron(Lap(ny, hy, boudary_cond=uby), sparse.eye(nx-1)))
Lv = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny-1), Lap(nx, hx, boudary_cond=vbx)) +
sparse.kron(Lap(ny-1, hy, boudary_cond=vby), sparse.eye(nx)))
Lp = sparse.kron(sparse.eye(ny), Lap(nx, hx, boudary_cond=pbx)) + \
sparse.kron(Lap(ny, hy, boudary_cond=pby), sparse.eye(nx))
Lq = sparse.kron(sparse.eye(ny-1), Lap(nx-1, hx, boudary_cond=qbx)) + \
sparse.kron(Lap(ny-1, hy, boudary_cond=qby), sparse.eye(nx-1))
return Lu, Lv, Lp, Lq
#velocity:
rhs = np.ravel(U)
U = spsolve(Lu, rhs, permc_spec='COLAMD') #colamd is one of preparation method used in spsolve
U = np.reshape(U, (nx-1, ny))
rhs = np.ravel(V)
V = spsolve(Lv, rhs, permc_spec='COLAMD')
V = np.reshape(V, (nx, ny-1))
#pressure:
rhs = np.ravel((np.diff(op.cncat((uW, U, uE)), axis=0)/hx + np.diff(op.cncat((vS.T, V, vN.T), dim='h'))/hy)) #gradient
P = spsolve(Lp, -rhs, permc_spec='COLAMD')
P = np.reshape(P[:], (nx, ny)) # 2.2a
U = U - np.diff(P, axis=0)/hx # 2.2b
V = V - np.diff(P)/hy
Where in formulas 2.2a 2.2b $\Delta t$ canceled out.
Real - time visualisation is possible due to incompressible flow argument. It means divergence of vector field $V$, $U$ is 0 so in analogy with electrodynamics we can introduce vector potential function call stream function $Q$ which is defined as follows:
$$\frac {\partial Q}{\partial x} = u \quad \frac {\partial Q}{\partial y} = -v$$
There is a operation called orthogonal gradient which is defined to match equations above, using it we find:
$$(\nabla Q)^o = \bigg[\frac{\partial Q}{\partial x},-\frac{\partial Q}{\partial y}\bigg] = {\bf U}$$
Applying curl on both sides yells:
$$\nabla \times (\nabla Q)^o = -\Delta Q^o = {\nabla \times{ \bf U} = \frac {\partial u}{\partial y} - \frac {\partial v}{\partial x}}$$
Which means calculating stream function Q turns into Poisson Equation. We are assuming Q as zero at the boundaries so final implementation has a following form:
rhs = np.diff(U)/hy - np.diff(V, axis=0)/hx
rhs = np.ravel(rhs)
q = spsolve(Lq, rhs, permc_spec='COLAMD')
Q = np.zeros((nx+1, ny+1))
Q[1:-1, 1:-1] = np.reshape(q, (nx-1, ny-1))
Visualisation should not be difficult! For the last time we connect interior points of velocity field using quiver plot which combinated with counter plot of stream function and coulored countur plot of pressure should result in beautifull animation.
cxf = plt.contourf(x_avg, y_avg, P.T)
cx = plt.contour(X, Y, Q.T)
Ue = op.cncat((op.vec2col(uS), avg(op.cncat((uW, U, uE)).T).T, op.vec2col(uN)), dim='h')
Ve = np.concatenate((op.vec2row(vW), avg(op.cncat((op.vec2col(vS), V, op.vec2col(vN)), dim='h')), op.vec2row(vE)))
Len = (Ue**2 + Ve**2 + eps)**0.5 #normalising factor
ax = plt.quiver(x, y, (Ue/Len).T, (Ve/Len).T, scale=50, headwidth=4, headlength=3, headaxislength=3)#normalizing factor
plt.draw()
plt.clf()