Skip to content

Instantly share code, notes, and snippets.

@jkfurtney
Created May 19, 2022 13:38
Show Gist options
  • Save jkfurtney/e21b8e33ad943468b78a81e1fa328e71 to your computer and use it in GitHub Desktop.
Save jkfurtney/e21b8e33ad943468b78a81e1fa328e71 to your computer and use it in GitHub Desktop.
Python program to calculate fluid temperature in borehole
"""
This file is the same as borehole.py but has an non-uniform cartesian mesh
Coupled thermal flow model for brine temperature and wall
temperature. Rock is an axisymmetric 2D mesh
A two-dimensional axi-symmetric domain represents the heat in the rock
surrounding the borehole. The rock is initially at the geo-thermal
gradient and is heated as hot brine is injected. The rock is heated as
the brine cools.
\[
\frac{\partial T_s}{\partial t} = -\frac{Q}{\pi r^2}\frac{\partial T_s}{\partial y} - \frac{h}{\rho C_p}\frac{2 \pi r}{\pi r^2}\left(T_s-T_w\right)
\]
\[
\frac{\partial T_w}{\partial t} = -\frac{\alpha}{x}\nabla\cdot\left( x \nabla T_w \right)
\]
\[
\left.\frac{\partial T_w}{\partial x} \right|_{x=r}=-\frac{h}{k}\left(T_s-T_w\right)
\]
\[T_s(y,0)=T0+\textrm{geotherm}\ y\]
\[T_w(x,y,0)=T0+\textrm{geotherm}\ y\]
\[T_s(0,t)=T_\infty\]
where $T_s$ is brine temperature, $T_w$ is the wall rock temperature
(axisymmetric about y=0), $t$ is time, $r$ is the borehole radius, $Q$
is the volumetric flow-rate and $y$ is the length along the borehole
and x is the radial direction.
\begin{verbatim}
+--> x
|
|
y
+--+--+--+--+--+
|s |w |w |w |w |
+--+--+--+--+--+
|s |w |w |w |w |
+--+--+--+--+--+
|s |w |w |w |w |
+--+--+--+--+--+
^ ^
| |
0 r0
\end{verbatim}
"""
import fipy as fp
import pylab as pl
from math import pi
from scipy.constants import foot, inch, day, week, hour, gallon, minute, day
from scipy.constants import F2C,C2F
import numpy as np
r0 = 0.08
rock_k = 2.8
rock_Cp = 740.0
rhom = 2200.0
alpha_rock = rock_k/rock_Cp/rhom
h = 978.0
rhow = 958.0
water_Cp = 4181.3
print "alpha_rock {}".format(alpha_rock)
ly = 1000.0
ny = 100
Q = 100.0/day
u = Q/pi/r0**2
Tf = 100.0
T0 = 20.0
geotherm = 0.03 # deg C / m
Tbed = T0 + geotherm*ly
print "Tbed {}".format(Tbed)
nx = 25
dx = np.logspace(np.log10(2.5e-3),np.log10(2.5),nx)
lx = dx.sum()
dy = ly/ny
time_flush = ly/u
print "flush time", time_flush/60., "minutes"
A = 2*pi*r0*dy # ymesh element surface area
V = pi*r0**2*dy # ymesh element volume
xymesh = fp.Grid2D(nx=nx, dx=dx, ny=ny, dy=dy)
ymesh = fp.Grid1D(nx=ny, dx=dy)
# for cylindrical geometry coefficients
rc = fp.CellVariable(mesh=xymesh, value=(xymesh.getCellCenters()[0]+r0))
faces = np.array([[r0+n*dx for n in range(nx)] for nny in range(ny)])
from itertools import chain
faces = np.array(list(chain(*faces))) # flatten list hack
cellVolumes = pi*((faces+dx)**2-faces**2)*dy
cellVolumes = xymesh.getCellVolumes()*2*pi*rc.value
vol = ly*((r0+lx)**2*pi-pi*r0**2)
assert abs(vol - cellVolumes.sum())/max(abs(vol),
abs(cellVolumes.sum())) < 2e-2
Tw = fp.CellVariable(mesh=xymesh, value=T0, hasOld=0, name='Tw')
Ts = fp.CellVariable(mesh=ymesh, value=T0, hasOld=0, name='Ts')
# for visualization
X = xymesh.getCellCenters().value[0].reshape(ny,nx)
Y = xymesh.getCellCenters().value[1].reshape(ny,nx)
Twa = Tw.value.reshape(ny,nx)
#apply geothermal gradient to initial conditions
Tw.setValue(xymesh.getCellCenters().value[1]*geotherm+T0)
Ts.setValue(ymesh.getCellCenters().value[0]*geotherm+T0)
alpha = fp.CellVariable(mesh=xymesh)
alpha.setValue(alpha_rock)
#mask for borehole wall surface
mask = fp.CellVariable(xymesh, value=0)
mask.setValue(1, where=xymesh.getCellCenters()[0]<dx[0])
# solid BC
# noflux on top and bottom, fixed T on outer radius flux term
# is volumetric
BC1 = ( fp.FixedValue(faces=xymesh.getFacesRight(), value = Ts.value))
# fluid BC # incoming fluid is constant
BC2 = ( fp.FixedValue(faces=ymesh.getFacesLeft(), value=Tf),
fp.FixedValue(faces=ymesh.getFacesRight(), value=T0))
# energy balance
e_0_rock = (Tw.value*cellVolumes*rhom*rock_Cp).sum()
e_0_fluid = (Ts.value*rhow*water_Cp*V).sum()
e_in = 0
e_out = 0
dt = 2.5
time = 0.0
#tend = 4946.4 * hour
tend = 4.1*day
#tend = 3.5 * hour
fluid_dt = 0.75 * dy/abs(u)
dt = fluid_dt
#assert dt < fluid_dt
print "fluid_dt", fluid_dt
#an array to apply the fluid temperature to the rock
Tboundxy = fp.CellVariable(mesh=xymesh, value=-1e100, name="Tboundxy")
#an array to apply the wall temperature to the fluid
Tboundy = fp.CellVariable(mesh=ymesh, value=-1e100, name="Tboundy")
# source term from fluid to rock
c_in = A*h*(Tboundxy-Tw)/rock_Cp/rhom/cellVolumes[0]
# sink term from fluid to rock
c_out = -A*h*(Ts-Tboundy)/water_Cp/rhow/V
# axisymmetric heat domain (solid rock)
eq1 = fp.TransientTerm(coeff=rc) == fp.DiffusionTerm(coeff=alpha*rc) \
+ mask*c_in*rc
#1d fluid flow domain
eq2 = fp.TransientTerm() == fp.VanLeerConvectionTerm(coeff=(-u,)) + c_out
# plot recording interval
i=0
ic=0
yy=ymesh.getCellCenters().value[0]
timel=[0]
templ=[Ts.value[-1]]
dt = 10.0
interval = int(tend/dt/5.0)
while time<=tend:
#for jj in range(2):
# copy fluid temperature to rock mesh
Tboundxy.value.reshape(ny, nx)[:, 0] = Ts.value
#copy rock wall temperature to fluid mesh
Tboundy.value = Tw.value.reshape(ny,nx)[:,0]
eq1.solve(var=Tw, boundaryConditions=BC1, dt=dt)
eq2.solve(var=Ts, boundaryConditions=BC2, dt=dt)
e_in += Q * dt * Tf * rhow * water_Cp
e_out += Q * dt * Ts.value[-1] * rhow * water_Cp
e_rock = (Tw.value * cellVolumes * rhom * rock_Cp).sum()
e_fluid = (Ts.value * rhow * water_Cp * V).sum()
if i % interval == 0:
print "output"
timel.append(time)
templ.append(Ts.value[-1])
pl.subplot2grid((2, 2), (0, 0), rowspan=2)
pl.contourf(np.log10(X/foot), Y[::-1] / foot,
Twa,np.linspace(T0 - 1, Tf, 10))
pl.ylabel("Depth (feet)")
pl.xlabel("Distance from borehole Log(feet)")
pl.subplots_adjust(wspace=0.3)
pl.subplot2grid((2, 2), (0, 1))
pl.plot(yy / foot,C2F(Ts.value))
pl.plot(yy / foot,C2F(Tboundy.value))
pl.ylabel("Temperature (F)")
pl.ylim(C2F((T0, Tf)))
pl.subplot2grid((2, 2), (1, 1))
pl.plot(np.array(timel) / hour, C2F(templ))
pl.ylabel("Exit Temperature (F)")
pl.xlabel("Time (hours)")
pl.ylim(C2F((T0, Tf)))
pl.savefig("aname%05i.png" % (ic))
ic += 1
i += 1
time += dt
e_0 = e_0_fluid + e_0_rock
deltaE = e_in - e_out
de_fluid = e_fluid - e_0_fluid
de_rock = e_rock - e_0_rock
print (e_0 + e_in - e_out -(e_rock + e_fluid)) / (e_fluid + e_rock)
pl.cla()
sand_x, sand_z = np.loadtxt("./data1.csv", skiprows=6, delimiter=",").T
pl.plot(Ts.value[::-1],yy[::-1])
pl.plot(sand_x, sand_z, "o", linewidth=3)
pl.xlabel("Temperature [$^o$C]")
pl.ylabel("Depth [m]")
pl.plot([20,50],[0,1000])
pl.legend(("Numerical Model","Exact Solution","Geotherm"), loc=4)
pl.gca().invert_yaxis()
pl.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment