# /usr/bin/env python
import copy
from math import *
import numpy as np
from numpy.linalg import inv, eig
import matplotlib.pyplot as plt
import matplotlib.animation as ani
[docs]def draw(frame, mesh, T_trajectory, Tlims):
"""
Draws the temperature profile.
Args:
frame (int): index of the frame to be drawn
mesh (array): x-coordinates
T_trajectory (list): list of arrays of temperatures (at different times)
Tlims (list): the bottom and top values of the temperature axis: [Tmin, Tmax]
"""
plt.clf()
ax = plt.axes()
ax.set_xlim([mesh[0],mesh[-1]])
ax.set_ylim(Tlims)
plt.plot(mesh, T_trajectory[frame], 'o-', markersize=2 )
[docs]def animate(mesh, T_trajectory):
"""
Animates the development of the temperature profile using :meth:`draw`.
Args:
mesh (array): x-coordinates
T_trajectory (list): list of arrays of temperatures (at different times)
"""
nframes = len(T_trajectory)
print("animating "+str(nframes)+" frames")
Tmax = T_trajectory[0][0]
Tmin = T_trajectory[0][0]
for Ts in T_trajectory:
for t in Ts:
if t > Tmax:
Tmax = t
if t < Tmin:
Tmin = t
Tmax = round(Tmax * 0.12, 0)*10
Tmin = round(Tmin * 0.08, 0)*10
fig = plt.figure()
motion = ani.FuncAnimation(fig, draw, nframes, interval=100, fargs=(mesh, T_trajectory, [Tmin,Tmax]) )
plt.show()
[docs]def show_history(mesh, T_trajectory, last=0):
"""
Plots several temperature profiles, from different simulation times, in the same plot.
Args:
mesh (array): x-coordinates
T_trajectory (list): list of arrays of temperatures (at different times)
last (int): only plot this many temperature profiles from the end of the simulation
(by default, all data is plotted)
"""
for values in T_trajectory[-last:]:
plt.plot(mesh, values, 'o-', markersize=2)
plt.show()
[docs]def read_mesh_from_file(filename):
"""
Read the mesh and initial conditions from a file.
Each line in the file defines one mesh point containing
the position, temperature, heating power and heat conductance,
separated by a space.
That is, the file format is:
.. code-block ::
x1 T1 rho1 k1to2
x2 T2 rho2 k2to3
x3 T3 rho3 k3to4
...
where
* x is the x coordinate
* T is the initial temperature at that point
* rho is the local heating power density at that point
* k is the heat conductance *between this point and the next point*
The first and last values of rho are meaningless, since they will
always be overwritten by the boundary conditions.
The last k value is also meaningless, because there is one less
segment between points than there are points.
The points must be ordered so that the x values always increase.
Args:
filename (str): the file to read
Returns:
tuple: 4 data arrays: mesh, temperatures, power, conductance
"""
f = open(filename)
lines = f.readlines()
f.close()
mesh = []
temperatures = []
power = []
conductance = []
for line in lines:
parts = line.split()
if len(parts) > 0:
mesh.append( float(parts[0]) )
temperatures.append( float(parts[1]) )
power.append( float(parts[2]) )
conductance.append( float(parts[3]) )
return np.array( mesh ), np.array( temperatures ),\
np.array( power ), np.array( conductance )
[docs]def write_temperatures_to_file(mesh, temperatures, filename):
"""
Writes the calculated temperatures to a file.
The data is written so that each line represents one
datapoint, listing first the position and then the temperature.
Args:
mesh (array): x coordinates
temperatures (array): temperatures
filename (str): name of the file to write
"""
writelines = ""
for value, point in zip(temperatures, mesh):
writelines += str(point)+" "+str(value)+"\n"
f = open(filename,'w')
f.write(writelines)
f.close()
[docs]def calculate_stiffness_matrix(mesh, cond,
left_boundary_fixed = True,
right_boundary_fixed = True,
printout=False):
"""
Calculates and returns the stiffness matrix S.
S is a matrix whose elements
are the overlap integrals of the derivatives of the FEM basis functions.
We define :math:`\\varphi_n(x)` to be the basis function that is
1 at mesh point :math:`x_n` and zero everywhere else.
We also denote thermal conductivity by :math:`k(x)`.
The elements of the stiffness matrix are then defined as
.. math ::
S_{i,j} = \\int_0^L k(x) \\varphi'_i(x) \\varphi'_j(x) d x
except where boundary conditions override this definition.
.. note ::
This function is incomplete!
Args:
mesh (array): x coordinates of the mesh points
cond (array): heat diffusivities k for each interval between mesh points
left_boundary_fixed (bool): if True, the temperature is fixed at x = 0,
if False, temperature gradient is fixed
right_boundary_fixed (bool): if True, the temperature is fixed at x = L,
if False, temperature gradient is fixed
printout (bool) : If True, the matrix is printed on screen.
Returns:
array: the S matrix :math:`S`
"""
n_mesh = len(mesh)
S = np.zeros( [n_mesh, n_mesh] )
# todo
if printout:
# print S with rounded values
# choose rounding precision according to the smallest non-zero element
maximum = np.amax( np.abs(S) )
minimum = np.amin( np.where( np.abs(S)>0, np.abs(S), maximum ) )
decimals = int( np.ceil( -np.log10( minimum ) ) )
print( "" )
print( "stiffness matrix" )
print( np.round( S, decimals ) )
print( "" )
return S
[docs]def calculate_mass_matrix(mesh,
left_boundary_fixed = True,
right_boundary_fixed = True,
printout=False):
"""
Calculates and returns the mass matrix M.
M is a matrix whose elements
are the overlap integrals of the FEM basis functions.
We define :math:`\\varphi_n(x)` to be the basis function that is
1 at mesh point :math:`x_n` and zero everywhere else.
The elements of the stiffness matrix are then defined as
.. math ::
M_{i,j} = \\int_0^L \\varphi_i(x) \\varphi_j(x) d x
except where boundary conditions override this definition.
.. note ::
This function is incomplete!
Args:
mesh (array): x coordinates of the mesh points
left_boundary_fixed (bool): if True, the temperature is fixed at x = 0,
if False, temperature gradient is fixed
right_boundary_fixed (bool): if True, the temperature is fixed at x = L,
if False, temperature gradient is fixed
printout (bool) : If True, the matrix is printed on screen.
Returns:
array: the M matrix :math:`M`
"""
n_mesh = len(mesh)
M = np.zeros( [n_mesh, n_mesh] )
# todo
if printout:
# print M with rounded values
# choose rounding precision according to the smallest non-zero element
maximum = np.amax( np.abs(M) )
minimum = np.amin( np.where( np.abs(M)>0, np.abs(M), maximum ) )
decimals = int( np.ceil( -np.log10( minimum ) ) )
print( "" )
print( "mass matrix" )
print( np.round( M, decimals ) )
print( "" )
return M
[docs]def calculate_power_vector(mesh, temperatures, heat, cond, M,
left_boundary_fixed, right_boundary_fixed,
printout = False):
"""
Calculates the power vector :math:`p`.
For a fixed boundary (temperature is kept constant):
* Replace the boundary node values of the heating power density vector
:math:`\\rho` with the temperature values at the boundary.
You get :math:`\\rho = [T(0), \\rho(x_1), \\rho(x_2), \\ldots, \\rho(x_{N}), T(L)]`.
* Calculate the power vector as :math:`p = M \\rho`.
For a fixed gradient boundary (heat flux is kept constant):
* Replace the boundary node values of the heating power density vector
:math:`\\rho` with zeros.
You get :math:`\\rho = [0, \\rho(x_1), \\rho(x_2), \\ldots, \\rho(x_{N}), 0]`.
* Calculate temperature derivatives at the boundaries, :math:`T'(0)` and :math:`T'(L)`.
This is done with :meth:`calculate_left_temperature_gradient` and
:meth:`calculate_right_temperature_gradient`.
* Define the heat flux vector :math:`d` so that it has zeros everywhere
except at the boundaries where it contains the value of the heat
flux flowing out of the system. At :math:`x = 0`, this value is
:math:`k_{0,1} T'(0)`, and at :math:`x = L`, it is :math:`-k_{N,N+1}T'(L)`.
You get :math:`d = [k_{0,1}T'(0), 0, 0, \\ldots, 0, -k_{N,N+1}T'(L))]`.
* Calculate the power vector as :math:`p = M \\rho - d`.
For a mix of different boundary conditions, the correct procedure is applied
separately at both boundaries (both ends of the vector).
.. note ::
This function is incomplete!
Args:
mesh (array): x coordinates
temperatures (array): temperatures
heat (array): heating power densities :math:`\\rho`
cond (array): heat conductivities :math:`k`
M (array): the mass matrix :math:`M`
left_boundary_fixed (bool): Specifiest the boundary condition at the left edge.
If True, a constant temperature is applied.
If False, a contant temperature gradient is applied.
right_boundary_fixed (bool): Specifiest the boundary condition at the right edge.
If True, a constant temperature is applied.
If False, a contant temperature gradient is applied.
printout (bool) : If True, the matrix is printed on screen.
Returns:
array: power vector :math:`p`
"""
n_mesh = len(mesh)
power = np.zeros( n_mesh )
# todo
if printout:
print( "" )
print( "power vector" )
print( np.round( power, 3 ) )
print( "" )
return power
[docs]def calculate_propagator(S, M, dt):
"""
Calculates the propagator matrix for a dynamic simulation.
Time is advanced according to the equation
.. math ::
T(t + \\Delta t) = P T(t) + L.
This function calculates the matrix :math:`P`.
Args:
S (array): the stiffness matrix :math:`S`
M (array): the mass matrix :math:`M`
dt (float): time step :math:`\\Delta t`
Returns:
array: the propagator matrix :math:`P`
"""
return inv(M+0.5*dt*S) @ (M-0.5*dt*S)
[docs]def calculate_load(S, M, p, dt):
"""
Calculates the load vector for a dynamic simulation.
Time is advanced according to the equation
.. math ::
T(t + \\Delta t) = P T(t) + L.
This function calculates the vector :math:`L`.
Args:
S (array): the stiffness matrix :math:`S`
M (array): the mass matrix :math:`M`
p (array): the power vector :math:`p`
dt (float): time step :math:`\\Delta t`
Returns:
array: the load vector :math:`L`
"""
return inv(M+0.5*dt*S) @ (dt*p)
[docs]def calculate_left_temperature_gradient(mesh, temperatures):
"""
Calculates and returns the temperature derivative at the left boundary, :math:`T'(0)`.
The derivative is calculated using a second degree numeric estimate.
That is, the function calculates the slope of the parabola that
passes through points :math:`(x_0, T(x_0)), (x_1, T(x_1)), (x_2, T(x_2))`,
evaluated at :math:`x_0`.
Args:
temperatures (array): temperatures
mesh (array): x coordinates
Returns:
float: numeric approximation for the derivative
"""
t3 = temperatures[2]
t2 = temperatures[1]
t1 = temperatures[0]
x3 = mesh[2]
x2 = mesh[1]
x1 = mesh[0]
return (x3**2 * (t1-t2) + 2*x1*(x3*(-t1+t2) + x2*(t1-t3)) + \
x2**2 * (-t1+t3) + x1**2 * (-t2+t3))/((x1-x2)*(x1-x3)*(x2-x3))
[docs]def calculate_right_temperature_gradient(mesh, temperatures):
"""
Calculates and returns the temperature derivative at the right boundary, :math:`T'(L)`.
The derivative is calculated using a second degree numeric estimate.
That is, the function calculates the slope of the parabola that
passes through points :math:`(x_{N-1}, T(x_{N-1})), (x_{N}, T(x_{N})), (x_{N+1}, T(x_{N+1}))`,
evaluated at :math:`x_{N+1}` (assuming :math:`x_{N+1} = L`).
Args:
temperatures (array): temperatures
mesh (array): x coordinates
Returns:
float: numeric approximation for the derivative
"""
t3 = temperatures[-1]
t2 = temperatures[-2]
t1 = temperatures[-3]
x3 = mesh[-1]
x2 = mesh[-2]
x1 = mesh[-3]
return (x3**2 * (-t1+t2) + 2*x2*x3*(t1-t3) + x1**2*(t2-t3) + \
x2**2*(-t1+t3) + 2*x1*x3*(-t2+t3))/((x1-x2)*(x1-x3)*(x2-x3))
[docs]def calculate_total_heating_power(mesh, heat):
"""
Integrates the heating power density with respect to position.
At continuum level, heating power density is a continuous quantity
and the total heating power is given by
.. math::
P = \\int_0^L \\rho(x) dx.
We only know the heating power density at discrete points, so the
function calculates the approximation
.. math::
P = \\sum_{i = 1}^{N} \\rho(x_i) \\frac{1}{2}(x_{i-1} - x_{i+1}).
Args:
mesh (array): x coordinates
heat (array): heating power densities :math:`\\rho`
Returns:
float: the total heating power
"""
heat_power = 0.0
# Integrate over the mesh.
# Ignore the end points of the mesh since the power
# density is modified there due to boundary conditions.
for i in range(1,len(mesh)-1):
# For each mesh point, take the power density at that
# point and calculate the area of the triangle whose
# height is the power density and width is the mesh
# width at that point.
# This is the heating power given by the power density
# at that mesh point.
heat_power += 0.5*heat[i]*(mesh[i+1] - mesh[i-1])
return heat_power
[docs]def run_simulation(temperatures, propagator, load, n_steps):
"""
Runs a time dependent heat simulation.
Time is advanced according to the equation
.. math ::
T(t + \\Delta t) = P T(t) + L.
The function does not change the original temperatures array.
Instead, the final temperatures are returned as a new array.
Args:
temperatures (array): temperatures
propagator (array): the propagator matrix :math:`P`
load (array): the load vector :math:`L`
n_steps (int): simulation length in number of time steps
Returns:
array: temperatures at the end of the simulation
"""
for i in range(n_steps):
temperatures = propagator @ temperatures + load
return temperatures
[docs]def calculate_heat_flux(mesh, temperatures, cond, printout=False):
"""
Calculates the total heat flux into the system at the system boundaries.
Heat flux measures the power of energy transfer via heat conductance.
It is proportional to thermal conductivity and the temperature gradient.
In a 1-dimensional system, it is simply
.. math ::
\\Phi_Q = -k T'(x).
Heat flux is a signed quantity. If energy flows in positive x direction,
flux is positive. If energy flows in negative x direction, flux is negative.
If there is non-zero heat flux at the system boundary, energy flows between
the system and the environment.
.. note ::
This function is incomplete!
Args:
mesh (array): x coordinates
temperatures (array): temperatures
cond (array): conductivities :math:`k`
Returns:
float, float: heat flux into the system at left boundary,
heat flux into the system at right boundary
"""
# Calculate the total heat flux *into* the system as
# P_x = -k dT/dx. Note that if this is positive, heat flows
# to positive x-direction. At the left boundary, this is into
# the system but at the right boundary this is out of the system.
# todo
heat_flux_left = 0.0
heat_flux_right = 0.0
heat_flux = heat_flux_left - heat_flux_right
if printout:
print( "heat flux from left ", np.round(heat_flux_left,4) )
print( "heat flux from right ", np.round(-heat_flux_right,4) )
print( "total heat flux ", np.round(heat_flux,4) )
print( "" )
return heat_flux_left, -heat_flux_right
[docs]def main(filename='mesh1.txt', dt=0.001, recording_dt=0.01, simulation_time=1.0,
left_boundary_fixed=True, right_boundary_fixed=True, printout=False):
"""
Run a complete simulation and save results at regular intervals.
Args:
filename (str): system file to read
dt (float): time step :math:`\\Delta t`
recording_dt (float): time between data recording
simulation_time (float): total time to simulate
"""
#
# Read mesh and system parameters from file.
#
# The resulting vectors are:
# * mesh: x-coordinates of mesh points
# * temperatures: initial temperatures T(x,0) at mesh points
# * heat: heating power density rho(x) at mesh points
# * cond: thermal conductivity *between* mesh points k(x,x+dx)
#
# Note that 'cond' contains data for elements between mesh
# points so cond[0] is the conductivity between points
# 0 and 1, cond[1] is the conductivity between points 1 and 2
# and so on.
#
mesh, temperatures, heat, cond = read_mesh_from_file(filename)
# simulation parameters
interval_steps = int(recording_dt / dt)
records = int(simulation_time / recording_dt)
# calculate matrices
S = calculate_stiffness_matrix(mesh, cond,
left_boundary_fixed,
right_boundary_fixed,
printout)
M = calculate_mass_matrix(mesh,
left_boundary_fixed,
right_boundary_fixed,
printout)
p = calculate_power_vector(mesh, temperatures, heat, cond, M,
left_boundary_fixed, right_boundary_fixed,
printout)
# calculate propagator and load
propagator = calculate_propagator(S, M, dt)
load = calculate_load(S, M, p, dt)
# list for storing the temperature profiles at different times
history = [temperatures]
# run the dynamic simulation
#
# the simulation is run in short sequences
# and temperatures are saved between each sequence
for i in range(records):
temperatures = run_simulation(temperatures, propagator, load, interval_steps)
history.append(temperatures)
print( "final temperatures ", np.round(temperatures,2) )
# calculate energy transfer
phiL, phiR = calculate_heat_flux(mesh, temperatures, cond, printout)
P = calculate_total_heating_power(mesh, heat)
print( "total energy change rate "+str( round(phiL + phiR + P, 3) ) )
# visualize
show_history(mesh, history)
animate(mesh, history)
if __name__ == "__main__":
# choose simulation settings
system = 'A5_rod_uniform.txt'
time = 2
left_edge_constant_T = True
right_edge_constant_T = True
# Suppress scientific notation of numbers
# This makes printed matrices easier to read.
np.set_printoptions(suppress=True)
main(printout=True, filename=system, simulation_time = time,
left_boundary_fixed = left_edge_constant_T, right_boundary_fixed = right_edge_constant_T )