机械示例

如果你对机械系统比较熟悉,这个例子可能有助于你理解前面所介绍的一些概念。我们期望建模的系统如下图所示:

Plant with pulse counter sensor

值得指出的是以示意图的形式表达模型的意图是多么容易。假设通过适当的图示表达,专家就可以很快的理解系统的组成并对系统可能具有的特性有一定的了解。虽然现在我们专注于方程和变量的定义,我们最终将转换工作方式(在本书后续章节组件将会讲解),在开始建模的时候就以系统示意图的形式来搭建.

现在,我们依然将重点放在如何用方程来表示这个简单的机械系统。每个与惯量相关的角位移\(\varphi\),及角速度\(\omega\)之间有这样的关系式:\(\omega = \dot{\varphi}\)。对于每个惯量,根据角动量守恒可以得到如下表达式:

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

换句话说,即施加在惯量元素上的扭矩\(\tau\)总和等于转动惯量\(J\)和角加速度\(\dot{\omega}\)的乘积。

在上述方程中,我们唯一缺少的是扭矩值\(\tau_i\)。从上面的示意图上我们可以看到,上述机械系统包含两个弹簧和两个阻尼器。对于弹簧,我们可以根据胡克定律来表示扭矩和角位移之间的关系,如下所示:

\[\tau = k \Delta \varphi\]

对每个阻尼器,我们可以用如下方程表示其扭矩和相对角速度之间的关系:

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

如果我们将所有的表达式放在一起,可以得到如下的系统方程:

\[\begin{split}\omega_1 &= \dot{\varphi}_1 \\ J_1 \dot{\omega}_1 &= k_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 &= k_1 (\varphi_1-\varphi_2) + d_1 \frac{\mathrm{d} (\varphi_1-\varphi_2)}{\mathrm{d}t} - k_2 \varphi_2 - d_2 \dot{\varphi}_2\end{split}\]

我们假设系统具有以下的初始条件:

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

这些初始化条件意味着系统的开始状态惯量元素没有转动(即\(\omega=0\)),但是在两个弹簧之间有一个非零的偏转。

将上述所有变量和方程放在一起,我们就可以用Modelica语言来描述这个问题,如下所示:

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 k1=11 "Spring constant for spring 1";
  parameter Stiffness k2=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) = k1*(phi2-phi1)+d1*der(phi2-phi1);
  // Equations for inertia 2
  omega2 = der(phi2);
  J2*der(omega2) = k1*(phi1-phi2)+d1*der(phi1-phi2)-k2*phi2-d2*der(phi2);
end SecondOrderSystem;

像我们在低通滤波器RLC1例子中讲解的那样,让我们一步步的来讲解上述模型。

像往常一样,我们先从模型的名称开始:

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 k1=11 "Spring constant for spring 1";
  parameter Stiffness k2=5 "Spring constant for spring 2";
  parameter Damping d1=0.2 "Damping for damper 1";
  parameter Damping d2=1.0 "Damping for damper 2";

对于这个系统,有四个非parameter变量。定义如下:

  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) = k1*(phi2-phi1)+d1*der(phi2-phi1);
  // Equations for inertia 2
  omega2 = der(phi2);
  J2*der(omega2) = k1*(phi1-phi2)+d1*der(phi1-phi2)-k2*phi2-d2*der(phi2);

最后,定义模型的结束。

end SecondOrderSystem;

这个模型的唯一缺点是我们所有的初始化条件已经被“硬编码”到模型中。这也意味着,我们将不能指定任何其他组调用该模型的初始条件。我们可以克服这个问题,就像在Newton cooling examples例子中通过定义parameter变量来表示初始条件,如下所示:

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 k1=11 "Spring constant for spring 1";
  parameter Stiffness k2=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) = k1*(phi2-phi1)+d1*der(phi2-phi1);
  J2*der(omega2) = k1*(phi1-phi2)+d1*der(phi1-phi2)-k2*phi2-d2*der(phi2);
end SecondOrderSystemInitParams;

通过这种方式,该参数值即可以在仿真环境中更改(参数值通常会被用户编辑)。另外,我们也可以通过所谓的修改语句(modifications)来改变参数。

在这个最新版本的模型中你会看到,新设置的参数和之前硬编码参数的数值是一样的。因此,默认的初始条件和之前的完全一样。但是现在,我们有充分的自由去探索其他初始化条件的方法。例如,我们仿真SecondOrderSystemInitParams模型,可以得到角位移及角速度解的轨迹,如下图:

../../../_images/SOSIP1.png

但是,如果将参数phi1_init的值修改为1,可以得到以下的仿真结果:

../../../_images/SOSIP11.png

机械示例的扩展

如果你想了解本实例更深的扩展,可以直接跳到包含更多旋转系统实例的速度的测量一节里。

否则,你可以继续下面涉及种群动力学方面的实例。