Homework: Staying warm

Theory

In one spatial dimension, heat equation is the partial differential equation

\[\frac{\partial u}{\partial t} = k \frac{\partial^2 u}{\partial x^2} + \rho.\]

Here \(u\) may be either thermal energy density or the temperature (these are directly proportional to each other), \(k\) is thermal conductivity or heat diffusivity (depending on what \(u\) represents) and \(\rho\) is internal heating power density (or a quantity proportional to it). Roughly speaking, the equation says that temperatures change if there is heating (or cooling) or if there is a changing temperature gradient.

Temperature gradient affects temperature because it generates a heat flux, i.e., spontaneous flow of energy as heat. In one dimension, heat flux is defined as the power transported by the flowing energy,

\[P = -k \frac{\partial T}{\partial x}.\]

Flux is positive, if energy flows in positive \(x\)-direction, which happens if \(\frac{\partial T}{\partial x} < 0\), hence the minus sign.

If heat flux is the same everywhere in a system, energy simply flows through the system at a constant rate. That means energy density does not change and neither does temperature. However, if the heat flux changes as one moves from one place to another, it means that energy flows at different rates at different parts of the system. If there is a point where more energy flows in than out, there is energy accumulation at that point, and temperature rises. (And similarly temperature drops if energy flowa out.) Mathematically, such conditions arise if the heat flux changes as a function of position, i.e., if

\[\frac{\partial P}{\partial x} = -k \frac{\partial^2 T}{\partial x^2} \ne 0.\]

This is the second degree term of the heat equation.

In this task, we wish to simulate the heat equation using the finite element method (FEM). We pick a mesh of points \(x_i\) and represent the temperature at these points as the vector

\[\boldsymbol{u}(t) = [ u(x_0, t), u(x_1, t), \ldots, u(x_{N+1}, t) ]^T.\]

As derived during the lectures, the time evolution of this vector is given by

\[\boldsymbol{u}(t + \Delta t) = \boldsymbol{P}\boldsymbol{u}(t) + \boldsymbol{L},\]

where the propagator \(\boldsymbol{P}\) and load \(\boldsymbol{L}\) depend on the so called stiffness and mass matrices \(\boldsymbol{S}\) and \(\boldsymbol{M}\). Generally, the stiffness matrix elements are

\[S_{i,j} = \int k(x) \varphi_i'(x) \varphi_j'(x) \mathrm{d} x\]

and mass matrix elements

\[M_{i,j} = \int \varphi_i(x) \varphi_j(x) \mathrm{d} x\]

where \(\varphi_i(x)\) is the piecewise linear function that is 1 at \(x_i\) and 0 at all other mesh points. Boundary conditions may override these definitions for the elements of the first and last row and column.

Task

To complete this assignment, you need to complete the following tasks, submit your code as a plain .py file, and answer the questions. In the computing tasks, you must implement the incomplete functions so that they operate exactly as described in the documentation.

  • Computing task:

    1. Implement the calculation of matrices \(\boldsymbol{S}\), \(\boldsymbol{M}\) and \(\boldsymbol{p}\). We derived these matrices for a uniform mesh during the lecture, but now you must calculate the more general form where the distance between mesh points may vary.

    2. Implement boundary conditions for fixed temperatures and fixed temperature gradients. Fix the values to be the same as in the initial temperature distribution.

    3. Implement calculation of heat flux at system boundaries.

    4. Simulate temperature evolution in the four given systems. Apply different boundary conditions and check how that affects the result. Calculate heat flux at system boundaries and the total heating power.

  • Questions to answer:

    1. Explain how you solved the necessary matrices.

    2. The final temperature profile of the homogeneous rod is a straight line. What kind of a profile do you get for the insulating layer? Explain in your own words, why the temperature behaves in this way.

    3. There are systems that never reach a steady state. What kinds of conditions can lead to this? Include plots of simulated temperatures in your answer.

    4. Study the boundary heat flux and total heating power in systems that reach a steady state and in systems that don’t. Explain, in your own words, what condition for the flux and power is necessary for a steady state, and give a physical explanation to why it must be so. Again, support your reasoning with simulation results.

A5_fem.py

A5_fem.animate(mesh, T_trajectory)[source]

Animates the development of the temperature profile using draw().

Parameters:

  • mesh (array) – x-coordinates

  • T_trajectory (list) – list of arrays of temperatures (at different times)

A5_fem.calculate_heat_flux(mesh, temperatures, cond, printout=False)[source]

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

\[\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!

Parameters:

  • mesh (array) – x coordinates

  • temperatures (array) – temperatures

  • cond (array) – conductivities \(k\)

Returns:

heat flux into the system at left boundary,

heat flux into the system at right boundary

Return type:

float, float

A5_fem.calculate_left_temperature_gradient(mesh, temperatures)[source]

Calculates and returns the temperature derivative at the left boundary, \(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 \((x_0, T(x_0)), (x_1, T(x_1)), (x_2, T(x_2))\), evaluated at \(x_0\).

Parameters:

  • temperatures (array) – temperatures

  • mesh (array) – x coordinates

Returns:

numeric approximation for the derivative

Return type:

float

A5_fem.calculate_load(S, M, p, dt)[source]

Calculates the load vector for a dynamic simulation.

Time is advanced according to the equation

\[T(t + \Delta t) = P T(t) + L.\]

This function calculates the vector \(L\).

Parameters:

  • S (array) – the stiffness matrix \(S\)

  • M (array) – the mass matrix \(M\)

  • p (array) – the power vector \(p\)

  • dt (float) – time step \(\Delta t\)

Returns:

the load vector \(L\)

Return type:

array

A5_fem.calculate_mass_matrix(mesh, left_boundary_fixed=True, right_boundary_fixed=True, printout=False)[source]

Calculates and returns the mass matrix M.

M is a matrix whose elements are the overlap integrals of the FEM basis functions.

We define \(\varphi_n(x)\) to be the basis function that is 1 at mesh point \(x_n\) and zero everywhere else. The elements of the stiffness matrix are then defined as

\[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!

Parameters:

  • 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:

the M matrix \(M\)

Return type:

array

A5_fem.calculate_power_vector(mesh, temperatures, heat, cond, M, left_boundary_fixed, right_boundary_fixed, printout=False)[source]

Calculates the power vector \(p\).

For a fixed boundary (temperature is kept constant):

  • Replace the boundary node values of the heating power density vector \(\rho\) with the temperature values at the boundary. You get \(\rho = [T(0), \rho(x_1), \rho(x_2), \ldots, \rho(x_{N}), T(L)]\).

  • Calculate the power vector as \(p = M \rho\).

For a fixed gradient boundary (heat flux is kept constant):

  • Replace the boundary node values of the heating power density vector \(\rho\) with zeros. You get \(\rho = [0, \rho(x_1), \rho(x_2), \ldots, \rho(x_{N}), 0]\).

  • Calculate temperature derivatives at the boundaries, \(T'(0)\) and \(T'(L)\). This is done with calculate_left_temperature_gradient() and calculate_right_temperature_gradient().

  • Define the heat flux vector \(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 \(x = 0\), this value is \(k_{0,1} T'(0)\), and at \(x = L\), it is \(-k_{N,N+1}T'(L)\). You get \(d = [k_{0,1}T'(0), 0, 0, \ldots, 0, -k_{N,N+1}T'(L))]\).

  • Calculate the power vector as \(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!

Parameters:

  • mesh (array) – x coordinates

  • temperatures (array) – temperatures

  • heat (array) – heating power densities \(\rho\)

  • cond (array) – heat conductivities \(k\)

  • M (array) – the mass matrix \(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:

power vector \(p\)

Return type:

array

A5_fem.calculate_propagator(S, M, dt)[source]

Calculates the propagator matrix for a dynamic simulation.

Time is advanced according to the equation

\[T(t + \Delta t) = P T(t) + L.\]

This function calculates the matrix \(P\).

Parameters:

  • S (array) – the stiffness matrix \(S\)

  • M (array) – the mass matrix \(M\)

  • dt (float) – time step \(\Delta t\)

Returns:

the propagator matrix \(P\)

Return type:

array

A5_fem.calculate_right_temperature_gradient(mesh, temperatures)[source]

Calculates and returns the temperature derivative at the right boundary, \(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 \((x_{N-1}, T(x_{N-1})), (x_{N}, T(x_{N})), (x_{N+1}, T(x_{N+1}))\), evaluated at \(x_{N+1}\) (assuming \(x_{N+1} = L\)).

Parameters:

  • temperatures (array) – temperatures

  • mesh (array) – x coordinates

Returns:

numeric approximation for the derivative

Return type:

float

A5_fem.calculate_stiffness_matrix(mesh, cond, left_boundary_fixed=True, right_boundary_fixed=True, printout=False)[source]

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 \(\varphi_n(x)\) to be the basis function that is 1 at mesh point \(x_n\) and zero everywhere else. We also denote thermal conductivity by \(k(x)\). The elements of the stiffness matrix are then defined as

\[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!

Parameters:

  • 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:

the S matrix \(S\)

Return type:

array

A5_fem.calculate_total_heating_power(mesh, heat)[source]

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

\[P = \int_0^L \rho(x) dx.\]

We only know the heating power density at discrete points, so the function calculates the approximation

\[P = \sum_{i = 1}^{N} \rho(x_i) \frac{1}{2}(x_{i-1} - x_{i+1}).\]
Parameters:

  • mesh (array) – x coordinates

  • heat (array) – heating power densities \(\rho\)

Returns:

the total heating power

Return type:

float

A5_fem.draw(frame, mesh, T_trajectory, Tlims)[source]

Draws the temperature profile.

Parameters:

  • 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]

A5_fem.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)[source]

Run a complete simulation and save results at regular intervals.

Parameters:

  • filename (str) – system file to read

  • dt (float) – time step \(\Delta t\)

  • recording_dt (float) – time between data recording

  • simulation_time (float) – total time to simulate

A5_fem.read_mesh_from_file(filename)[source]

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:

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.

Parameters:

filename (str) – the file to read

Returns:

4 data arrays: mesh, temperatures, power, conductance

Return type:

tuple

A5_fem.run_simulation(temperatures, propagator, load, n_steps)[source]

Runs a time dependent heat simulation.

Time is advanced according to the equation

\[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.

Parameters:

  • temperatures (array) – temperatures

  • propagator (array) – the propagator matrix \(P\)

  • load (array) – the load vector \(L\)

  • n_steps (int) – simulation length in number of time steps

Returns:

temperatures at the end of the simulation

Return type:

array

A5_fem.show_history(mesh, T_trajectory, last=0)[source]

Plots several temperature profiles, from different simulation times, in the same plot.

Parameters:

  • 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)

A5_fem.write_temperatures_to_file(mesh, temperatures, filename)[source]

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.

Parameters:

  • mesh (array) – x coordinates

  • temperatures (array) – temperatures

  • filename (str) – name of the file to write