一维热传导

我们前面关于状态空间的讨论引入了矩阵和向量。上面讨论的焦点主要是数组的数学性质。在本节中,我们将考虑如何用数组来表述更实际的问题,即变量的一维空间分布。我们将观察Modelica与数组有关的语言特性,以及这些特性如何帮助我们更简略地表达系统行为。

我们的问题关于一个简单的热传导问题。考虑如下所示的一维棒:

Discretization of a one-dimensional bar

推导方程

让我们考虑此杆在每个离散截面的热平衡。首先,我们可以得到第\(i\)个截面的热容。这可以表示为:

\[m_i C T_i\]

其中\(m\)为第\(i\)节的质量,\(C\)为比热容(材料属性), 而\(T\)i则是第\(i\)节的温度,我们可以进一步将质量描述为:

\[m_i = \rho V_i\]

其中\(\rho\)是材料的密度,而\(V_i\)则是第\(i\)节的体积。最后,第\(i\)节的体积如下:

\[V_i = A_i L_i\]

其中\(A_i\)为第\(i\)节的截面积(假设其为恒定),而\(L_i\)为第\(i\)节的长度。对于这个例子,我们假定该杆是由相等大小的部分组成。在这种情况下,我们可以定义各段长度\(L_i\)如下:

\[L_i = \frac{L}{n}\]

我们也将假定全杆的截面积为恒定。因此,每节的质量可以写为:

\[m = \rho A L_i\]

在这种情况下,各部分的热容量会是:

\[\rho A L_i C T_i\]

反过来,这意味着在任何时间该节中的热量增长是:

\[\rho A L_i C \frac{\mathrm{d} T_i}{\mathrm{d}t}\]

其中我们假定\(A\)\(L_i\)还有\(C\)不会随着时间变化。

完成了以上热容。此外,我们将考虑两种不同形式的热传递。热传递的第一种形式,我们将考虑每节在特定环境温度\(T_{amb}\)下的对流。

\[q_h = -h A (T_i-T_{amb})\]

其中\(h\)为对流系数。另一个形式的热传递是相邻部分之间的热传导。这里有两个影响因素,其中之一是与\({i-1}\)元素的热传导(如果这个元素存在),另外则是与\({i-1}\)元素的热传导(前提也是该元素存在)。这些因素可以分别表示为:

\[q_{k_{i \rightarrow {i-1}}} = -k A \frac{T_i-T_{i-1}}{L_i}\]
\[q_{k_{i \rightarrow {i+1}}} = -k A \frac{T_i-T_{i+1}}{L_i}\]

利用上述关系,我们对首个元素的热平衡方程有:

\[\begin{split}\rho A L_i C \frac{\mathrm{d} T_1}{\mathrm{d}t} &= -k A \frac{T_1-T_2}{L_i} - h A (T_1-T_{amb})\end{split}\]

类似地,最后一个元素的热平衡方程为:

\[\begin{split}\rho A L_i C \frac{\mathrm{d} T_n}{\mathrm{d}t} &= -k A\frac{T_n-T_{n-1}}{L_i} -h A (T_n-T_{amb})\end{split}\]

最后,所有其他元素的热平衡有:

\[\begin{split}\rho A L_i C \frac{\mathrm{d} T_i}{\mathrm{d}t} &= -k A\frac{T_i-T_{i-1}}{L_i} -k A\frac{T_i-T_{i+1}}{L_i} -h A (T_i-T_{amb})\end{split}\]

实现

我们首先定义各物理量的类型。定义类型可以让变量有正确的单位。而在特定的软件会支持对方程的进行量纲检测。我们的类型定义如下:

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

我们也将定义几个参数以描述所模拟的棒。

  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 = pi*R^2;
  parameter Volume V = A*L/n;

最后,这个问题里唯一的数组是各部分的温度(因为实际上这是唯一沿杆长而改变的量):

  Temperature T[n];

以上就是我们需要做出的所有声明。现在让我们考虑所需的各种方程。首先,我们需要指定杆的初始条件。我们将假定\(T_1(0)=200\)\(T_n(0)=300\),而其他所有元素的初始温度则为上述条件的线性内插值。下列方程便体现了上述的条件:

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

其中linspace操作符用于产生一个n值的数组。而这n个值介于200300之间线性改变。回想之前的状态空间例子,我们可以加入其两边表达式均为向量的方程。本例的方程便是此类方程的又一个例子。

最后,我们介绍每部分的温度随时间变化的方程:

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

第一个方程对应的第\(1\)个元素的热平衡,而最后一个方程则对应第\(n\)个。中间的方程则对应了其他所有元素。注意这里使用了end作为下标。如果某维度的下标为一个表达式时,表达式内的end表示该维度的大小。在这里,我们用end表示最后一个元素。当然,我们在这里也可以用n来表示,但一般而言,当某个维度的大小不曾与某个特定变量相关时,end可能会非常有用。

另外,请注意模型内使用的for循环。for循环让循环变量在一定范围的值里变化。在本例里,循环变量为i,而值的变化范围则是从\(2\)\(n-1\)for循环的基本语法如下:

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

其中,<range>是一个含有不同值的向量。而一种方便的产生数列的方式是使用范围操作符。范围操作符前的值是数列的初值,而在操作符后的值则是数列的终值。因此,以表达式5:10例,表达式会生成一个由5678910组成的向量。注意,生成的向量包括了用于指定范围的两个值。

for循环用在等式区域时,循环内每个方程都会在for循环每次迭代时产生一个新的方程。因此,在本例对应每个在2n-1之间的i,我们会一共生成\(n-2\)个方程。

综合以上内容,完整的模型如下:

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 = pi*R^2;
  parameter Volume V = A*L/n;

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

Note

请注意,我们用常量形式在模型里包含了pi。在这本书的后面,我们将讨论如何正确导入常用的常数

对模型进行仿真就会得到各节点温度的解,如下:

../../../_images/RFL.png

注意温度在开始是如何线性分布的(正如我们在initial equation区域所指定的一样)。

其他可选方法

实际上,有好几种方法可以生成我们需要的方程。而每种方法在不同的情景下有其优点与缺点。我们将在下面一一介绍这些可能的方法。而选择哪种方法可以让方程感觉上最容易理解,则是取决于模型开发者自己。

我们可以用一个数组特性让这些方程变得更为简单。这个特性叫做数组解析(译注:原文为array comprehension,类似Python的list comprehension)。数组解析将for循环倒了过来。也就是说,我们在单一的等式后加上了循环变量在循环时如何取不同值的信息。在我们的例子中,我们可以使用数组解析将方程表达为如下方式:

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

我们还可以在数组解析内加上一些if表达式以删去不存在的热平衡效应。在这种情况下,我们可以把equation简化为一个(虽然占用了多行的)方程:

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

回顾前述几个例子,Modelica语言支持矢量方程。在这些情况下,只要方程左右均是大小同样的向量,我们就可以使用一个(向量)方程来表示多个标量方程。我们可以利用这个特性来简化方程式,结果如下:

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

注意,倘若大家使用一定范围的下标去存取向量变量,如T,那么结果就会得到这些下标对应元素所组成的向量。例如,表达式T[2:4]等价于{T[2], T[3], T[4]}。下标表达式并不需要是一个范围表达式。例如,表达式T[{2,5,9}]等价于{T[2], T[5], T[9]}

最后,让我们考虑重构这些方程的最后一个方法。想象一下,我们引入了另外两个向量变量:

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

然后我们就可以写出以下的两个方程(再次使用矢量方程),用以定义到前后两节分别的热损失。

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

这使我们可以为用矢量方程表达每个部分的热平衡,而无须包含任何下标。

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

结论

在本节中,我们看到了使用向量变量以及向量方程式来表示一维传热的多种方法。当然,这些向量的相关功能可用于类型广泛的不同问题。本节目的是为大家介绍几种特性,以演示开发者在使用向量时拥有的不同选项。