Clone wiki

Tis / Home


Tis is a C++ library implementing explicit finite-difference scheme for the solution of 1d heat conduction equation. Interval between \(x_{left}\) and \(x_{right}\) with specified boundary conditions is considered. Additional internal boundaries inside the interval can be added, e.g. contacts of different materials. Several kinds of end and middle boundary conditions are available, including given value, flux, linear, radiation exchange. The material properties can vary with position and temperature. The solution is obtained with the time step necessary for the stability of the scheme.


The library integrates 1d heat conduction equation:

\begin{equation*} \rho c \frac{\partial T}{\partial t} = \frac{\partial}{\partial x} \left(K\frac{\partial T}{\partial x}\right) + S = K\frac{\partial^2 T}{\partial x^2} + \frac{\partial K}{\partial x} \cdot \frac{\partial T}{\partial x} + S \end{equation*}

in the interval:

\begin{align*} x_{left} \leq x \leq x_{right}\\ 0 \leq t \leq t_{fin} \end{align*}

with appropriate boundary and initial conditions. Additional boundaries \(x_{m,i}\) inside the integration interval can be presented with contact conditions specified.

\begin{equation*} x_{left} < x_{m,i} < x_{right},\ i=0, \ldots N-1 \end{equation*}

A contact condition returns two boundary values for the both side of the contact. Material constants \(\rho,\ c, \ K\) are functions of position and temperature. The source term \(S\) depends on time, position and temperature.

The equation is integrated with finite differences on a uniform \(t-x\) grid. The solution is \(T_{i,j}\), where \(i\) is time index and \(j\) is space index. The explicit finite difference scheme is obtained by using the following approximations of the derivatives (\(\tau\) - time step, \(h\) - space step):

\begin{align*} \left(\frac{\partial T}{\partial t}\right)_{i,j} &= \frac{T_{i+1,j} - T_{i,j}}{\tau} \\ \left(\frac{\partial^2 T}{\partial x^2}\right)_{i,j} &= \frac{T_{i,j-1} - 2T_{i,j}+ T_{i,j+1}}{h^2} \\ \left(\frac{\partial K}{\partial x}\right)_{i,j} &= \frac{K_{i,j+1} - K_{i,j-1}}{2h} \\ \left(\frac{\partial T}{\partial x}\right)_{i,j} &= \frac{T_{i,j+1}^t - T_{i,j-1}^t}{2h} \\ \end{align*}

As a result the scheme is:

\begin{equation*} T_{i+1,j} = T_{i,j} + \frac{\tau}{h^2} \cdot \left( T_{i,j-1} \cdot \left(\kappa_0 - \frac{\Delta\kappa}{4}\right) - T_{i,j}\cdot 2\kappa_0 + T_{i,j+1}^t\cdot\left(\kappa_0 + \frac{\Delta\kappa}{4}\right) \right) + \tau\cdot\frac{S_{i,j}}{(\rho c)_{i,j}} \end{equation*}


\begin{align*} \kappa_0 &\equiv \frac{K_{i,j}}{(\rho c)_{i,j}}\\ \Delta\kappa &\equiv \frac{K_{i,j+1} - K_{i,j-1}}{(\rho c)_{i,j}}\\ \end{align*}

The scheme is stable provided:

\begin{equation*} \left(\frac{K}{\rho c}\right)_{max} \cdot \frac{\tau}{h^2} < \frac{1}{2} \end{equation*}

and \(\Delta\kappa\) is sufficiently small, i.e. \(h\) is fine enough. The scheme is used to obtain \(T_{i+1,j}\) values at the new time slice for all the points from \(j=1\) to \(j=N-2\). The most left and right values are found from the corresponding boundary conditions. The process is continued until the desired time point.


The top level class solving the 1d heat conduction equation is ExplicitNSolver. To construct an instance of the class one has to specify left and right \(x\) boundaries with boundary conditions and the final time to integrate to.

#include <Tis/ExplicitNSolver.h>
#include <Tis/BoundaryCondition.h>
#include <Tis/Source.h>
#include <Tis/Material.h>

FixedValueLeft cond_left(1.0);
FixedValueRight cond_right(1.0);
ExplicitNSolver solver(0.0, cond_left, 1e-2, cond_right, 1.0);

Here the integration is to be done between \(0\) and \(1\) cm, from \(0\) to \(1\) second. On both sides of the interval a fixed temperature of \(1\) is given. Before integrating it is necessary to specify the material, i.e. \((K, \rho, c)\), source function, spatial resolution and initial values. The time resolution is chosen from the requirement of the scheme stability with given safety margin \(\tau = f \cdot h^2 \cdot \left(\frac{K}{\rho c}\right)_{max}^{-1}\), by default \(f=0.25\) and can not be more than \(0.5\).

ZeroSource source; // simply returns 0
ConstantMaterial material(1, 1, 1); // K = 1, rho = 1, c = 1
vector<double> initvalues(101, 0.0);
solver.solve(); // do the integration

To add a contact between two materials inside the interval it is necessary to specify the boundary position and the contact condition.

IdealContact contact; // equal flux and temperatures
solver.addBoundary(0.5, contact);

In addition to that it is now necessary to add source, material, resolution and initial values for the second sub-interval. The previously set values are attributed to the region from \(0\) to \(0.5\) only.

ConstantMaterial material2(2, 1, 1); // K = 2, rho = 1, c = 1
vector<double> initvalues2(51, 0.0);
solver.solve(); // do the integration

The solution of the problem is stored as a vector of solutions on sub-intervals. The latter is of type IntervalSolution, which contains \(x\) vector, \(t\) vector and a 2d matrix of \(T_{ij}\). For the matrix time is constant in a row, i.e. \(T_{0,j}\) selects time step \(0\).

const IntervalSolution &solution = solver.getSolutions()[0]; // first sub-interval
solution.getXValue(50); // get x value at index 50
solution.getTimeValue(10); //get time value at index 10
solution.getValue(10, 50); //get temperature at time 10 and x 50
solution(10, 50); //the same as above

There are other methods to investigate the interval, e.g. finding the size, resolution, etc. It is also possible to find the temperature value from the solver.

solver(0, 10, 50); // equivalent to solution(10, 50) for the first sub-interval

To save the solution to disk use dumpSolution method. This saves to a specified directory \(x\), \(t\) vectors and solutions matrices for all sub-intervals.

solver.dumpSolution("run1"); //save every time step
solver.dumpSolution("run1_short", 10); // save every tenth time step

To reset any integrated solutions and any set material, sources, etc. use reset method.

Boundary conditions

There are several left/right boundary conditions available.

FixedValueLeft, FixedValueRight
Given temperature on the left/right.
FixedFluxLeft, FixedFluxRight
Given flux \(-K\frac{\partial T}{\partial x}\) on the left/right.
LinearExchangeLeft, LinearExchangeRight
Linear heat exchange with a medium at fixed temperature. \(-K\frac{\partial T}{\partial x} = H\cdot(T-T_0)\).
RadiationLeft, RadiationRight
Radiation exchange with a medium at fixed temperature. \(-K\frac{\partial T}{\partial x} = \varepsilon\sigma\cdot(T^4-T_0^4)\).

It is possible to add a new boundary condition by inheriting LeftBoundaryCondition or RightBoundaryCondition and implementing getBoundaryValue method. The method is called with the time step index, an appropriate Material instance and an IntervalSolution instance. The method has to return the found boundary value.

In the available boundary conditions the temperature derivatives are approximated with the second order accuracy.

\begin{align*} \left(\frac{\partial T}{\partial x}\right)_{0} &= \frac{4 T_{0+1} - T_{0+2} -3T_{0}}{2h}\\ \left(\frac{\partial T}{\partial x}\right)_{N-1} &= \frac{T_{N-3} - 4T_{N-2} +3T_{N-1}}{2h}\\ \end{align*}

There are also several available contact conditions.

Fixed temperature at the contact.
Equal flux and equal temperatures. \(-K_1\frac{\partial T_1}{\partial x} = -K_2\frac{\partial T_2}{\partial x}\), \(T_1 = T_2\)
Equal flux and flux is proportional to the temperature difference. \(-K_1\frac{\partial T_1}{\partial x} = -K_2\frac{\partial T_2}{\partial x}\), \(-K_1\frac{\partial T_1}{\partial x} = H\cdot(T_1 - T_2)\)
Equal flux and flux is given by Stefan-Boltzmann law. \(-K_1\frac{\partial T_1}{\partial x} = -K_2\frac{\partial T_2}{\partial x}\), \(-K_1\frac{\partial T_1}{\partial x} = \epsilon_{12}\sigma\cdot(T_1^4 - T_2^4)\), \(\epsilon_{12} = \frac{\epsilon_1\epsilon_2}{\epsilon_1+\epsilon_2 -\epsilon_1\epsilon_2}\)
Radiation exchange between surfaces in limited solid angle and in the remaining solid angle radiative exchange with an outside medium. \(-K_1\frac{\partial T_1}{\partial x} = f\cdot\epsilon_{12}\sigma\cdot(T_1^4 - T_2^4)+ (1-f)\cdot\epsilon_1\sigma\cdot(T_1^4 - T_0^4)\), \(-K_2\frac{\partial T_2}{\partial x} = f\cdot\epsilon_{12}\sigma\cdot(T_1^4 - T_2^4)- (1-f)\cdot\epsilon_2\sigma\cdot(T_2^4 - T_0^4)\)

To add a new contact condition one has to inherit from MidBoundaryCondition and implement method getBoundaryValue. The method is called with the index of the current time step, materials on the left and right of the contact, and with solutions on the left and on the right. The method has to return a pair of temperatures: one for the left sub-interval and one for the right one.


Available source classes are:

No source.
Imitates radiation exchange along sides of a rod. This is equivalent to the volume source \(\epsilon\sigma\cdot(T_0^4 - T^4)\cdot \frac{L}{S}\), where \(L\) is perimeter of the rod cross-section and \(S\) is its area.

To add a new source function inherit Source and implement getSource method, which is called with current time, position and temperature.


A material class has to inherit Material and to implement the corresponding methods getHeatConductivity, getHeatCapacity and getDensity. These methods are called with current position and temperature.

Available are:

Fixed values independent on anything.
Data for stainless steel, graphite, CuCrZr
Available in MaterialSteel.h, MaterialGraphite.h, MaterialCuCrZr.h.

Build instructions

Installation is standard for cmake projects.

>>> mkdir build && cd build
>>> cmake ..
>>> make
>>> make install

To build tests Boost.Tests is required.

>>> make tests
>>> tests/tests --log_level=test_suite

The library has an extensive test suite that validates it against known solutions.