One-Dimensional Heat Transfer¶

Deriving Equations¶

Before getting into the math, there is an important point worth making. Modelica is a language for representing lumped systems. What this means is that the behavior must be expressed in terms of ordinary differential equations (ODEs) or in some cases differential-algebraic equations (DAEs). But Modelica does not include any means for describing partial differential equations (i.e., equations that involve the gradient of variables in spatial directions). As such, in this section we will derive the ordinary differential equations that result when we discretize a rod into discrete pieces and then model the behavior of each of these discrete (lumped) pieces.

With that out of the way, let us consider the heat balance for each discrete section of this rod. First, we have the thermal capacitance of the section. This can expressed as:

where is the mass of the section, is the capacitance (a material property) and is the temperature of the section. We can further describe the mass as:

where is the material density and is the volume of the section. Finally, we can express the volume of the section as:

where is the cross-sectional area of the section, which is assumed to be uniform, and is the length of the section. For this example, we will assume the rod is composed of equal size pieces. In this case, we can define the segment length, , to be:

We will also assume that the cross-sectional area is uniform along the length of the rod. As such, the mass of each segment can be given as:

In this case, the thermal capacitance of each section would be:

This, in turn, means that the net heat gained in that section at any time will be:

where we assume that , and don’t change with respect to time.

That covers the thermal capacitance. In addition, we will consider two different forms of heat transfer. The first form of heat transfer we will consider is convection from each section to some ambient temperature, . We can express the amount of heat lost from each section as:

where is the convection coefficient and is the surface area of the section. The other form of heat transfer is conduction to neighboring sections. Here there will be two contributions, one lost to the section, if it exists, and the other lost to the section, if it exists. These can be represented, respectively, as:

Using these relations, we know that the heat balance for the first element would be:

Similarly, the heat balance for the last element would be:

Finally, the heat balance for all other elements would be:

Implementation¶

We start by defining types for the various physical quantities. This will give us the proper units and, depending on the tool, allows us to do unit checking on our equations. Our type definitions are as follows:

  type Temperature=Real(unit="K", min=0);
  type ConvectionCoefficient=Real(unit="W/K", min=0);
  type ConductionCoefficient=Real(unit="W.m-1.K-1", min=0);
  type Mass=Real(unit="kg", min=0);
  type SpecificHeat=Real(unit="J/(K.kg)", min=0);
  type Density=Real(unit="kg/m3", min=0);
  type Area=Real(unit="m2");
  type Volume=Real(unit="m3");
  type Length=Real(unit="m", min=0);
  type Radius=Real(unit="m", min=0);

We will also define several parameters to describe the rod we are simulating:

  parameter Integer n=10;
  parameter Length L=1.0;
  parameter Radius R=0.1;
  parameter Density rho=2.0;
  parameter ConvectionCoefficient h=2.0;
  parameter ConductionCoefficient k=10;
  parameter SpecificHeat C=10.0;
  parameter Temperature Tamb=300 "Ambient temperature";

Given these parameters, we can compute the areas and volume for each section in terms of the parameters we have already defined using the following declarations:

  parameter Area A_c = pi*R^2, A_s = 2*pi*R*L;
  parameter Volume V = A_c*L/n;

Finally, the only array in this problem is the temperature of each section (since this is the only quantity that actually varies along the length of the rod):

  Temperature T[n];

This concludes all the declarations we need to make. Now let’s consider the various equations required. First, we need to specify the initial conditions for the rod. We will assume that , and the initial temperatures of all other sections can be linearly interpolated between these two end conditions. This is captured by the following equation:

initial equation
  T = linspace(200,300,n);

where the linspace operator is used to create an array of n values that vary linearly between 200 and 300. Recall from our State Space examples that we can include equations where the left hand side and right hand side expressions are vectors. This is another example of such an equation.

Finally, we come to the equations that describe how the temperature in each section changes over time:

equation
  rho*V*C*der(T[1]) = -h*A_s*(T[1]-Tamb)-k*A_c*(T[1]-T[2])/(L/n);
  for i in 2:(n-1) loop
    rho*V*C*der(T[i]) = -h*A_s*(T[i]-Tamb)-k*A_c*(T[i]-T[i-1])/(L/n)-k*A_c*(T[i]-T[i+1])/(L/n);
  end for;
  rho*V*C*der(T[end]) = -h*A_s*(T[end]-Tamb)-k*A_c*(T[end]-T[end-1])/(L/n);

The first equation corresponds to the heat balance for section , the last equation corresponds to the heat balance for section and the middle equation covers all other sections. Note the use of end as a subscript. When an expression is used to evaluate a subscript for a given dimension, end represents the size of that dimension. In our case, we use end to represent the last section. Of course, we could use n in this case, but in general, end can be very useful when the size of a dimension is not already associated with a variable.

Also note the use of a for loop in this model. A for loop allows the loop index variable to loop over a range of values. In our case, the loop index variable is i and the range of values is through . The general syntax for a for loop is:

for <var> in <range> loop
  // statements
end for;

where <range> is a vector of values. A convenient way to generate a sequence of values is to use the range operator, :. The value before the range operator is the initial value in the sequence and the value after the range operator is the final value in the sequence. So, for example, the expression 5:10 would generate a vector with the values 5, 6, 7, 8, 9 and 10. Note that this includes the values used to specify the range.

When a for loop is used in an equation section, each iteration of the for loop generates a new equation for each equation inside the for loop. So in our case, we will generate equations corresponding to values of i between 2 and n-1.

Putting all this together, the complete model would be:

model Rod_ForLoop "Modeling heat conduction in a rod using a for loop"
  type Temperature=Real(unit="K", min=0);
  type ConvectionCoefficient=Real(unit="W/K", min=0);
  type ConductionCoefficient=Real(unit="W.m-1.K-1", min=0);
  type Mass=Real(unit="kg", min=0);
  type SpecificHeat=Real(unit="J/(K.kg)", min=0);
  type Density=Real(unit="kg/m3", min=0);
  type Area=Real(unit="m2");
  type Volume=Real(unit="m3");
  type Length=Real(unit="m", min=0);
  type Radius=Real(unit="m", min=0);

  constant Real pi = 3.14159;

  parameter Integer n=10;
  parameter Length L=1.0;
  parameter Radius R=0.1;
  parameter Density rho=2.0;
  parameter ConvectionCoefficient h=2.0;
  parameter ConductionCoefficient k=10;
  parameter SpecificHeat C=10.0;
  parameter Temperature Tamb=300 "Ambient temperature";

  parameter Area A_c = pi*R^2, A_s = 2*pi*R*L;
  parameter Volume V = A_c*L/n;

  Temperature T[n];
initial equation
  T = linspace(200,300,n);
equation
  rho*V*C*der(T[1]) = -h*A_s*(T[1]-Tamb)-k*A_c*(T[1]-T[2])/(L/n);
  for i in 2:(n-1) loop
    rho*V*C*der(T[i]) = -h*A_s*(T[i]-Tamb)-k*A_c*(T[i]-T[i-1])/(L/n)-k*A_c*(T[i]-T[i+1])/(L/n);
  end for;
  rho*V*C*der(T[end]) = -h*A_s*(T[end]-Tamb)-k*A_c*(T[end]-T[end-1])/(L/n);
end Rod_ForLoop;

Simulating this model yields the following solution for each of the nodal temperatures:

Note how the temperatures are initially distributed linearly (as we specified in our initial equation section).

Alternatives¶

It turns out that there are several ways we can generate the equations we need. Each has its own advantages and disadvantages depending on the context. We’ll present them here just to demonstrate the possibilities. The choice of which one they feel leads to the most understandable equations is up to the model developer.

One array feature we can use to make these equations slightly simpler is called an array comprehension. An array comprehension flips the for loop around so that we take a single equation and add some information at the end indicating that the equation should be evaluated for different values of the loop index variable. In our case, we can represent our equation section using array comprehensions as follows:

equation
  rho*V*C*der(T[1]) = -h*A_s*(T[1]-Tamb)-k*A_c*(T[1]-T[2])/(L/n);
  rho*V*C*der(T[2:n-1]) = {-h*A_s*(T[i]-Tamb)-k*A_c*(T[i]-T[i-1])/(L/n)-k*A_c*(T[i]-T[i+1])/(L/n) for i in 2:(n-1)};
  rho*V*C*der(T[end]) = -h*A_s*(T[end]-Tamb)-k*A_c*(T[end]-T[end-1])/(L/n);

We could also combine the array comprehension with some if expressions to nullify contributions to the heat balance that don’t necessarily apply. In that case, we can simplify the equation section to the point where it contains one (admittedly multi-line) equation:

equation
  rho*V*C*der(T) = {-h*A_s*(T[i]-Tamb)
                    -(if i==1 then 0 else k*A_c/(L/n)*(T[i]-T[i-1]))
                    -(if i==n then 0 else k*A_c/(L/n)*(T[i]-T[i+1])) for i in 1:n};

Recall, from several previous examples, that Modelica supports vector equations. In these cases, when the left hand and right hand side are vectors of the same size, we can use a single (vector) equation to represent many scalar equations. We can use this feature to simplify our equations as follows:

equation
  rho*V*C*der(T[1]) = -h*A_s*(T[1]-Tamb)-k*A_c*(T[1]-T[2])/(L/n);
  rho*V*C*der(T[2:n-1]) = -h*A_s*(T[i]-Tamb)-k*A_c*(T[2:n-1]-T[1:n-2])/(L/n)-k*A_c*(T[2:n-1]-T[3:n])/(L/n);
  rho*V*C*der(T[end]) = -h*A_s*(T[end]-Tamb)-k*A_c*(T[end]-T[end-1])/(L/n);

Note that when a vector variable like T has a range of subscripts applied to it, the result is a vector containing the components indicated by the values in the subscript. For example, the expression T[2:4] is equivalent to {T[2], T[3], T[4]}. The subscript expression doesn’t need to be a range. For example, T[{2,5,9}] is equivalent to {T[2], T[5], T[9]}.

Finally, let us consider one last way of refactoring these equations. Imagine we introduced three additional vector variables:

  Heat Qconv[n];
  Heat Qleft[n];
  Heat Qright[n];

Then we can write these two equations (again using vector equations) to define the heat lost to the ambient, previous section and next section in the rod:

  Qconv = {-h*A_s*(T[i]-Tamb) for i in 1:n};
  Qleft = {(if i==1 then 0 else -k*A_c*(T[i]-T[i-1])/(L/n)) for i in 1:n};
  Qright = {(if i==n then 0 else -k*A_c*(T[i]-T[i+1])/(L/n)) for i in 1:n};

This allows us to express the heat balance for each section using a vector equation that doesn’t include any subscripts:

  rho*V*C*der(T) = Qconv+Qleft+Qright;

Conclusion¶

In this section, we’ve seen various ways that we can use vector variables and vector equations to represent one-dimensional heat transfer. Of course, this vector related functionality can be used for a wide range of different problem types. The goal of this section was to introduce several features to demonstrate the various options that are available to the developer when working with vectors.