Skip to content
Snippets Groups Projects
Commit 13201b19 authored by Miro Pütz's avatar Miro Pütz
Browse files

Replace 2d_heat_conduction_v2.py

parent 012fe9bf
Branches master
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
2d_heat_conduction_v2.py
+ 76
− 30
View file @ 13201b19
......@@ -2,97 +2,143 @@
# for a non-moving medium with constant conductivity
import numpy as np
import matplotlib.pyplot as plt
import strength_calc
plt.close('all')
# -------------------- Numerical model parameters -----------------------------#
# Boundary conditions = constant temperature (gradient)
# Model size [m]
X_SIZE = 100000 # Profile length [m]
Y_SIZE = 150000 # Depth [m]
X_SIZE = 200e3 # Profile length [m]
Y_SIZE = 150e3 # Depth [m]
# Number of nodes
X_NUM = 201
Y_NUM = 201
X_NUM = 301
Y_NUM = 301
# Grid step
X_STEP = float(X_SIZE)/(X_NUM-1) # [m]
Y_STEP = float(Y_SIZE)/(Y_NUM-1) # [m]
Y_STEP2 = Y_STEP**2 # to save cpu
X_STEP2 = X_STEP**2
TIMESTEPS = 1000 # Number of timesteps
TIMESTEPS = 5000 # Number of timesteps
# Material properties
K_uc = 2.75 # Thermal conductivity: upper crust [W/(m*K)]
K_m = 2.75 # mantle
CP = 1000 # Heat cpacity [J/(kg*K)]
RHO = 3200 # Density [kg/m**3]
RHOCP = RHO*CP # Volumetric heat capacity [J/m**3]
KAPPA = float(K_uc)/RHOCP # Thermal diffusivity [m**2/s]
RHO_uc = 2750 # Density [kg/m**3]
RHO_m = 3300
RHOCP_uc = RHO_uc*CP # Volumetric heat capacity [J/m**3]
RHOCP_m = RHO_m*CP
KAPPA_uc = float(K_uc)/RHOCP_uc # Thermal diffusivity [m**2/s]
KAPPA_m = float(K_m)/RHOCP_m
S_uc = 0.2/1.e6 # Heat Production: upper crust [W/m**3]
S_m = 0.05/1.e6 # mantle
S_m = 0.1/1.e6 # mantle
QMAX = 20./1000 # Basal heat flow [W/m**2]
ZC = 35000. # Crustal thickness [m]
ZL = 125000. # Lithospheric thickness [m]
ZC = 45e3 # Crustal thickness [m]
ZL = 125e3 # Lithospheric thickness [m]
# Vectors for nodal points positions
cols = np.linspace(0, X_SIZE, X_NUM)
rows = np.linspace(0, Y_SIZE, Y_NUM)
# Timestep
dtmax = min(X_STEP, Y_STEP)**2/(3*KAPPA) # timestep limitation for stability
dt = 0.9*dtmax # timestep [s]
dtmax = min(X_STEP, Y_STEP)**2/(3*max(KAPPA_uc, KAPPA_m)) # timestep limitation for stability
dt = 0.8*dtmax # timestep [s]
time = (TIMESTEPS*dt)/(10**6*365.25*24*3600) # time in Myr
# ----------------- initial conditions (temperature profile) ------------------#
# create temperature arrays
index_zc = round(ZC/Y_STEP)+1 # index for crust, mantle boundary
temp_gradient = np.zeros((Y_NUM, X_NUM))
temp_gradient[:47] = (-(S_uc * rows[:47]**2.)/(2.*K_uc) + (QMAX + S_uc * ZL)/K_uc * rows[:47])[:,np.newaxis]
temp_gradient[47:] = (-(S_m * rows[47:]**2.)/(2.*K_m) + (QMAX + S_m * ZL)/K_m * rows[47:])[:,np.newaxis]
temp_gradient[:] = (-(S_m * rows[:]**2.)/(2.*K_m) + (QMAX + S_m * ZL)/K_m * rows[:])[:,np.newaxis]
#temp_gradient[:index_zc] = (-(S_uc * rows[:index_zc]**2.)/(2.*K_uc) + (QMAX + S_uc * ZL)/K_uc * rows[:index_zc])[:,np.newaxis]
#temp_gradient[index_zc:] = (-(S_m * rows[index_zc:]**2.)/(2.*K_m) + (QMAX + S_m * ZL)/K_m * rows[index_zc:])[:,np.newaxis]
told = temp_gradient*np.ones((len(temp_gradient), X_NUM))
tnew = told # new temperature
# create logical array for magma chambers
magma_chamber = 900 # temperature in grad
chambers = [{'x': 140, 'y': 30, 'width': 35, 'hight': 15},
{'x': 50, 'y':20, 'width': 15, 'hight': 5},]
# create logical array for magma chambers
mask = np.zeros((Y_NUM, X_NUM), dtype=bool)
for i in range(X_NUM):
for j in range(Y_NUM):
if ( ( (i*X_STEP-70000)**2+(j*Y_STEP-40000)**2 <= 13000**2)
or ((i*X_STEP-25000)**2+(j*Y_STEP-50000)**2 <= 10000**2) ):
mask[j,i] = True # True if there is anomaly (magma chamber)
for c in chambers:
if ( (i*X_STEP-c['x']*1000)**2/(c['width']*500)**2+(j*Y_STEP-c['y']*1000)**2/(c['hight']*500)**2 <= 1):
mask[j,i] = True # True if there is anomaly (magma chamber)
told[mask] = magma_chamber
tplot_initial = np.copy(told) # for plotting the initial temperature profile
t_initial = np.copy(told) # for plotting the initial temperature profile
# ----------------------------- Explicit solving ------------------------------#
for i in xrange(TIMESTEPS):
tnew[1:-1, 1:-1] = dt*KAPPA*( (told[2:, 1:-1]-2*told[1:-1, 1:-1]+told[:-2, 1:-1])/Y_STEP2 +
(told[1:-1, 2:]-2*told[1:-1, 1:-1]+told[1:-1,:-2])/X_STEP2)+told[1:-1, 1:-1]+S_uc
#tnew[1:-1, 1:-1] = dt*KAPPA*( (told[2:, 1:-1]-2*told[1:-1, 1:-1]+told[:-2, 1:-1])/Y_STEP2 +
# (told[1:-1, 2:]-2*told[1:-1, 1:-1]+told[1:-1,:-2])/X_STEP2)+told[1:-1, 1:-1]+S_uc
# for crust
tnew[1:index_zc, 1:-1] = dt*KAPPA_uc*( (told[2:index_zc+1, 1:-1]-2*told[1:index_zc, 1:-1]+told[:index_zc-1, 1:-1])/Y_STEP2 +
(told[1:index_zc, 2:]-2*told[1:index_zc, 1:-1]+told[1:index_zc,:-2])/X_STEP2)+told[1:index_zc, 1:-1]+S_uc
#for mantle
tnew[index_zc:-1, 1:-1] = dt*KAPPA_m*( (told[index_zc+1:, 1:-1]-2*told[index_zc:-1, 1:-1]+told[index_zc-1:-2, 1:-1])/Y_STEP2 +
(told[index_zc:-1, 2:]-2*told[index_zc:-1, 1:-1]+told[index_zc:-1,:-2])/X_STEP2)+told[index_zc:-1, 1:-1]+S_m
told = np.copy(tnew)
if (i*dt)/(10**6*365.25*24*3600) < 2.5: # time after magma chamber is gone
if (i*dt)/(10**6*365.25*24*3600) < 4.: # time after magma chamber is gone
told[mask] = magma_chamber
np.save('temp_array', told) # save the final temp. array in a numpy format
#np.save('temp_array', told) # save the final temp. array in a numpy format
#np.save('initial_array', tplot_initial)
# for the strength calculation
# ------------------------------- Plots ---------------------------------------#
# --------------------------- strength calculation ----------------------------#
strength_vector_initial = np.zeros(X_NUM)
strength_vector_final = np.zeros(X_NUM)
total_strength_initial = 0
total_strength_final = 0
for i in range(X_NUM):
strength_vector_initial[i] = strength_calc.strength(t_initial[:,i], rows, ZC)
strength_vector_final[i] = strength_calc.strength(told[:,i], rows, ZC)
total_strength_initial += strength_calc.strength(t_initial[:,i], rows, ZC)
total_strength_final += strength_calc.strength(told[:,i], rows, ZC)
percent = 100-(total_strength_final*100)/total_strength_initial
print 'After %.3f Myr the whole profile is %.3f percent weaker' % (time, percent)
# ------------------------------- Plots ---------------------------------------#
# Initital conditions
plt.figure(1)
X, Y = np.meshgrid(cols,rows)
im = plt.contourf(X/1000, Y/1000, tplot_initial, 100, cmap=plt.cm.hot, interpolation='nearest', origin='lower')
plt.contourf(X/1000, Y/1000, t_initial, 100, cmap=plt.cm.hot, interpolation='nearest', origin='lower')
plt.title('Initial Temperature Profile')
plt.xlabel('Profile Distance [km]')
plt.ylabel('Depth [km]')
plt.gca().invert_yaxis()
plt.colorbar(plt.clim(0,np.max(told)))
plt.show(block=False)
plt.show()
# Final plot
plt.figure(2)
myplot = plt.contourf(X/1000., Y/1000., told, 100, cmap=plt.cm.hot, interpolation='nearest', origin='lower')
plt.contourf(X/1000., Y/1000., told, 100, cmap=plt.cm.hot, interpolation='nearest', origin='lower')
#plt.plot([cols[200]/1000, cols[200]/1000], [0,150], 'w')
plt.title('After %.3f Myr' %time)
plt.xlabel('Profile Distance [km]')
plt.ylabel('Depth [km]')
plt.gca().invert_yaxis()
plt.colorbar(plt.clim(0,np.max(told)))
plt.show(block=False)
\ No newline at end of file
plt.show()
# strength plot
plt.figure(3)
plt.plot(cols/1000, strength_vector_initial/1000, label='Initial Strength Profile')
plt.plot(cols/1000, strength_vector_final/1000, label='After %.3f Myr' %time)
plt.title('Differential Stress Profile')
plt.xlabel('Profile Distance [km]')
plt.ylabel('Differential stress [GPa]')
plt.legend(loc=4)
plt.show()
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment