Example: Thermal Simulation¶
Objective¶
We will simulate the temporal change in temperature of the copper rectangle whose length is 50 cm (= LX), width is 30 cm (= LY) using a finite-difference method.
Physics¶
Governing Equation¶
\[\frac{\partial}{\partial t} T=\frac{k}{\rho c}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)T\]
Variable |
Description |
---|---|
\(T\) |
Temperature [K] |
\(t\) |
Time [s] |
\(x\) |
x coordinate [m] |
\(y\) |
y coordinate [m] |
Constant |
Description |
Value of Copper |
---|---|---|
\(k\) |
Thermal conductivity |
398.0 [W/m K] |
\(\rho\) |
Density |
8960.0 [kg/m 3] |
\(c\) |
Specific heat capacity |
385.0 [J/kg K] |
Boundary Conditions¶
\[\begin{split}\begin{eqnarray*}
T &=& T_1 ~ (x=0, x={\rm LX}, y=0) \\
T &=& T_1+T_2 \sin \frac{\pi x}{\mathrm LX} ~ (y={\rm LY}) \\
\frac{\partial}{\partial x}T &=& 0 ~ (x=0, x={\rm LX}) \\
\frac{\partial}{\partial y}T &=& 0 ~ (y=0, y={\rm LY}) \\
T_1 &=& 40 ~ [^\circ C] \\
T_2 &=& 60 ~ [^\circ C] \\
\end{eqnarray*}\end{split}\]
Initial Conditions¶
\[\begin{split}\begin{eqnarray*}
T &=& T_0 \\
T_0 &=& 20 ~ [^\circ C]
\end{eqnarray*}\end{split}\]
Program¶
import nlcpy as vp
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
from matplotlib.animation import PillowWriter
LX = 50e-2
LY = 30e-2
T0 = 20.0
T1 = 40.0
T2 = 60.0
HC = 398.0 / (8960.0 * 385.0)
def initialize(grid):
grid.fill(T0)
grid[:, 0] = T1
grid[:, -1] = T1
grid[0] = T1
grid[-1] = T1 + T2 * \
vp.sin(vp.pi * vp.linspace(0, LX, grid.shape[1]) / LX)
def compute(grid, coef, temp, mx, my):
vp.multiply(grid[2:my, 1:mx-1], coef[0], out=temp)
temp += grid[1:my-1, 2:mx] * coef[1]
temp += grid[1:my-1, 1:mx-1] * coef[2]
temp += grid[0:my-2, 1:mx-1] * coef[0]
temp += grid[1:my-1, 0:mx-2] * coef[1]
grid[1:my-1, 1:mx-1] = temp[...]
def plot(grid, ax, ims):
x = vp.linspace(0, LX, grid.shape[1])
y = vp.linspace(0, LY, grid.shape[0])
try:
# meshgrid is not implemented yet.
xx, yy = vp.meshgrid(x, y)
except AttributeError:
import numpy as np
xx, yy = np.meshgrid(x, y)
im = ax.plot_wireframe(xx, yy, grid, rstride=10, cstride=10)
ims.append([im])
def create_animation(fig, ims):
ani = animation.ArtistAnimation(fig, ims, interval=200)
ani.save("thermal_simulation.gif", writer='pillow')
def heatequation(
nx, # The number of grid points in X-direction.
ny, # The number of grid points in Y-direction.
dt, # The time step interval.
mt, # The maximum number of time steps.
kp, # The number of time steps for drawing interval.
):
mx = nx + 2
my = ny + 2
grid = vp.empty((my, mx), dtype='f8')
temp = vp.empty((ny, nx), dtype='f8')
dx = LX / (nx + 1)
dy = LY / (ny + 1)
coef = [
(HC * dt) / (dx * dx),
(HC * dt) / (dy * dy),
1.0 - HC * dt * (2.0 / (dx * dx) + 2.0 / (dy * dy)),
]
print("initializing grid...", end="", flush=True)
initialize(grid)
print("done", flush=True)
ims = []
fig = plt.figure()
ax = Axes3D(fig, rect=(0., 0., .95, 1.))
ax.set_xlabel("x[m]")
ax.set_ylabel("y[m]")
ax.set_zlabel("T[$^{\circ}$C]")
ax.zaxis.set_rotate_label(False)
print("computing difference method...", end="", flush=True)
for i in range(int(mt/dt)):
compute(grid, coef, temp, mx, my)
if i % int(kp/dt) == 0:
plot(grid, ax, ims)
print("done", flush=True)
print("creating animation...", end="", flush=True)
create_animation(fig, ims)
print("done", flush=True)
if __name__ == "__main__":
heatequation(500, 300, 0.001, 30, 1.)