Python program to calculate fluid temperature in borehole

borehole_temperature.py

"""
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()