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.

Meader

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.

Carolinecarolingian

Thank you Karim! The unknown things are records and types, respectively, that help making the whole thing unit-consistent. I'll update the question later.Also, the equation got lost during copy-paste, you already guessed the proper place correctly. – Meader

Using OMEdit, how do I see this switch to the total pivot solver? I only see the hint about the singular matrix shown above before the sim concludes successfully. – Meader

Unfortunately we have to update some of the warnings and error messages. Right now the default ones do not report correctly. The message you see sometimes indicates that the newton solver did not converge and therefore the initial values for the newton solver are varied (by a small amount like 1%). Afterwards the newton sometimes converges, if it does not the fallback total pivot solver is used. In this case varying the initial values was sufficient. You can have more information on that using Simulation -> Simulation Setup -> Simulation Flags and check the box for LOG_NLS_V. – Carolinecarolingian

I incorrectly assumed it was total pivot, as it is oftentimes the case. – Carolinecarolingian

Thanks for your help in troubleshooting! Do you known about Q1 in the question, i.e. is there a better way to provide the missing equation for getting the eigenvalue than making w constant (either by parameter or by der(w)=0)? – Meader

No problem! Since w is a parameter and what you are doing kind of is a parameter identification, declaring it as a parameter should be the best and most natural way to do it. – Carolinecarolingian

Answer 2

I want to add one more thing about the failing nonlinear solver! If you are working on the newest nightly build you should get following message:

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

Every variable reported here should have a start value, as you can see w already has a start value, but the imaginary part of U_b is missing one. When declaring it you could change it to U_b(im(start=10)). Although you know the result will be rather small, it has to be quite big to avoid singularities.

Carolinecarolingian