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

What are the classes ComplexU, ComplexPressure, Compliance and Inertance? I tried running your model but these seem to be part of another library you are using. I replaced them with MSL or primitive types.

Furthermore i don't really understand how that model is supposed to work, you only defined an initial equation block and no actual equations. I tried following model:

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;

I this what you wanted?

You can compare w with real_w. One is computed by solving the system and one is computed by the equation.

As you can see the standard solver struggles, but the total pivot solver manages to solve the system. It converges to the other side (p_d.re = -1e+2;) so maybe this is the correct value anyways?

EDIT: I changed the model to the correct one, i fiddled around with the initial equation, now everything works fine! The main solver still fails, but total pivot finds the solution.