In [1]:
import math
import numpy
from matplotlib import pyplot
In [2]:
#Q = 2000/3 #strength of the source-sheet,stb/d
Dx=9.8425
Dy=9.8425
h=26.25 #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
Nimg=40 #Number of image well
rab=0.00003 #distance betwwen node pairs
Qwell_1=2000 #Flow rate of well 1
In [3]:
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.
x_a,y_a -- node pair A, inner point
y_b,
"""
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
#unit vector (x2-x1)/length (y2-y1)/length
#normal unit vector (y1-y2)/length (x2-x1)/length
#point with a given distance along a line: x?=x_Start+unit_vector_x*distance y?=y_Start+unit_vector_y*distance
self.x_a=self.xc-rab/2*(ya-yb)/self.length
self.y_a=self.yc-rab/2*(xb-xa)/self.length
self.x_b=self.xc+rab/2*(ya-yb)/self.length
self.y_b=self.yc+rab/2*(xb-xa)/self.length
# 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 # well radius
In [4]:
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
In [5]:
#Create the array
N=80
Nbd=20
x_ends = numpy.linspace(0, 1, N) # computes a 1D-array for x
x_ends_img= numpy.linspace(0, 1, N)
y_ends = numpy.linspace(0, 1, N) # computes a 1D-array for y
y_ends_img= numpy.linspace(0, 1, N)
interval=cosspace(0,Dx,Nbd)
rinterval=cosspace(Dx,0,Nbd)
#interval=numpy.linspace(0,Dx,Nbd+1)
#rinterval=numpy.linspace(Dx,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/node pairs
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 boundary conditions
for i in range(Nbd):
panels[i].Q=0.
panels[i+Nbd].Q=1000.
panels[i+Nbd*2].Q=1000.
panels[i+Nbd*3].Q=0.
# Generalize the image well position/Define the image panel around the computation domain
R=Dx
x_ends = Dx/2+R*numpy.cos(numpy.linspace(0, 2*math.pi, Nimg+1))
y_ends = Dy/2+R*numpy.sin(numpy.linspace(0, 2*math.pi, Nimg+1))
panels_img = numpy.empty(Nimg, dtype=object)
for i in range(Nimg):
panels_img[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 [6]:
#for i in range(N):
#print("%s Panel Coordinate (%s,%s) sina,cosa (%s,%s) Q:%s" % (i+1,panels[i].xc,panels[i].yc,panels[i].sinalpha,panels[i].cosalpha,panels[i].Q))
# print("%s NodePair Coordinate A(%s,%s) B(%s,%s) Q:%s" % (i+1,panels[i].x_a,panels[i].y_a,panels[i].x_b,panels[i].y_b,panels[i].Q))
#print("Well Location (%s,%s) radius: %s Flow rate:%s " % (wells[0].xw,wells[0].yw,wells[0].rw,wells[0].Q))
In [19]:
#Plot the panel
%matplotlib inline
val_x, val_y = 1.2, 1.2
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 = 7
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=0, color='#CD2305');
pyplot.plot(numpy.append([panel.xa for panel in panels_img], panels_img[0].xa),
numpy.append([panel.ya for panel in panels_img], panels_img[0].ya),
linestyle='-', linewidth=1, marker='o', markersize=6, color='r');
#pyplot.scatter([w.xw for w in Image_wells], [w.yw for w in Image_wells], color='#CD2305', s=40)
pyplot.scatter([p.x_a for p in panels], [p.y_a for p in panels], color='b', s=40)
pyplot.scatter([p.x_b for p in panels], [p.y_b for p in panels], color='b', s=40)
pyplot.scatter(wells[0].xw,wells[0].yw,s=100,alpha=0.5)
pyplot.legend(['Interest domain','panels','Boundary node pair A','Boundary node pair B','Wells'],
loc=1, prop={'size':10})
Out[19]:
In [497]:
#Panel infuence factor
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 = -1/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
In [498]:
#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)
Cp=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 [499]:
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, Nimg), dtype=float)
#numpy.fill_diagonal(A, 0.5)
Coeff=-70.6*miu/h/math.sqrt(kx*ky)
for i, p_i in enumerate(panels): #target nodes
for j, p_j in enumerate(panels_img): #BE source
#if i>=0 and i<Nbd or i>=3*Nbd and i<4*Nbd:
A[i,j] = InflueceP(p_i.x_a,p_i.y_a,p_j)-InflueceP(p_i.x_b,p_i.y_b,p_j)
#A[i,j] /= Coeff
#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, p_i in enumerate(panels):
#b[i]=Qwell_1*( InflueceP_W(p_i.x_b,p_i.y_b,wells[0])-InflueceP_W(p_i.x_a,p_i.y_a,wells[0]) )+887.1857*p_i.Q*miu*rab/Dx/h/kx
#b[i] /= Coeff
b[i]=Qwell_1*( InflueceP_W(p_i.x_b,p_i.y_b,wells[0])-InflueceP_W(p_i.x_a,p_i.y_a,wells[0]) )-4*numpy.pi*p_i.Q*rab*math.sqrt(kx*ky)/Dx/kx
#b[i]=887.2*p_i.Q*miu*rab/Dx/h
return b
In [500]:
A = build_matrix(panels) # computes the singularity matrix
b = build_rhs(panels) # computes the freestream RHS
In [501]:
4*numpy.pi*1000*rab*math.sqrt(kx*ky)/Dx/kx
Out[501]:
In [502]:
70.6*miu/h/math.sqrt(kx*ky)
Out[502]:
In [503]:
b
Out[503]:
In [504]:
A
Out[504]:
In [505]:
# solves the linear system
Q = numpy.linalg.lstsq(A, b)
for i in range(Nimg):
panels_img[i].Q =Q[0][i]
In [506]:
Q[0]
Out[506]:
In [507]:
#for i in range(Nimg):
#print("%s Image well flow rate:%s" %(i+1,Image_wells[i].Q))
In [508]:
#Visulize the pressure and velocity field
#Define meshgrid
Nx, Ny = 50, 50 # number of points in the x and y directions
val_x, val_y = 0, 0
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)
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)
Coeff=-70.6*miu/h/math.sqrt(kx*ky)
for i in range(Nx):
for j in range(Ny):
p[i,j] =sum([w.Q*Coeff*InflueceP(X[i,j], Y[i,j], w) for w in panels_img])+Qwell_1*Coeff*InflueceP_W(X[i,j], Y[i,j], wells[0])-160
u[i,j] =sum([w.Q*InflueceU(X[i,j], Y[i,j], w) for w in panels_img])+Qwell_1*InflueceU_W(X[i,j], Y[i,j], wells[0])
v[i,j] =sum([w.Q*InflueceV(X[i,j], Y[i,j], w) for w in panels_img])+Qwell_1*InflueceV_W(X[i,j], Y[i,j], wells[0])
In [509]:
# plots the streamlines
%matplotlib inline
size = 6
pyplot.figure(figsize=(size, size))
pyplot.grid(True)
pyplot.title('Streamline Plot')
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');
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 [510]:
size = 7
pyplot.figure(figsize=(size, size-1))
pyplot.title('Pressure field (psi)')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.contour(X, Y, p, 15,linewidths=0.5, colors='k')
pyplot.contourf(X, Y, p, 15,cmap='rainbow',)
pyplot.colorbar() # draw colorbar
Out[510]:
In [511]:
size = 7
pyplot.figure(figsize=(size, size-1.5))
pyplot.title('Total Velocity field (ft/day)')
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
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=110, vmin=0)
pyplot.colorbar() # draw colorbar
Out[511]:
In [512]:
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[512]:
In [513]:
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[513]:
In [ ]:
In [ ]:
In [ ]: