A Mechanical Example

If you are more familiar with mechanical systems, this example might help reinforce some of the concepts we’ve covered so far. The system we wish to model is the one shown in the following figure:

Plant with pulse counter sensor

It is worth pointing out how much easier it is to convey the intention of a model by presenting it in schematic form. Assuming appropriate graphical representations are used, experts can very quickly understand the composition of the system and develop an intuition about how it is likely to behave. While we are currently focusing on equations and variables, we will eventually work our way up to an approach (in the upcoming section of the book on Components) where models will be built in this schematic form from the beginning.

For now, however, we will focus on how to express the equations associated with this simple mechanical system. Each inertia has a rotational position, \varphi , and a rotational speed, \omega where \omega = \dot{\varphi} . For each inertia, the balance of angular momentum for the inertia can be expressed as:

J \dot{\omega} = \sum_i \tau_i

In other words, the sum of the torques, \tau , applied to the inertia should be equal to the product of the moment of inertia, J , and the angular acceleration, \dot{\omega} .

At this point, all we are missing are the torque values, \tau_i . From the previous figure, we can see that there are two springs and two dampers. For the springs, we can use Hooke’s law to express the relationship between torque and angular displacement as follows:

\tau = c \Delta \varphi

For each damper, we express the relationship between torque and relative angular velocity as:

\tau = d \Delta \dot{\varphi}

If we pull together all of these relations, we get the following system of equations:

\begin{split}\omega_1 &= \dot{\varphi}_1 \\ J_1 \dot{\omega}_1 &= c_1 (\varphi_2-\varphi_1) + d_1 \frac{\mathrm{d} (\varphi_2-\varphi_1)}{\mathrm{d}t} \\ \omega_2 &= \dot{\varphi}_2 \\ J_2 \dot{\omega}_2 &= c_1 (\varphi_1-\varphi_2) + d_1 \frac{\mathrm{d} (\varphi_1-\varphi_2)}{\mathrm{d}t} - c_2 \varphi_2 - d_2 \dot{\varphi}_2\end{split}

Let’s assume our system has the following initial conditions as well:

\begin{split}\varphi_1 &= 0 \\ \omega_1 &= 0 \\ \varphi_2 &= 1 \\ \omega_2 &= 0\end{split}

These initial conditions essentially mean that the system starts in a state where neither inertia is actually moving (i.e., \omega=0 ), but there is a non-zero deflection across both springs.

Pulling all of these variables and equations together, we can express this problem in Modelica as follows:

model SecondOrderSystem "A second order rotational system"
  type Angle=Real(unit="rad");
  type AngularVelocity=Real(unit="rad/s");
  type Inertia=Real(unit="kg.m2");
  type Stiffness=Real(unit="N.m/rad");
  type Damping=Real(unit="N.m.s/rad");
  parameter Inertia J1=0.4 "Moment of inertia for inertia 1";
  parameter Inertia J2=1.0 "Moment of inertia for inertia 2";
  parameter Stiffness c1=11 "Spring constant for spring 1";
  parameter Stiffness c2=5 "Spring constant for spring 2";
  parameter Damping d1=0.2 "Damping for damper 1";
  parameter Damping d2=1.0 "Damping for damper 2";
  Angle phi1 "Angle for inertia 1";
  Angle phi2 "Angle for inertia 2";
  AngularVelocity omega1 "Velocity of inertia 1";
  AngularVelocity omega2 "Velocity of inertia 2";
initial equation
  phi1 = 0;
  phi2 = 1;
  omega1 = 0;
  omega2 = 0;
equation
  // Equations for inertia 1
  omega1 = der(phi1);
  J1*der(omega1) = c1*(phi2-phi1)+d1*der(phi2-phi1);
  // Equations for inertia 2
  omega2 = der(phi2);
  J2*der(omega2) = c1*(phi1-phi2)+d1*der(phi1-phi2)-c2*phi2-d2*der(phi2);
end SecondOrderSystem;

As we did with the low-pass filter example, RLC1, let’s walk through this line by line.

As usual, we start with the name of the model:

model SecondOrderSystem "A second order rotational system"

Next, we introduce physical types for a rotational mechanical system, namely:

  type Angle=Real(unit="rad");
  type AngularVelocity=Real(unit="rad/s");
  type Inertia=Real(unit="kg.m2");
  type Stiffness=Real(unit="N.m/rad");
  type Damping=Real(unit="N.m.s/rad");

Then we define the various parameters used to represent the different physical characteristics of our system:

  parameter Inertia J1=0.4 "Moment of inertia for inertia 1";
  parameter Inertia J2=1.0 "Moment of inertia for inertia 2";
  parameter Stiffness c1=11 "Spring constant for spring 1";
  parameter Stiffness c2=5 "Spring constant for spring 2";
  parameter Damping d1=0.2 "Damping for damper 1";
  parameter Damping d2=1.0 "Damping for damper 2";

For this system, there are four non-parameter variables. These are defined as follows:

  Angle phi1 "Angle for inertia 1";
  Angle phi2 "Angle for inertia 2";
  AngularVelocity omega1 "Velocity of inertia 1";
  AngularVelocity omega2 "Velocity of inertia 2";

The initial conditions (which we will revisit shortly) are then defined with:

initial equation
  phi1 = 0;
  phi2 = 1;
  omega1 = 0;
  omega2 = 0;

Then come the equations describing the dynamic response of our system:

equation
  // Equations for inertia 1
  omega1 = der(phi1);
  J1*der(omega1) = c1*(phi2-phi1)+d1*der(phi2-phi1);
  // Equations for inertia 2
  omega2 = der(phi2);
  J2*der(omega2) = c1*(phi1-phi2)+d1*der(phi1-phi2)-c2*phi2-d2*der(phi2);

And finally, we have the closing of our model definition.

end SecondOrderSystem;

The only drawback of this model is that all of our initial conditions have been “hard-coded” into the model. This means that we will be unable to specify any alternative set of initial conditions for this model. We can overcome this issue, as we did with our Newton cooling examples, by defining parameter variables to represent the initial conditions as follows:

model SecondOrderSystemInitParams
  "A second order rotational system with initialization parameters"
  type Angle=Real(unit="rad");
  type AngularVelocity=Real(unit="rad/s");
  type Inertia=Real(unit="kg.m2");
  type Stiffness=Real(unit="N.m/rad");
  type Damping=Real(unit="N.m.s/rad");
  parameter Angle phi1_init = 0;
  parameter Angle phi2_init = 1;
  parameter AngularVelocity omega1_init = 0;
  parameter AngularVelocity omega2_init = 0;
  parameter Inertia J1=0.4 "Moment of inertia for inertia 1";
  parameter Inertia J2=1.0 "Moment of inertia for inertia 2";
  parameter Stiffness c1=11 "Spring constant for spring 1";
  parameter Stiffness c2=5 "Spring constant for spring 2";
  parameter Damping d1=0.2 "Damping for damper 1";
  parameter Damping d2=1.0 "Damping for damper 2";
  Angle phi1 "Angle for inertia 1";
  Angle phi2 "Angle for inertia 2";
  AngularVelocity omega1 "Velocity of inertia 1";
  AngularVelocity omega2 "Velocity of inertia 2";
initial equation
  phi1 = phi1_init;
  phi2 = phi2_init;
  omega1 = omega1_init;
  omega2 = omega2_init;
equation
  omega1 = der(phi1);
  omega2 = der(phi2);
  J1*der(omega1) = c1*(phi2-phi1)+d1*der(phi2-phi1);
  J2*der(omega2) = c1*(phi1-phi2)+d1*der(phi1-phi2)-c2*phi2-d2*der(phi2);
end SecondOrderSystemInitParams;

In this way, the parameter values can be changed either in the simulation environment (where parameters are typically editable by the user) or, as we will see shortly, via so-called “modifications”.

You will see in this latest version of the model that the values for the newly introduced parameters are the same as the hard-coded values used earlier. As a result, the default initial conditions will be exactly the same as they were before. But now, we have the freedom to explore other initial conditions as well. For example, if we simulate the SecondOrderSystemInitParams model as is, we get the following solution for the angular positions and velocities:

/static/_images/SOSIP1.svg

However, if we modify the phi1_init parameter to be 1 at the start of our simulation, we get this solution instead:

/static/_images/SOSIP11.svg

Expanding on this mechanical example

If you would like to see this example further developed, you may wish to jump to the set of examples involving rotational systems found in the section on Speed Measurement.

Otherwise, you can continue to the next set of examples which involve population dynamics.