In [794]:

import math
import numpy
from matplotlib import pyplot



# BEM method



In [795]:

#Q = 2000/3    #strength of the source-sheet,stb/d
h=25.26        #thickness of local gridblock,ft
phi=0.2        #porosity
kx=200         #pemerability in x direction,md
ky=200         #pemerability in y direction,md
kr=kx/ky       #pemerability ratio
miu=1          #viscosity,cp

Nw=1           #Number of well
Qwell_1=2000   #Flow rate of well 1
Boundary_V=-400 #boundary velocity ft/day



## Boundary Discretization

we will create a discretization of the body geometry into panels (line segments in 2D). A panel's attributes are: its starting point, end point and mid-point, its length and its orientation. See the following figure for the nomenclature used in the code and equations below.

Figure 1. Nomenclature of the boundary element in the local coordinates

### Create panel and well class



In [796]:

class Panel:
"""Contains information related to a panel."""
def __init__(self, xa, ya, xb, yb):
"""Creates a panel.

Arguments
---------
xa, ya -- Cartesian coordinates of the first end-point.
xb, yb -- Cartesian coordinates of the second end-point.
"""
self.xa, self.ya = xa, ya
self.xb, self.yb = xb, yb

self.xc, self.yc = (xa+xb)/2, (ya+yb)/2       # control-point (center-point)
self.length = math.sqrt((xb-xa)**2+(yb-ya)**2)     # length of the panel

# orientation of the panel (angle between x-axis and panel)
self.sinalpha=(yb-ya)/self.length
self.cosalpha=(xb-xa)/self.length

self.Q = 0.                                # source strength
self.U = 0.                                # velocity component
self.V = 0.                                # velocity component
self.P = 0.                                # pressure coefficient

class Well:
"""Contains information related to a panel."""
def __init__(self, xw, yw,rw,Q):
"""Creates a panel.

Arguments
---------
xw, yw -- Cartesian coordinates of well source.
Q      -- Flow rate of well source.
rw     -- radius of well source.
"""
self.xw, self.yw = xw, yw

self.Q = Q                               # source strength
self.rw = rw                                # velocity component



We create a node distribution on the boundary that is refined near the corner with cosspace function



In [797]:

def cosspace(st,ed,N):
N=N+1
AngleInc=numpy.pi/(N-1)
CurAngle = AngleInc
space=numpy.linspace(0,1,N)
space[0]=st
for i in range(N-1):
space[i+1] = 0.5*numpy.abs(ed-st)*(1 - math.cos(CurAngle));
CurAngle += AngleInc
if ed<st:
space[0]=ed
space=space[::-1]
return space



### Discretize boundary element along the boundary

Here we implement BEM in a squre grid



In [798]:

N=80     #Number of boundary element
Nbd=20   #Number of boundary element in each boundary
Dx=1.    #Grid block length in X direction
Dy=1.    #Gird block lenght in Y direction

#Create the array
x_ends = numpy.linspace(0, Dx, N)    # computes a 1D-array for x
y_ends = numpy.linspace(0, Dy, N)    # computes a 1D-array for y
interval=cosspace(0,Dx,Nbd)
rinterval=cosspace(Dx,0,Nbd)
#interval=numpy.linspace(0,1,Nbd+1)
#rinterval=numpy.linspace(1,0,Nbd+1)

#Define the rectangle boundary

for i in range(Nbd):
x_ends[i]=0
y_ends[i]=interval[i]

for i in range(Nbd):
x_ends[i+Nbd]=interval[i]
y_ends[i+Nbd]=Dy

for i in range(Nbd):
x_ends[i+Nbd*2]=Dx
y_ends[i+Nbd*2]=rinterval[i]

for i in range(Nbd):
x_ends[i+Nbd*3]=rinterval[i]
y_ends[i+Nbd*3]=0

x_ends,y_ends=numpy.append(x_ends, x_ends[0]), numpy.append(y_ends, y_ends[0])

#Define the panel
panels = numpy.empty(N, dtype=object)
for i in range(N):
panels[i] = Panel(x_ends[i], y_ends[i], x_ends[i+1], y_ends[i+1])

#Define the well
wells = numpy.empty(Nw, dtype=object)

wells[0]=Well(Dx/2,Dy/2,0.025,Qwell_1)




In [799]:

#for i in range(N):
#   print("Panel Coordinate (%s,%s) sina,cosa (%s,%s) " % (panels[i].xc,panels[i].yc,panels[i].sinalpha,panels[i].cosalpha))
#print("Well Location (%s,%s) radius: %s Flow rate:%s  " % (wells[0].xw,wells[0].yw,wells[0].rw,wells[0].Q))



### Plot boundary elements and wells



In [800]:

#Plot the panel
%matplotlib inline

val_x, val_y = 0.3, 0.3
x_min, x_max = min(panel.xa for panel in panels), max(panel.xa for panel in panels)
y_min, y_max = min(panel.ya for panel in panels), max(panel.ya for panel in panels)
x_start, x_end = x_min-val_x*(x_max-x_min), x_max+val_x*(x_max-x_min)
y_start, y_end = y_min-val_y*(y_max-y_min), y_max+val_y*(y_max-y_min)

size = 5
pyplot.figure(figsize=(size, (y_end-y_start)/(x_end-x_start)*size))
pyplot.grid(True)
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)

pyplot.plot(numpy.append([panel.xa for panel in panels], panels[0].xa),
numpy.append([panel.ya for panel in panels], panels[0].ya),
linestyle='-', linewidth=1, marker='o', markersize=6, color='#CD2305');
pyplot.scatter(wells[0].xw,wells[0].yw,s=100,alpha=0.5)

pyplot.legend(['panels', 'Wells'],
loc=1, prop={'size':12})




Out[800]:

<matplotlib.legend.Legend at 0x19cb8278>



## Boundary element implementation

Figure 2. Representation of a local gridblock with boundary elements

Generally, the influence of all the j panels on the i BE node can be expressed as follows:

\begin{matrix} {{c}_{ij}}{{p}_{i}}+{{p}_{i}}\int_{{{s}_{j}}}{{{H}_{ij}}d{{s}_{j}}}=({{v}_{i}}\cdot \mathbf{n})\int_{{{s}_{j}}}{{{G}_{ij}}}d{{s}_{j}} \end{matrix}

Where,

${{c}_{ij}}$ is the free term, cased by source position.

${{c}_{ij}}=\left\{ \begin{matrix} \begin{matrix} 1 & \text{source j on the internal domain} \\ \end{matrix} \\ \begin{matrix} 0.5 & \text{source j on the boundary} \\ \end{matrix} \\ \begin{matrix} 0 & \text{source j on the external domain} \\ \end{matrix} \\ \end{matrix} \right.$

$\int_{{{s}_{j}}}{{{H}_{ij}}d{{s}_{j}}\text{ }}$ is the integrated effect of the boundary element source i on the resulting normal flux at BE node j.

$\int_{{{s}_{j}}}{{{G}_{ij}}}d{{s}_{j}}$ is the is the integrated effect of the boundary element source i on the resulting pressure at BE node j

### Line segment source solution for pressure and velocity (Derived recently)

The integrated effect can be formulated using line segment source solution, which givs:

$$\int_{{{s}_{j}}}{{{G}_{ij}}}d{{s}_{j}}=B{{Q}_{w}}=P({{{x}'}_{i}},{{{y}'}_{i}})=-\frac{70.60\mu }{h\sqrt{{{k}_{x}}{{k}_{y}}}}\int_{t=0}^{t={{l}_{j}}}{\ln \left\{ {{({x}'-t\cos {{\alpha }_{j}})}^{2}}+\frac{{{k}_{x}}}{{{k}_{y}}}{{({y}'-t\sin {{\alpha }_{j}})}^{2}} \right\}dt}\cdot {{Q}_{w}}$$$$\int_{{{s}_{j}}}{{{H}_{ij}}d{{s}_{j}}\text{ }}={{v}_{i}}(s)\cdot {{\mathbf{n}}_{i}}=-{{u}_{i}}\sin {{\alpha }_{i}}+{{v}_{i}}\cos {{\alpha }_{i}}$$

Where,

$$u\left( {{{{x}'}}_{i}},{{{{y}'}}_{i}} \right)={{A}_{u}}{{Q}_{j}}=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}_{x}}}{{{k}_{y}}}}\int_{t=0}^{t={{l}_{j}}}{\frac{{{{{x}'}}_{i}}-t\cos {{\alpha }_{j}}}{{{\left( {{{{x}'}}_{i}}-t\cos {{\alpha }_{j}} \right)}^{2}}+\frac{{{k}_{x}}}{{{k}_{y}}}{{({{{{y}'}}_{i}}-t\sin {{\alpha }_{j}})}^{2}}}dt}\cdot {{Q}_{j}}$$$$v\left( {{{{x}'}}_{i}},{{{{y}'}}_{i}} \right)={{A}_{v}}{{Q}_{j}}=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}_{x}}}{{{k}_{y}}}}\int_{t=0}^{t={{l}_{j}}}{\frac{{{{{y}'}}_{i}}-t\sin {{\alpha }_{j}}}{{{\left( {{{{x}'}}_{i}}-t\cos {{\alpha }_{j}} \right)}^{2}}+\frac{{{k}_{x}}}{{{k}_{y}}}{{({{{{y}'}}_{i}}-t\sin {{\alpha }_{j}})}^{2}}}dt}\cdot {{Q}_{j}}$$

## Line segment source Integration function (Bij and Aij)



In [801]:

#Panel infuence factor Bij
def InflueceP(x, y, panel):
"""Evaluates the contribution of a panel at one point.

Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.

Returns
-------
Integral over the panel of the influence at one point.
"""
#Transfer global coordinate point(x,y) to local coordinate
x=x-panel.xa
y=y-panel.ya
L1=panel.length

#Calculate the pressure and velocity influence factor
a=panel.cosalpha**2+kr*panel.sinalpha**2
b=x*panel.cosalpha+kr*panel.sinalpha*y
c=y*panel.cosalpha-x*panel.sinalpha
dp=70.6*miu/h/math.sqrt(kx*ky)
Cp = dp/a*(
(
b*math.log(x**2-2*b*L1+a*L1**2+kr*y**2)
-L1*a*math.log((x-L1*panel.cosalpha)**2+kr*(y-L1*panel.sinalpha)**2)
+2*math.sqrt(kr)*c*math.atan((b-a*L1)/math.sqrt(kr)/c)
)
-
(
b*math.log(x**2+kr*y**2)
+2*math.sqrt(kr)*c*math.atan((b)/math.sqrt(kr)/c)
)
)
#debug
#print("a: %s b:%s c:%s  " % (a,b,c))
#angle=math.atan((b-a*L1)/math.sqrt(kr)/c)*180/numpy.pi
#print("Magic angle:%s"% angle)
return Cp

def InflueceU(x, y, panel):
"""Evaluates the contribution of a panel at one point.

Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.

Returns
-------
Integral over the panel of the influence at one point.
"""
#Transfer global coordinate point(x,y) to local coordinate
x=x-panel.xa
y=y-panel.ya
L1=panel.length

#Calculate the pressure and velocity influence factor
a=panel.cosalpha**2+kr*panel.sinalpha**2
b=x*panel.cosalpha+kr*panel.sinalpha*y
c=y*panel.cosalpha-x*panel.sinalpha
dv=-0.4468/h/phi*math.sqrt(kx/ky)
Cu = dv/a*(
(
panel.cosalpha*math.log(x**2-2*b*L1+a*L1**2+kr*y**2)+ 2*math.sqrt(kr)*panel.sinalpha*math.atan((a*L1-b)/math.sqrt(kr)/c)
)
-
(
panel.cosalpha*math.log(x**2+kr*y**2)+2*math.sqrt(kr)*panel.sinalpha*math.atan((-b)/math.sqrt(kr)/c)
)
)
#print("a: %s b:%s c:%s  " % (a,b,c))
#angle=math.atan((b-a*L1)/math.sqrt(kr)/c)*180/numpy.pi
#print("Magic angle:%s"% angle)
return Cu

def InflueceV(x, y, panel):
"""Evaluates the contribution of a panel at one point.

Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.

Returns
-------
Integral over the panel of the influence at one point.
"""
#Transfer global coordinate point(x,y) to local coordinate
x=x-panel.xa
y=y-panel.ya
L1=panel.length

#Calculate the pressure and velocity influence factor
a=panel.cosalpha**2+kr*panel.sinalpha**2
b=x*panel.cosalpha+kr*panel.sinalpha*y
c=y*panel.cosalpha-x*panel.sinalpha
dv=-0.4468/h/phi*math.sqrt(kx/ky)
Cv = dv/a*(
(
panel.sinalpha*math.log(x**2-2*b*L1+a*L1**2+kr*y**2)+ 2*math.sqrt(1/kr)*panel.cosalpha*math.atan((b-a*L1)/math.sqrt(kr)/c)
)
-
(
panel.sinalpha*math.log(x**2+kr*y**2)+2*math.sqrt(1/kr)*panel.cosalpha*math.atan((b)/math.sqrt(kr)/c)
)
)
#print("a: %s b:%s c:%s  " % (a,b,c))
#angle=math.atan((b-a*L1)/math.sqrt(kr)/c)*180/numpy.pi
#print("Magic angle:%s"% angle)

return Cv



## Well source function

### Line source solution for pressure and velocity (Datta-Gupta, 2007)

$$P(x,y)=B{{Q}_{w}}=-\frac{70.60\mu }{h\sqrt{{{k}_{x}}{{k}_{y}}}}\ln \left\{ {{(x-{{x}_{w}})}^{2}}+\frac{{{k}_{x}}}{{{k}_{y}}}{{(y-{{y}_{w}})}^{2}} \right\}{{Q}_{w}}+{{P}_{avg}}$$$$\frac{\partial P}{\partial x}=u=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}_{x}}}{{{k}_{y}}}}\sum\limits_{k=1}^{{{N}_{w}}}{{{Q}_{k}}}\frac{x-{{x}_{k}}}{{{\left( x-{{x}_{k}} \right)}^{2}}+\frac{{{k}_{x}}}{{{k}_{y}}}{{(y-{{y}_{k}})}^{2}}}$$$$\frac{\partial P}{\partial y}=v=\frac{0.8936}{h\phi }\sqrt{\frac{{{k}_{x}}}{{{k}_{y}}}}\sum\limits_{k=1}^{{{N}_{w}}}{{{Q}_{k}}}\frac{y-{{y}_{k}}}{{{\left( x-{{x}_{k}} \right)}^{2}}+\frac{{{k}_{x}}}{{{k}_{y}}}{{(y-{{y}_{k}})}^{2}}}$$


In [802]:

#Well influence factor
def InflueceP_W(x, y, well):
"""Evaluates the contribution of a panel at one point.

Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.

Returns
-------
Integral over the panel of the influence at one point.
"""
dp=-70.6*miu/h/math.sqrt(kx*ky)
Cp=dp*math.log((x-well.xw)**2+kr*(y-well.yw)**2)
return Cp

def InflueceU_W(x, y, well):
"""Evaluates the contribution of a panel at one point.

Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.

Returns
-------
Integral over the panel of the influence at one point.
"""
dv=0.8936/h/phi*math.sqrt(kx/ky)
Cu=dv*(x-well.xw)/((x-well.xw)**2+kr*(y-well.yw)**2)
return Cu

def InflueceV_W(x, y, well):
"""Evaluates the contribution of a panel at one point.

Arguments
---------
x, y -- Cartesian coordinates of the point.
panel -- panel which contribution is evaluated.

Returns
-------
Integral over the panel of the influence at one point.
"""
dv=0.8936/h/phi*math.sqrt(kx/ky)
Cv=dv*(y-well.yw)/((x-well.xw)**2+kr*(y-well.yw)**2)
return Cv




In [803]:

#InflueceV(0.5,1,panels[3])
#InflueceP(0,0.5,panels[0])
#InflueceU(0,0.5,panels[0])



## BEM function solution

Generally, the influence of all the j panels on the i BE node can be expressed as follows:

\begin{matrix} {{c}_{ij}}{{p}_{i}}+{{p}_{i}}\int_{{{s}_{j}}}{{{H}_{ij}}d{{s}_{j}}}=({{v}_{i}}\cdot \mathbf{n})\int_{{{s}_{j}}}{{{G}_{ij}}}d{{s}_{j}} \end{matrix}

Applying boundary condition along the boundary on above equation, a linear systsem can be constructed as follows:

\begin{matrix} \left[ {{{{H}'}}_{ij}} \right]\left[ {{P}_{i}} \right]=\left[ {{G}_{ij}} \right]\left[ {{v}_{i}}\cdot \mathbf{n} \right] \end{matrix}

!!!!!MY IMPLEMENTATION MAY HAS SOME PROBLEM HERE!!!!!!

All the integration solution can be evaluated except on itself. Where,

$\left[ {{{{H}'}}_{ij}} \right]=\left\{ \begin{matrix} \begin{matrix} {{H}_{ij}} & i\ne j \\ \end{matrix} \\ \begin{matrix} {{H}_{ij}}+\frac{1}{2} & i=j \\ \end{matrix} \\ \end{matrix} \right.$

Figure 3. Representation of coordinate systems and the principle of superstition with well source and boundary element source

As shown in Fig.3, the pressure and velocity at any point i in the local gridblock can be determined using Eqs. below. Applying principle of superposition for each BE node along the boundary (Fig. 3), boundary condition can be written as follows:

\begin{matrix} {{P}_{i}}(s)=\sum\limits_{j=1}^{M}{{{B}_{ij}}{{Q}_{j}}} & \text{constant pressure boundary} \\ \end{matrix}\begin{matrix} {{v}_{i}}(s)\cdot {{\mathbf{n}}_{i}}=\sum\limits_{j=1}^{M}{{{A}_{ij}}{{Q}_{j}}} & \text{constant flux boundary} \\ \end{matrix}

The Pi and v ·n are the konwn boundary codition. The flow rate(strength) of boundary elements in Hij and Gij are the only unknown terms. So we could rearrange the matrix above as linear system:

${{\left[ \begin{matrix} {{A}_{ij}} \\ {{B}_{ij}} \\ \end{matrix} \right]}_{N\times N}}{{\left[ \begin{matrix} {{Q}_{j}} \\ {{Q}_{j}} \\ \end{matrix} \right]}_{N\times 1}}={{\left[ \begin{matrix} -{{u}_{i}}\sin {{\alpha }_{i}}+{{v}_{i}}\cos {{\alpha }_{i}} \\ {{P}_{i}} \\ \end{matrix} \right]}_{N\times 1}}$


In [804]:

def build_matrix(panels):
"""Builds the source matrix.

Arguments
---------
panels -- array of panels.

Returns
-------
A -- NxN matrix (N is the number of panels).
"""
N = len(panels)
A = numpy.empty((N, N), dtype=float)
#numpy.fill_diagonal(A, 0.5)

for i, p_i in enumerate(panels): #target nodes
for j, p_j in enumerate(panels): #BE source
#if i != j:   ###Matrix construction
if i>=0 and i<Nbd or i>=3*Nbd and i<4*Nbd:
A[i,j] = -p_j.sinalpha*InflueceU(p_i.xc, p_i.yc, p_j)+p_j.cosalpha*InflueceV(p_i.xc, p_i.yc, p_j)
#A[i,j] = InflueceP(p_i.xc, p_i.yc, p_j)
if i>=Nbd and i<2*Nbd or i>=2*Nbd and i<3*Nbd:
A[i,j] = -p_j.sinalpha*InflueceU(p_i.xc, p_i.yc, p_j)+p_j.cosalpha*InflueceV(p_i.xc, p_i.yc, p_j)
#A[i,j] = InflueceP(p_i.xc, p_i.yc, p_j)

return A

def build_rhs(panels):
"""Builds the RHS of the linear system.

Arguments
---------
panels -- array of panels.

Returns
-------
b -- 1D array ((N+1)x1, N is the number of panels).
"""
b = numpy.empty(len(panels), dtype=float)

for i, panel in enumerate(panels):
V_well=( -panel.sinalpha*Qwell_1*InflueceU_W(panel.xc, panel.yc, wells[0])+panel.cosalpha*Qwell_1*InflueceV_W(panel.xc, panel.yc, wells[0]) )
if i>=0 and i<Nbd:
b[i]=0+V_well
#b[i]=4000
#b[i]=84
if i>=Nbd and i<2*Nbd:
b[i]=-V_well
#b[i]=-42
if i>=2*Nbd and i<3*Nbd:
b[i]=-V_well
#b[i]=-42
if i>=3*Nbd and i<4*Nbd:
b[i]=0+V_well
#b[i]=84
return b




In [805]:

#Qwell_1=300   #Flow rate of well 1
#Boundary_V=-227 #boundary velocity ft/day

A = build_matrix(panels)                    # computes the singularity matrix
b = build_rhs(panels)                       # computes the freestream RHS




d:\Anaconda3\lib\site-packages\ipykernel\__main__.py:66: RuntimeWarning: divide by zero encountered in double_scalars
d:\Anaconda3\lib\site-packages\ipykernel\__main__.py:70: RuntimeWarning: divide by zero encountered in double_scalars
d:\Anaconda3\lib\site-packages\ipykernel\__main__.py:102: RuntimeWarning: divide by zero encountered in double_scalars
d:\Anaconda3\lib\site-packages\ipykernel\__main__.py:106: RuntimeWarning: divide by zero encountered in double_scalars




In [806]:

# solves the linear system
Q = numpy.linalg.solve(A, b)

for i, panel in enumerate(panels):
panel.Q = Q[i]



# Plot results



In [807]:

#Visulize the pressure and velocity field

#Define meshgrid
Nx, Ny = 50, 50                  # number of points in the x and y directions
x_start, x_end = -0.01, 1.01            # x-direction boundaries
y_start, y_end = -0.01, 1.01            # y-direction boundaries
x = numpy.linspace(x_start, x_end, Nx)    # computes a 1D-array for x
y = numpy.linspace(y_start, y_end, Ny)    # computes a 1D-array for y
X, Y = numpy.meshgrid(x, y)               # generates a mesh grid

#Calculate the velocity and pressure field
p = numpy.empty((Nx, Ny), dtype=float)
u = numpy.empty((Nx, Ny), dtype=float)
v = numpy.empty((Nx, Ny), dtype=float)

#for i, panel in enumerate(panels):
#panel.Q = 0.

#panels[0].Q=100
#panels[5].Q=100
#Qwell_1=400

for i in range(Nx):
for j in range(Ny):
p[i,j] =sum([p.Q*InflueceP(X[i,j], Y[i,j], p) for p in panels])+Qwell_1*InflueceP_W(X[i,j], Y[i,j], wells[0])
u[i,j] =sum([p.Q*InflueceU(X[i,j], Y[i,j], p) for p in panels])+Qwell_1*InflueceU_W(X[i,j], Y[i,j], wells[0])
v[i,j] =sum([p.Q*InflueceV(X[i,j], Y[i,j], p) for p in panels])+Qwell_1*InflueceV_W(X[i,j], Y[i,j], wells[0])
#p[i,j] =sum([p.Q*InflueceP(X[i,j], Y[i,j], p) for p in panels])
#u[i,j] =sum([p.Q*InflueceU(X[i,j], Y[i,j], p) for p in panels])
#v[i,j] =sum([p.Q*InflueceV(X[i,j], Y[i,j], p) for p in panels])
#p[i,j] =Qwell_1*InflueceP_W(X[i,j], Y[i,j], wells[0])
#u[i,j] =Qwell_1*InflueceU_W(X[i,j], Y[i,j], wells[0])
#v[i,j] =Qwell_1*InflueceV_W(X[i,j], Y[i,j], wells[0])




In [808]:

# plots the streamlines
%matplotlib inline

size = 6
pyplot.figure(figsize=(size, size))
pyplot.grid(True)
pyplot.title('Streamline field')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(-0.2, 1.2)
pyplot.ylim(-0.2, 1.2)

pyplot.plot(numpy.append([panel.xa for panel in panels], panels[0].xa),
numpy.append([panel.ya for panel in panels], panels[0].ya),
linestyle='-', linewidth=1, marker='o', markersize=6, color='#CD2305');
stream =pyplot.streamplot(X, Y, u, v,density=2, linewidth=1, arrowsize=1, arrowstyle='->') #streamline
#cbar=pyplot.colorbar(orientation='vertical')

#equipotential=pyplot.contourf(X, Y, p1, extend='both')







In [809]:

size = 7
pyplot.figure(figsize=(size, size-1))
pyplot.title('Pressure field')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(0, 1)
pyplot.ylim(0, 1)

pyplot.contour(X, Y, p, 15, linewidths=0.5, colors='k')
pyplot.contourf(X, Y, p, 15, cmap='rainbow',
vmax=abs(p).max(), vmin=-abs(p).max())
pyplot.colorbar()  # draw colorbar




Out[809]:

<matplotlib.colorbar.Colorbar at 0x22b3b5f8>




In [810]:

size = 7
pyplot.figure(figsize=(size, size-1))
pyplot.title('Total Velocity field')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(0, 1)
pyplot.ylim(0, 1)

Vtotal= numpy.sqrt(u**2+v**2)
#Vtotal= numpy.abs(v)
pyplot.contour(X, Y, Vtotal, 15, linewidths=0.5, colors='k')
pyplot.contourf(X, Y, Vtotal, 15, cmap='rainbow')
#vmax=50, vmin=0)
pyplot.colorbar()  # draw colorbar




Out[810]:

<matplotlib.colorbar.Colorbar at 0x1c9c09b0>




In [811]:

pyplot.title('Darcy velocity on the outflow boundary, x component (ft/day)')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)

pyplot.plot(y, u[49,:], '--', linewidth=2)
pyplot.plot(9.8425+y, u[:,49], '--', linewidth=2)
u[:,49]




Out[811]:

array([ 3991.53526449,  4697.23214482,  4427.94528366,  4213.10279637,
4056.60881842,  3935.0524005 ,  3835.61636196,  3758.7042195 ,
3678.83215419,  3624.28803091,  3574.58365565,  3512.86339586,
3477.97481989,  3446.47687207,  3393.90113729,  3365.99608357,
3346.67520352,  3314.63917316,  3278.95525446,  3263.89905598,
3249.45519765,  3219.13221798,  3199.8855574 ,  3188.04140691,
3175.46727089,  3157.2693224 ,  3144.22991707,  3132.72142655,
3123.92635749,  3118.96690841,  3108.32799206,  3098.56287193,
3101.19357179,  3101.62102519,  3094.55662065,  3093.65095622,
3110.17941018,  3111.68331515,  3115.93189462,  3141.31882334,
3156.43775968,  3175.3078    ,  3215.96547724,  3250.99279038,
3305.71980505,  3375.05945207,  3470.19304338,  3601.64257033,
3743.87730196,  3061.23226415])




In [812]:

pyplot.title('Darcy velocity on the outflow boundary, y component (ft/day)')

pyplot.plot(y, v[:,49], '--', linewidth=2)
pyplot.plot(9.8425+y, v[49,:], '--', linewidth=2)
v[49,:]




Out[812]:

array([ 3991.53526449,  4697.23214482,  4427.94528366,  4213.10279637,
4056.60881842,  3935.0524005 ,  3835.61636196,  3758.7042195 ,
3678.83215419,  3624.28803091,  3574.58365565,  3512.86339586,
3477.97481989,  3446.47687207,  3393.90113729,  3365.99608357,
3346.67520352,  3314.63917316,  3278.95525446,  3263.89905598,
3249.45519765,  3219.13221798,  3199.8855574 ,  3188.04140691,
3175.46727089,  3157.2693224 ,  3144.22991707,  3132.72142655,
3123.92635749,  3118.96690841,  3108.32799206,  3098.56287193,
3101.19357179,  3101.62102519,  3094.55662065,  3093.65095622,
3110.17941018,  3111.68331515,  3115.93189462,  3141.31882334,
3156.43775968,  3175.3078    ,  3215.96547724,  3250.99279038,
3305.71980505,  3375.05945207,  3470.19304338,  3601.64257033,
3743.87730196,  3061.23226415])




In [ ]:




In [ ]: