How to find a resonance frequency of an oscillator?

Question

I am currently trying to simulate acoustic resonators using OpenModelica, and I am wondering how to robustly/nicely calculate their resonance frequency.

As a simplified example (without Media, etc), I have implemented a dual Helmholtz resonator, in essence two volumes (compliances) connected by a pipe (inertance). The real system consists of more components connected together. The oscillation of pressure and volume flow (both complex values) follow a sinusoidal expression, with a resonance angular frequency w. This yields 8 equations for 4 pressures and 4 volume flows (at the end points and compliance-inertance connection points).

Here a Modelica code I am solving in OpenModelica nightly:

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;

with the additional definitions

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");

Calculating on paper, one arrives at a resonance angular frequency of w=\sqrt{\frac{2}{LC}} (in this case ~1871 1/s) for the system to have a non-zero solution.

To avoid the solver going to the uninteresting solution of zero everything, I have to add some stimulation at one point, hence the initial equation p_a.re = 1e+2.

Now, to simulate this, since the w is an additional variable, I need to introduce an additional equation, choosing der(w) = 0; as the resonance frequency is constant in this case. Unfortunately, this makes it impossible to go to a more complex/realistic situation, where the resonance frequency changes over time e.g. with temperature or other changing values.

Q1: Is there a better way to supply the additional equation for the resonance frequency/calculate this eigenvalue of the system?

In addition, the success of simulation depends on the value of the initial stimulation (in some ranges this fails or I get singular equations at every time step). Also, in effect the problem is being solved in the initialisation phase. In the best case, I get the output

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: Is there a way to avoid the singularity and/or deal with this initialisation cleanly (e.g. with homotopy)? While this works adequately in the simplified example (and results in correct values for w), I am concerned that for more complex/realistics models, I could encountered more problematic numeric difficulties. I have looked into homotopy, but I couldn't really see how to apply this here. I thought applying this to w somehow, but the Fritzson book even seems to warn explicitly against using this on the derivative expression, and aside of that only the w.start value seems to present itself.

Answer 1

ComplexU、ComplexPressure和Compliance是什么类Inertance？我尝试运行您的模型，但这些似乎是您正在使用的另一个库的一部分。我用 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;

我这就是你想要的吗？

您可以w与进行比较real_w。一种是通过求解系统来计算，一种是通过方程来计算。

正如您所看到的，标准求解器很困难，但总枢轴求解器设法解决了系统问题。它收敛到另一边 ( p_d.re = -1e+2;) 所以也许这就是正确的值？

编辑：我将模型更改为正确的模型，我摆弄了初始方程，现在一切正常！主求解器仍然失败，但总枢轴找到了解决方案。