化学系统

在本节中,我们会考虑几种不同的方式去描述化学系统的行为。首先,我们会在不使用数组功能的前提下进行建模。然后,我们会用向量描述相同的行为。最后,我们将使用枚举再次实现相同的模型。

在我们所有的例子里,我们建立的模型均是基于以下的反应系统

\[\begin{split}A + B &\underset{1}{\rightarrow} X \\ A + B &\underset{2}{\leftarrow} X \\ X + B &\underset{3}{\rightarrow} R + S\end{split}\]

必须注意\(X\)不过是反应的中间结果。总反应可表示为:

\[A + 2B \rightarrow R + S\]

使用质量作用定律,我们可以将这些化学方程式转化为一下数学方程:

\[\begin{split}\frac{\mathrm{d}[A]}{\mathrm{d}t} &= -k_1 [A] [B] + k_2 [X] \\ \frac{\mathrm{d}[B]}{\mathrm{d}t} &= -k_1 [A] [B] + k_2 [X] -k_3 [B] [X] \\ \frac{\mathrm{d}[X]}{\mathrm{d}t} &= k_1 [A] [B] - k_2 [X] -k_3 [B] [X]\end{split}\]

其中\(k_1\)\(k_2\)\(k_3\)分别是第一、二、三个反应的反应系数。这些方程是通过在考虑每种物质的变化,以及涉及该物质的每个反应后所推导得到的。因此,例如第一个反应\(A + B \rightarrow X\)是把\(A\)分子和\(B\)分子转化为\(X\)分子。我们可以看到\(-k_1 [A] [B]\)这项出现在\(A\)的平衡方程里。这项表示了\(A\)应为该反应而减少的量。这些平衡方程的每一项均是以类似的方式推导出来的。

不使用数组

首先,让我们采取完全不使用数组的方法。在这里,我们直接把浓度\([A]\)\([B]\)以及\([X]\)用变量cAcBcX表示如下:

model Reactions_NoArrays "Modeling a chemical reaction without arrays"
  Real cA;
  Real cB;
  Real cX;
  parameter Real k1=0.1;
  parameter Real k2=0.1;
  parameter Real k3=10;
initial equation
  cA = 1;
  cB = 1;
  cX = 0;
equation
  der(cA) = -k1*cA*cB + k2*cX;
  der(cB) = -k1*cA*cB + k2*cX - k3*cB*cX;
  der(cX) = k1*cA*cB - k2*cX - k3*cB*cX;
end Reactions_NoArrays;

通过这种方法,我们建立了每种物质的平衡方程。对模型进行仿真后,我们会得到以下结果:

../../../_images/RNA.png

使用数组

另一种对化学系统进行建模的方法是使用向量。用这种方法,我们将化学物质\(A\)\(B\)\(X\)分别用索引\(1\)\(2\)\(3\)标记。浓度则被映射到向量变量C。我们同时也将反应系数储存在反应系数向量k里。

进行了上述变换后,所有的方程都转化为向量方程:

model Reactions_Array "Modeling a chemical reaction with arrays"
  Real C[3];
  parameter Real k[3] = {0.1, 0.1, 10};
initial equation
  C = {1, 1, 0};
equation
  der(C) = {-k[1]*C[1]*C[2] + k[2]*C[3],
            -k[1]*C[1]*C[2] + k[2]*C[3] - k[3]*C[2]*C[3],
            k[1]*C[1]*C[2] - k[2]*C[3] - k[3]*C[2]*C[3]};
end Reactions_Array;

反应方程是非线性的,所以这些方程不能被转换成完全线性的形式。但是,我们可以通过使用矩阵向量积进行进一步简化。换句话说,方程:

\[\begin{split}\frac{\mathrm{d}[A]}{\mathrm{d}t} &= -k_1 [A] [B] + k_2 [X] \\ \frac{\mathrm{d}[B]}{\mathrm{d}t} &= -k_1 [A] [B] + k_2 [X] -k_3 [B] [X] \\ \frac{\mathrm{d}[X]}{\mathrm{d}t} &= k_1 [A] [B] - k_2 [X] -k_3 [B] [X]\end{split}\]

可以转化为下面的形式:

\[\begin{split}\frac{\mathrm{d}}{\mathrm{d}t} \left\{ \begin{array}{c} [A] \\[0pt] [B] \\[0pt] [X] \end{array} \right\} = \left[ \begin{array}{rrr} -k_1 [B] & 0 & k_2 \\ -k_1 [B] & -k_3 [X] & k_2 \\ k_1 [B] & -k_3 [X] & -k_2 \end{array} \right] \left\{ \begin{array}{c} [A] \\[0pt] [B] \\[0pt] [X] \end{array} \right\}\end{split}\]

而上述方程可以用以下的Modelica形式进行表示:

der(C) = [-k[1]*C[2], 0,          k[2];
          -k[1]*C[2], -k[3]*C[3], k[2];
          k[1]*C[2],  -k[3]*C[3], -k[2]]*C;

这种方法的问题是,我们必须时刻留意各个索引(如123)分别对应哪种物质(如\(A\)\(B\)\(X\))。

使用枚举

为了解决数字和名称来回映射的问题,我们的第三种方法则是利用Modelica的enumeration类型。枚举类型可以让我们定义一个名称的集合。然后,我们可以用这个名称的几何对应数组的下表。我们将定义以下的枚举:

type Species = enumeration(A, B, X);

上述语句定义了一个特殊的类型叫做Species,而这个类型有三个可能值:ABX。然后,我们可以使用这个枚举作为数组的一个维度,如下:

Real C[Species];

由于Species类型只有三种可能的值,这意味着矢量C恰好有三个分量。然后,我们可以分别使用C[Species.A]C[Species.B]C[Species.X]去指代C的分量。

由于每次在物质名称前加入Species颇为不便,为方便起见,我们可以定义以下常量:

constant Species A = Species.A;
constant Species B = Species.B;
constant Species X = Species.X;

这样一来,我们现在可以用C[A]来指代物质\(A\)的浓度。综合以上的结果,我们可以使用如下方式用枚举描述我们的化学系统:

model Reactions_Enum "Modeling a chemical reaction with enums"
  type Species = enumeration(
      A,
      B,
      X);
  Real C[Species] "Species concentrations";
  parameter Real k[3] = {0.1, 0.1, 10};
  constant Species A = Species.A;
  constant Species B = Species.B;
  constant Species X = Species.X;
initial equation
  C[A] = 1.0;
  C[B] = 1.0;
  C[X] = 0.0;
equation
  der(C[A]) = -k[1]*C[A]*C[B] + k[2]*C[X];
  der(C[B]) = -k[1]*C[A]*C[B] + k[2]*C[X] - k[3]*C[B]*C[X];
  der(C[X]) = k[1]*C[A]*C[B] - k[2]*C[X] - k[3]*C[B]*C[X];
end Reactions_Enum;
../../../_images/RE.png

结论

在这一章中,我们介绍了如何在使用以及不使用数组的前提下描述一组化学方程式。我们还演示了如何在数组中使用enumeration类型。这样用名称代替数字索引,可以让产生的方程式更具可读性。此外,本节还表明了enumeration类型可以不仅用于索引数组,也能在变量声明中定义一个或多个维度。