diffusion.py
#
# Unsteady diffusion in cooling blades
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['toolbar'] = 'None'
myColorMap = matplotlib.cm.jet
# =========================================================================
def diffusionSolve(betaRequested,n):
#
# -1- Geometrical and Material Data
#
rho = 2707; c = 896; k = 204
alpha = k / (rho * c)
Ly = 0.015; Lx = Ly/6; Lt = 60
qin = 30000; Tf = 300; hx = 60; hy = 100
print('Geometrical and Physical Data :')
print('==========================================================')
print(' rho = %10.3e [W/(m2 K)] (density) ' % rho)
print(' c = %10.3e [J/(kg K)] (heat capacity) ' % c)
print(' k = %10.3e [W/(m K)] (thermal conductivity) ' % k)
print(' alpha = %10.3e [m2/s] (thermal diffusivity) ' % alpha)
print(' qin = %10.3e [W/(m2)] (heat flux density) ' % qin)
print(' hx = %10.3e [W/(m2 K)] (convection coefficient) ' %hx)
print(' hy = %10.3e [W/(m2 K)] (convection coefficient) ' % hy)
print(' Tf = %10.3e [K] ' % Tf)
print(' Lx = %10.3e [m] ' % Lx)
print(' Ly = %10.3e [m] ' % Ly)
print(' Lt = %10.3e [s] \n' % Lt)
#
# -2- Discretization Parameters
#
dx = Lx/n; nx = n + 1; ny = 6*n + 1
dt = (dx * dx) * betaRequested / alpha
nt = int(ceil(Lt/dt)); dt = Lt/nt
beta = alpha * dt / (dx * dx)
print('Discretization :')
print('==========================================================')
print(' dx = %10.3e [m] (nx = %i, ny = %i)' % (dx,nx,ny))
print(' dt = %10.3e [m] (nt = %i)' % (dt,nt))
print(' => beta requested = %14.7e [m] ' % betaRequested)
print(' => beta selected = %14.7e [m] \n' % beta)
#
# -3- Computation
#
T = Tf * ones((nx+2,ny+2))
Tcurve = Tf * ones((3,nt+1))
iPlot = set((array([0.2,1,30,60])/dt).astype(int))
for it in range(nt):
T[1:-1,0] = T[1:-1,2] - 2*dx*hy*(T[1:-1,1]-Tf) /k
T[1:-1,-1] = T[1:-1,-3] + 2*dx*qin /k
T[0,1:-1] = T[2,1:-1]
T[-1,1:-1] = T[-3,1:-1] - 2*dx*hx*(T[-2,1:-1]-Tf) /k
T[1:-1,1:-1] += beta*(T[2:,1:-1] + T[:-2,1:-1]
+ T[1:-1,2:] + T[1:-1,:-2] - 4*T[1:-1,1:-1])
Tcurve[0,it+1] = T[1,1]
Tcurve[1,it+1] = T[1,(ny+1)//2]
Tcurve[2,it+1] = T[1,-2]
if ((it+1) in iPlot):
#
# -4- Graphics and results
#
plt.figure("Cooling blade temperature at t = %4.1f [s]" % (it*dt))
X,Y = meshgrid(linspace(0,Lx,n+1),linspace(0,Ly,6*n+1))
Tplot = T[1:-1,1:-1].T - Tf
Trange = arange(floor(amin(Tplot)*10)/10,ceil(amax(Tplot)*10)/10+0.1,0.1)
plt.subplot(1,2,2)
plt.contourf( X,Y,Tplot,Trange,cmap=myColorMap)
plt.contourf(-X,Y,Tplot,Trange,cmap=myColorMap)
plt.colorbar(shrink=0.7)
plt.contour( X,Y,Tplot,Trange,colors='k',linewidths=1)
plt.contour(-X,Y,Tplot,Trange,colors='k',linewidths=1)
plt.plot([-Lx,Lx,Lx,-Lx,-Lx],[0,0,Ly,Ly,0],'-k',linewidth=1)
plt.plot([0,0,0],[0,Ly/2,Ly],'ok')
plt.plot([0,0],[-0.0005,Ly+0.0005],'ow')
plt.axis("equal");
plt.axis("off")
plt.subplot(1,2,1)
t = arange(0,it*dt+dt,dt)
plt.plot(t,Tcurve[0,0:len(t)]-Tf,'-b',
t,Tcurve[1,0:len(t)]-Tf,'-g',
t,Tcurve[2,0:len(t)]-Tf,'-r')
# plt.savefig("diffusion-time%d.pdf" % int(it*dt+dt))
print('Results at t = %4.1f [s] :' % (it*dt))
print('==========================================================')
print(' T(0,0 )-Tf = %14.6e [K] ' % (T[1,1]-Tf))
print(' T(0,Ly/2)-Tf = %14.6e [K] ' % (T[1,(ny+1)//2]-Tf))
print(' T(0,Ly )-Tf = %14.6e [K] \n' % (T[1,-2]-Tf))
plt.show()
# ============================= mainProgram ===============================
diffusionSolve(0.25,5)
#diffusionSolve(0.125,5)
# =========================================================================