如何找到振荡器的共振频率?
How to find a resonance frequency of an oscillator?
我目前正在尝试使用 OpenModelica 模拟声学谐振器,我想知道如何robustly/nicely 计算它们的谐振频率。
作为一个简化的例子(没有媒体等),我已经实现了一个双亥姆霍兹谐振器,本质上是两个由管道(惯性)连接的体积(顺应性)。实际系统由更多连接在一起的组件组成。
压力和体积流量(均为复值)的振荡遵循正弦表达式,具有共振 angular 频率 w
。这产生了 4 个压力和 4 个体积流量的 8 个方程(在端点和柔量-惯性连接点)。
这是我每晚在 OpenModelica 中求解的 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");
在纸上计算,系统的共振 angular 频率为 w=\sqrt{\frac{2}{LC}}
(在本例中为 ~1871 1/s),因此系统具有非零解。
为了避免解算器进入无趣的零解,我必须在一点添加一些刺激,因此初始方程 p_a.re = 1e+2
。
现在模拟这个,因为w
是一个附加变量,我需要引入一个附加方程,选择der(w) = 0;
作为共振在这种情况下频率是恒定的。
不幸的是,这使得无法达到更多 complex/realistic 的情况,其中共振频率随时间变化,例如随着温度或其他变化的值。
Q1: 有没有更好的方法来提供共振的附加方程 frequency/calculate 这个系统的特征值?
此外,模拟的成功取决于初始刺激的值(在某些范围内这会失败,或者我在每个时间步都得到奇异方程)。此外,实际上问题是在初始化阶段解决的。在最好的情况下,我得到输出
Simulation finishes, values reasonable, only during initialisation we get:
Matrix singular!
under-determined linear system not solvable!
The initialization finished successfully without homotopy method.
问题 2:有没有办法避免奇点 and/or 干净地处理这个初始化(例如 homotopy
)?
虽然这在简化的示例中足够有效(并导致 w
的正确值),但我担心对于更多 complex/realistics 模型,我可能会遇到更多有问题的数字困难。
我查看了 homotopy
,但我真的看不出如何在此处应用它。我想以某种方式将它应用于 w
,但 Fritzson 的书甚至似乎明确警告不要在导数表达式上使用它,除此之外,似乎只有 w.start
值会出现。
什么是类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;
) 所以也许这是正确的值?
编辑:
我将模型更改为正确的模型,我摆弄了初始方程式,现在一切正常!主求解器仍然失败,但 total pivot 找到了解决方案。
我想补充一件事关于失败的非线性求解器!
如果您正在处理最新的夜间构建,您应该会收到以下消息:
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))
。尽管您知道结果会很小,但它必须非常大才能避免奇点。
我目前正在尝试使用 OpenModelica 模拟声学谐振器,我想知道如何robustly/nicely 计算它们的谐振频率。
作为一个简化的例子(没有媒体等),我已经实现了一个双亥姆霍兹谐振器,本质上是两个由管道(惯性)连接的体积(顺应性)。实际系统由更多连接在一起的组件组成。
压力和体积流量(均为复值)的振荡遵循正弦表达式,具有共振 angular 频率 w
。这产生了 4 个压力和 4 个体积流量的 8 个方程(在端点和柔量-惯性连接点)。
这是我每晚在 OpenModelica 中求解的 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");
在纸上计算,系统的共振 angular 频率为 w=\sqrt{\frac{2}{LC}}
(在本例中为 ~1871 1/s),因此系统具有非零解。
为了避免解算器进入无趣的零解,我必须在一点添加一些刺激,因此初始方程 p_a.re = 1e+2
。
现在模拟这个,因为w
是一个附加变量,我需要引入一个附加方程,选择der(w) = 0;
作为共振在这种情况下频率是恒定的。
不幸的是,这使得无法达到更多 complex/realistic 的情况,其中共振频率随时间变化,例如随着温度或其他变化的值。
Q1: 有没有更好的方法来提供共振的附加方程 frequency/calculate 这个系统的特征值?
此外,模拟的成功取决于初始刺激的值(在某些范围内这会失败,或者我在每个时间步都得到奇异方程)。此外,实际上问题是在初始化阶段解决的。在最好的情况下,我得到输出
Simulation finishes, values reasonable, only during initialisation we get:
Matrix singular!
under-determined linear system not solvable!
The initialization finished successfully without homotopy method.
问题 2:有没有办法避免奇点 and/or 干净地处理这个初始化(例如 homotopy
)?
虽然这在简化的示例中足够有效(并导致 w
的正确值),但我担心对于更多 complex/realistics 模型,我可能会遇到更多有问题的数字困难。
我查看了 homotopy
,但我真的看不出如何在此处应用它。我想以某种方式将它应用于 w
,但 Fritzson 的书甚至似乎明确警告不要在导数表达式上使用它,除此之外,似乎只有 w.start
值会出现。
什么是类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;
) 所以也许这是正确的值?
编辑: 我将模型更改为正确的模型,我摆弄了初始方程式,现在一切正常!主求解器仍然失败,但 total pivot 找到了解决方案。
我想补充一件事关于失败的非线性求解器! 如果您正在处理最新的夜间构建,您应该会收到以下消息:
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))
。尽管您知道结果会很小,但它必须非常大才能避免奇点。