Chemical System¶

In this section, we’ll consider a few different ways to describe the behavior of a chemical system. We’ll start by building a model without using the array functionality. Then, we’ll implement the same behavior using vectors. Finally, we’ll implement the same model again using enumerations.

In all of our examples, we’ll be building a model for the following system of reactions:

It should be noted that is simply an intermediate result of this reaction. The overall reaction can be expressed as:

Using the law of mass action we can transform these chemical equations into the following mathematical ones:

where , and are the reaction coefficients for the first, second and third reactions, respectively. These equations are derived by considering the change in each species due to each reaction involving that species. So, for example, since the first reaction transforms molecules of and into molecules of , we see the term in the balance equation for , which represents the reduction in the amount of as a result of that reaction. Each term in these balance equations is derived in a similar fashion.

Without Arrays¶

Let us start with an approach that doesn’t utilize arrays at all. In this case, we simply represent the concentrations , and by the variables cA, cB and cX as follows:

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;

With this approach, we create an equation for the balance of each species. If we simulate this model, we get the following results:

Using Arrays¶

Another way to approach modeling of the chemical system is to use vectors. With this approach, we associated the species , and with the indices , and , respectively. The concentrations are mapped to the vector variable C. We can also cast the reaction coefficients into a vector of reaction coefficients, k.

With this transformation, all the equations are then transformed into vector equations:

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;

The reaction equations are non-linear, so they cannot be transformed into a completely linear form. But we could simplify them further by using a matrix-vector product. In other words, the equations:

can be transformed into the following form:

which can then be represented in Modelica as:

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;

The drawback of this approach is that we have to constantly keep track of which index (e.g., 1, 2, or 3) corresponds to which species (e.g., , , or ).

Using Enumerations¶

To address this issue of having to map back and forth from numbers to names, our third approach will utilize the enumeration type in Modelica. An enumeration allows us to define a set of names which we can then use to define the subscripts associated with an array. We’ll define our enumeration as follows:

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

This defines a special type named Species that has exactly three possible values, A, B and X. We can then use this enumeration as a dimension in an array as follows:

Real C[Species];

Since the Species type has only three possible values, this means that the vector C has exactly three components. We can then refer to the individual components of C as C[Species.A], C[Species.B] and C[Species.X].

Because it is awkward to constantly prefix each species name with Species, we can define a few convenient constants as follows:

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

In this way, we can now refer to the concentration of species as C[A]. Pulling all of this together we can represent our chemical system using enumerations as:

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;

Conclusion¶

In this chapter, we showed how a set of chemical equations could be represented with and without arrays. We also demonstrated how the enumeration type can be used in conjunction with arrays to make the resulting equations more readable by replacing numeric indices with names. Furthermore, this section also demonstrated how the enumeration type can be used not only to index the array, but also to define one or more dimensions in the declaration.