3D-heat-conduction.py

from fipy import CellVariable, FaceVariable, TransientTerm, DiffusionTerm, Viewer, Grid3D
from fipy.meshes.nonUniformGrid3D import NonUniformGrid3D
import numpy as np
import matplotlib.pyplot as plt

FT2M = 0.3048

uniformMesh = True
solveLinear = False
noViewer = False
plotTdependence = False

if uniformMesh:
    domainX = 30.0
    domainY = 27.5
    domainZ = 70.0

    nx = 33*2
    ny = 28*2
    nz = 70*2

    print "uniform mesh"
    print "number elements",nx,'*',ny,'*',nz,'=',nx*ny*nz
    m = Grid3D(nx=nx, ny=ny, nz=nz, Lx=domainX*FT2M, Ly=domainY*FT2M, Lz=domainZ*FT2M)

else:
    # grid spacing
    dz = (2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,
          1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 0.5,0.5,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5, 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
          2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0) 
    
    dx = (0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          1.0,1.0,1.0,1.0,1.0)  # symmetric about middle heater
    
    dx = dx[::-1]
    
    dy = (1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
          1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
          0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,
          0.5,0.5,0.5,0.5,0.25,0.25)  # symmetric about midline of heaters
    
    domainX = sum(dx)
    domainY = sum(dy)
    domainZ = sum(dz)
    
    print "non-uniform mesh"
    print "number elements",len(dx),'*',len(dy),'*',len(dz),'=',len(dx)*len(dy)*len(dz)
    
    # convert from feet (specified above) to meters (used in calculation)
    f = lambda x:x*FT2M
    dx = map(f,dx)
    dy = map(f,dy)
    dz = map(f,dz)
    
    m = NonUniformGrid3D(dx=dx, dy=dy, dz=dz)

print "domain dimensions (ft)",domainX,domainY,domainZ

x,y,z = m.cellCenters.value

# convert cell center coordinates (not used in calculation) 
# to feet (for specifying geometry) 
x /= FT2M
y /= FT2M
z /= FT2M

# flip domain in x for simplified default plotting in MayaVi2
x = domainX-x

# (10 ft high drift centered in 70 ft high domain)
driftFloor  = 30.0 # relative z coordinate from bottom (ft) 

driftHeight = 10.0 # ft (z-direction)
driftWidth =  15.0 # ft (y-direction)
# drift runs entire x-length of domain

driftRib = domainY - driftWidth/2.0 # relative y coordinate from far side (ft)

drift   =  (z>driftFloor) & (z<driftFloor+driftHeight) & (y>driftRib) 

spslope = 6.0/7.0  # 36 degree angle of repose
spintercept = 6 + spslope*9
saltPile = ((z>driftFloor) & (z<driftFloor+6.0) & 
            (z<driftFloor+spintercept-spslope*x) & (y>driftRib))

# unit circles centered at x={0,3,6}, z=31
# 2.5 ft of space on either end of heaters
htr1 = (np.sqrt(      x**2 + (z-(driftFloor+1))**2) < 1.0) & (y>driftRib+2.5) 
htr2 = (np.sqrt((x-3.0)**2 + (z-(driftFloor+1))**2) < 1.0) & (y>driftRib+2.5) 
htr3 = (np.sqrt((x-6.0)**2 + (z-(driftFloor+1))**2) < 1.0) & (y>driftRib+2.5) 
heaters = htr1 | htr2 | htr3

# top of marker bed 139 is 6.3 ft below bottom of drift (2.8 ft thick)
MB139 = (z>driftFloor-6.3-2.8) & (z<driftFloor-6.3)

# initial temp 300 K
T0 = 26.85

if solveLinear:
    T = CellVariable(name='linear T', mesh=m, value=T0)
else:
    T = CellVariable(name='non-linear T', mesh=m, value=T0, hasOld=1)

### for quickly debugging geometry issues (contour plots of T)
##T.setValue(50,  where=MB139)
##T.setValue(100, where=drift)
##T.setValue(75,  where=saltPile)
##T.setValue(125, where=heaters)

Q = CellVariable(name='source', mesh=m, value=0.0)

# heater canisters are cylinders 2 ft in diameter, 9 ft long.
# ~ 0.8 cubic meters
vol = np.pi*9.0* (FT2M**3)

heaterQ = 1500.0/vol  # volume source has units W/m^3 

print 'volumetric heater power:',heaterQ

# heaters 
Q.setValue(heaterQ, where=heaters)

# parameter value sources
#-------------------------------
# http://en.wikipedia.org/wiki/List_of_thermal_conductivities
# BAMBUS II report p. 45 (Figure 2.37)
# Munson (8/22/1989) "Proposed New Structural Reference Stratigraphy, Law, and Properties" 
# Phil ppt presentation 7/10/2013

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# density [gram/meter^3]

# Munson (1989)
intactSaltDensity = 2.163E6 
anhydriteDensity  = 2.963E6
polyhaliteDensity = 2.780E6

crushedSaltDensity = 1.50E6
sandDensity = 2.33E6
airDensity = 1.2E+3

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# specific heat capacity [Joules/(gram*K)]
# volumetric heat capacity [Joules/(meter^3*K)]
# Watt = Joule/second => t is seconds (since RHS of PDE is Watts)

# Munson (1989)
intactSaltRhoCp = intactSaltDensity*0.8628  
anhydriteRhoCp  =  anhydriteDensity*0.7333 
polyhaliteRhoCp = polyhaliteDensity*0.890 

crushedSaltRhoCp = crushedSaltDensity*0.8628  # intact salt?
sandRhoCp        = sandDensity*0.8
airRhoCp         = airDensity*1.012

print 'rho * Cp (J/(m^3 * K))'
print 'intact salt',intactSaltRhoCp
print 'crushed salt',crushedSaltRhoCp
print 'air',airRhoCp
print 'glass/sand',sandRhoCp

# volumetric heat capacity [J/(m^3 * K)]
rhocp = CellVariable(mesh=m, value=intactSaltRhoCp) 

rhocp.setValue(airRhoCp,         where=drift)
rhocp.setValue(crushedSaltRhoCp, where=saltPile)
rhocp.setValue(sandRhoCp,        where=heaters)
rhocp.setValue(anhydriteRhoCp,   where=MB139)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# thermal conductivity [Watts/(meter * Kelvin)] = [Joules/(meter * sec * Kelvin)]

# Munson (1989)
intactSaltKap = 5.4 
anhydriteKap  = 4.7
polyhaliteKap = 1.4

crushedSaltKap = 1.0 # 30% porosity (Bambus II)
sandKap = 1.1 # borosilicate glass (from Phil 7/10/13)
airKap = 14.0 # effective thermal conductivity 
# including approximate radiation and convection (Phil 7/10/13)

# airKap = 0.0262  #(1 atm, 300K)   <- conduction only

if not solveLinear:
    
    saltKexp = 1.14
    anhydriteKexp = 1.15
    polyhaliteKexp = 0.35

    # thermal conductivity exponent
    kexp = CellVariable(mesh=m, value=saltKexp)

    # air and heaters don't have nonlinear kappa 
    kexp.setValue(0.0,           where=drift)  
    kexp.setValue(saltKexp,      where=saltPile)
    kexp.setValue(0.0,           where=heaters)
    kexp.setValue(anhydriteKexp, where=MB139)

# thermal conductivity at 300 K [W/(m deg-C)]
kappa = CellVariable(mesh=m, value=intactSaltKap) 

kappa.setValue(airKap,         where=drift)
kappa.setValue(crushedSaltKap, where=saltPile)
kappa.setValue(sandKap,        where=heaters)
kappa.setValue(anhydriteKap,   where=MB139)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# heat conduction governing equation
if solveLinear:
    eq = TransientTerm(rhocp, var=T) == DiffusionTerm(kappa, var=T) + Q
else:
    # T dependence of thermal conductivity from Munson (1989) 
    eq = TransientTerm(rhocp, var=T) == DiffusionTerm(kappa*(300/(T+273.15))**kexp, var=T) + Q

    if plotTdependence:
        tp = np.linspace(10,300)
        plt.plot(tp,intactSaltKap*(300/(tp+273.15))**saltKexp,      'r-',label='salt')
        plt.plot(tp, anhydriteKap*(300/(tp+273.15))**anhydriteKexp, 'g-',label='anhydrite')
        plt.plot(tp,polyhaliteKap*(300/(tp+273.15))**polyhaliteKexp,'b-',label='polyhalite')
        plt.xlabel('$T$ (degrees C)')
        plt.ylabel('$\\kappa$ (W/(m $\cdot$ K))')
        plt.legend(loc=0)
        plt.grid()
        plt.savefig('temp-dependence-kappa.png')

if noViewer:
    #fh = open('output.dat','w')
    pass
else:
    viewer = Viewer(vars=(T),datamin=T0,datamax=200)
    viewer.plot()

t = 0.0
t1 = 86400.0 # 1 day

# linearly increasing time steps ~100 days total
# constant timesteps 2 years, 2-day timesteps

for j,step in enumerate(np.concatenate((np.linspace(1,2*t1,100),np.ones((365,))*2*t1),axis=0)):
    t += step

    if solveLinear:
        eq.solve(dt=step) # linear timestep

    else:
        T.updateOld() # advance timestep

        res = 1000.0
        counter = 0
        while res > 1.0E-4 and counter < 10:
            res = eq.sweep(dt=step) # iterate nonlinear timestep
            counter += 1
        
        print 'DEBUG %3i %.3g %i' % (j,res,counter)

    if j%10 == 0:
        print "step %i: %.4g days" % (j,t/t1)
        if noViewer:
            if j%100 == 0:
                pass
                #fh.write('\n# t=%.5g\n' % t)
                #np.savetxt(fh,T,fmt='%.3f')
        else:
            viewer.plot()

#fh.close()