如何找到振荡器的共振频率?

8

我目前正在尝试使用OpenModelica模拟声学共振器,并想知道如何稳健/优雅地计算它们的共振频率。

作为一个简化示例(不考虑介质等因素),我实现了一个双Helmholtz共振器,基本上是由两个容积(顺应度)通过一条管道(惯性)连接而成。真实系统由更多组件连接在一起。压力和体积流的振荡(均为复值)遵循正弦表达式,具有共振角频率w。这会产生8个方程式,用于计算4个末端和4个顺应度-惯性连接点处的压力和体积流。

以下是我正在OpenModelica Nightly中解决的Modelica代码:

model Helmholtz_test "Simple test of double Helmholtz resonator"
  constant Complex j = Modelica.ComplexMath.j;
  ComplexU U_a, U_b, U_c, U_d "Oscillating volume flow rate";
  ComplexPressure p_a, p_b, p_c, p_d "Oscillating pressure";
  Modelica.SIunits.AngularFrequency w(start=2000, fixed=false);
  Compliance C = 7.14e-9;
  Inertance L = 80;
initial equation
  p_a.re = 1e+2; //Simulation finishes, values reasonable, only during initialisation we get:
  //Matrix singular!
  //under-determined linear system not solvable!
  //The initialization finished successfully without homotopy method.
equation
//BCs on ends
  U_a = Complex(0);
  U_d = Complex(0);
//Left compliance a-b;
  p_a = p_b;
  p_a = -1 / (j * w * C) * (U_b - U_a);
//Inertance b-c
  U_b = U_c;
  p_c - p_b = -j * w * L * U_b;
//Right compliance c-d
  p_c = p_d;
  p_c = -1 / (j * w * C) * (U_d - U_c);
//Additional condition for Eigenvalue
  der(w) = 0;
//w^2 = 2/(L*C); //The "real" resonance frequency
  annotation(
    experiment(StartTime = 0, StopTime = 1, Tolerance = 1e-06, Interval = 0.002));
end Helmholtz_test;

加入额外的定义

operator record ComplexPressure =
  Complex(redeclare Modelica.SIunits.Pressure re,
           redeclare Modelica.SIunits.Pressure im)
  "Complex pressure";

operator record ComplexU = 
  Complex(redeclare Modelica.SIunits.VolumeFlowRate re,
           redeclare Modelica.SIunits.VolumeFlowRate im)
  "Complex volume flow rate";

type Compliance = Real(final quantity = "Compliance", final unit = "m3.Pa-1");

type Inertance = Real(final quantity="Inertance", final unit="kg.m-4");

在纸上计算后,我们可以得到一个共振角频率的公式w=\sqrt{\frac{2}{LC}}(在这个例子中约为1871 1/s),以获得非零解。

为了避免求解器陷入无聊的零解,我必须在某个位置添加一些刺激,因此最初的方程式为p_a.re=1e+2

现在为了模拟这个情况,由于w是一个额外的变量,我需要引入一个额外的方程,并选择der(w)=0;,因为在这种情况下共振频率是恒定的。不幸的是,这使得不可能转向更复杂/更现实的情况,例如随温度或其他变化值而改变的共振频率。

Q1:有没有更好的方法来提供共振频率的额外方程/计算系统的特征值?

此外,仿真的成功取决于初始刺激的值(在某些范围内会失败,或者我在每个时间步骤中得到奇异方程)。事实上,在初始化阶段解决了这个问题。在最好的情况下,我可以得到输出:

Simulation finishes, values reasonable, only during initialisation we get:
Matrix singular!
under-determined linear system not solvable!
The initialization finished successfully without homotopy method.

Q2: 有没有办法避免奇异点和/或清洁地处理这种初始化(例如,使用homotopy)? 虽然这在简化示例中可以正常工作(并且会产生正确的w值),但我担心对于更复杂/现实模型,我可能会遇到更多问题的数值困难。 我已经研究了homotopy,但我真的不知道如何在这里应用它。我想以某种方式将其应用于w,但是Fritzson的书甚至似乎明确警告不要在导数表达式上使用它,并且除此之外似乎只有w.start值存在。

2个回答

1
这些类 ComplexU, ComplexPressure, ComplianceInertance 是什么?我尝试运行您的模型,但似乎它们是另一个库的一部分。我用 MSL 或原始类型替换了它们。
此外,我并不真正理解该模型应该如何工作,您只定义了一个 initial equation 块,没有实际方程式。我尝试遵循以下模型:
model Helmholtz_test "Simple test of double Helmholtz resonator"
  constant Complex j = Modelica.ComplexMath.j;
  Complex U_a, U_b, U_c, U_d "Oscillating volume flow rate";
  Complex p_a, p_b, p_c, p_d "Oscillating pressure";
  parameter Modelica.SIunits.AngularFrequency w(start=2000, fixed=false);
  Modelica.SIunits.AngularFrequency real_w; //The "real" resonance frequency
  Real C = 7.14e-9;
  Real L = 80;
initial equation
  p_a.re = 1e+2;
equation
  U_a = Complex(0);
  U_d = Complex(0);
  p_a = p_b;
  p_a = -1 / (j * w * C) * (U_b - U_a);
  U_b = U_c;
  p_c - p_b = -j * w * L * U_b;
  p_c = p_d;
  p_c = -1 / (j * w * C) * (U_d - U_c);
  real_w = abs(sqrt(2/(L*C))); //The "real" resonance frequency
end Helmholtz_test;

我这样做是否符合您的要求?

您可以将wreal_w进行比较。一个是通过解决系统计算得出的,另一个是通过方程计算得出的。

正如您所看到的,标准求解器很难解决问题,但总枢轴求解器成功地解决了该系统。它收敛到另一侧(p_d.re = -1e+2;),因此也许这是正确的值?

编辑: 我更正了模型,调整了初始方程,现在一切正常!主求解器仍然失败,但总枢轴找到了解决方案。


谢谢Karim!未知的东西分别是记录和类型,它们有助于使整个事物具有单位一致性。我稍后会更新问题。此外,在复制粘贴过程中equation丢失了,你已经正确猜到了适当的位置。 - Christoph
使用OMEdit,我如何查看此切换到总枢轴求解器?在模拟成功结束之前,我只看到上面显示的奇异矩阵提示。 - Christoph
很遗憾,我们必须更新一些警告和错误信息。目前默认的信息报告不正确。有时您看到的消息表示牛顿求解器没有收敛,因此牛顿求解器的初始值会变化(如1%)。之后,如果牛顿有时收敛,如果没有,则使用备用的全量主元求解器。在这种情况下,改变初始值就足够了。您可以通过模拟->模拟设置->模拟标志并选中LOG_NLS_V框来获取更多信息。 - kabdelhak
我错误地假设它是完全的枢轴,因为这通常是情况。 - kabdelhak
感谢你的帮助排除问题!你是否了解问题中的Q1,即是否有更好的方法提供缺失的方程来获取特征值,而不是通过“parameter”或“der(w)=0”使“w”保持不变呢? - Christoph
没问题!由于w是一个参数,而你所做的有点像参数识别,将其声明为参数应该是最好和最自然的方法。 - kabdelhak

0
我想再补充一点关于失败的非线性求解器!如果您正在使用最新的夜间构建版本,则应该会收到以下消息:
Nonlinear iteration variables with default zero start attribute in NLSJac8. (1)
========================================
1: U_b.im:VARIABLE()  "Imaginary part of complex number" type: Real 
Nonlinear iteration variables with predefined start attribute in NLSJac8. (1)
========================================
1: w:VARIABLE(start = 2000.0 unit = "rad/s" fixed = false )  type: Real Info: Only non-linear iteration variables in non-linear eqation systems require start values. All other start values have no influence on convergence and are ignored. Use "-d=dumpLoops" to show all loops. In OMEdit Tools->Options->Simulation->OMCFlags, in OMNotebook call setCommandLineOptions("-d=dumpLoops") 

这里报告的每个变量都应该有一个起始值,如您所见w已经有了一个起始值,但是U_b的虚部缺少一个。在声明它时,您可以将其更改为U_b(im(start=10))。尽管您知道结果会相当小,但它必须相当大才能避免奇点。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接