*** Plane Poiseuille flow - BVP solve by shooting method
:PROPERTIES:
:categories: BVP
:date: 2013/02/15 09:00:00
:updated: 2013/03/08 10:08:42
:tags: fluids
:END:
[[http://matlab.cheme.cmu.edu/2011/09/08/plane-poiseuille-flow-bvp-solve-by-shooting-method/][Matlab post]]
One approach to solving BVPs is to use the shooting method. The reason we cannot use an initial value solver for a BVP is that there is not enough information at the initial value to start. In the shooting method, we take the function value at the initial point, and guess what the function derivatives are so that we can do an integration. If our guess was good, then the solution will go through the known second boundary point. If not, we guess again, until we get the answer we need. In this example we repeat the pressure driven flow example, but illustrate the shooting method.
In the pressure driven flow of a fluid with viscosity $\mu$ between two stationary plates separated by distance $d$ and driven by a pressure drop $\Delta P/\Delta x$, the governing equations on the velocity $u$ of the fluid are (assuming flow in the x-direction with the velocity varying only in the y-direction):
$$\frac{\Delta P}{\Delta x} = \mu \frac{d^2u}{dy^2}$$
with boundary conditions $u(y=0) = 0$ and $u(y=d) = 0$, i.e. the no-slip condition at the edges of the plate.
we convert this second order BVP to a system of ODEs by letting $u_1 = u$, $u_2 = u_1'$ and then $u_2' = u_1''$. This leads to:
$\frac{d u_1}{dy} = u_2$
$\frac{d u_2}{dy} = \frac{1}{\mu}\frac{\Delta P}{\Delta x}$
with boundary conditions $u_1(y=0) = 0$ and $u_1(y=d) = 0$.
for this problem we let the plate separation be d=0.1, the viscosity $\mu = 1$, and $\frac{\Delta P}{\Delta x} = -100$.
**** First guess
We need u_1(0) and u_2(0), but we only have u_1(0). We need to guess a value for u_2(0) and see if the solution goes through the u_2(d)=0 boundary value.
#+BEGIN_SRC python
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
d = 0.1 # plate thickness
def odefun(U, y):
u1, u2 = U
mu = 1
Pdrop = -100
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
u1_0 = 0 # known
u2_0 = 1 # guessed
dspan = np.linspace(0, d)
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0])
plt.plot([d],[0], 'ro')
plt.xlabel('d')
plt.ylabel('$u_1$')
plt.savefig('images/bvp-shooting-1.png')
#+END_SRC
#+RESULTS:
[[./images/bvp-shooting-1.png]]
Here we have undershot the boundary condition. Let us try a larger guess.
**** Second guess
:PROPERTIES:
:date: 2013/03/08 09:33:13
:updated: 2013/03/08 09:33:13
:END:
#+BEGIN_SRC python
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
d = 0.1 # plate thickness
def odefun(U, y):
u1, u2 = U
mu = 1
Pdrop = -100
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
u1_0 = 0 # known
u2_0 = 10 # guessed
dspan = np.linspace(0, d)
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0])
plt.plot([d],[0], 'ro')
plt.xlabel('d')
plt.ylabel('$u_1$')
plt.savefig('images/bvp-shooting-2.png')
#+END_SRC
#+RESULTS:
[[./images/bvp-shooting-2.png]]
Now we have clearly overshot. Let us now make a function that will iterate for us to find the right value.
**** Let fsolve do the work
#+BEGIN_SRC python
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
d = 0.1 # plate thickness
Pdrop = -100
mu = 1
def odefun(U, y):
u1, u2 = U
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
u1_0 = 0 # known
dspan = np.linspace(0, d)
def objective(u2_0):
dspan = np.linspace(0, d)
U = odeint(odefun, [u1_0, u2_0], dspan)
u1 = U[:,0]
return u1[-1]
u2_0, = fsolve(objective, 1.0)
# now solve with optimal u2_0
U = odeint(odefun, [u1_0, u2_0], dspan)
plt.plot(dspan, U[:,0], label='Numerical solution')
plt.plot([d],[0], 'ro')
# plot an analytical solution
u = -(Pdrop) * d**2 / 2 / mu * (dspan / d - (dspan / d)**2)
plt.plot(dspan, u, 'r--', label='Analytical solution')
plt.xlabel('d')
plt.ylabel('$u_1$')
plt.legend(loc='best')
plt.savefig('images/bvp-shooting-3.png')
#+END_SRC
#+RESULTS:
[[./images/bvp-shooting-3.png]]
You can see the agreement is excellent!
This also seems like a useful bit of code to not have to reinvent regularly, so it has been added to pycse as BVP_sh. Here is an example usage.
#+BEGIN_SRC python
from pycse import BVP_sh
import matplotlib.pyplot as plt
d = 0.1 # plate thickness
Pdrop = -100
mu = 1
def odefun(U, y):
u1, u2 = U
du1dy = u2
du2dy = 1.0 / mu * Pdrop
return [du1dy, du2dy]
x1 = 0.0; alpha = 0.0
x2 = 0.1; beta = 0.0
init = 2.0 # initial guess of slope at x=0
X,Y = BVP_sh(odefun, x1, x2, alpha, beta, init)
plt.plot(X, Y[:,0])
plt.ylim([0, 0.14])
# plot an analytical solution
u = -(Pdrop) * d**2 / 2 / mu * (X / d - (X / d)**2)
plt.plot(X, u, 'r--', label='Analytical solution')
plt.savefig('images/bvp-shooting-4.png')
plt.show()
#+END_SRC
#+RESULTS:
[[./images/bvp-shooting-4.png]]