为什么我的模型与 Matlab 中的 ode45 求解器不收敛？

Question

我正在尝试模拟一种电子设备，该设备可以通过具有附加非线性力的质量弹簧阻尼系统进行建模。系统平衡时的方程如下：

rac{ psilon_0 a V^2}{2 x^2} - Bx' + Kx = 0

我的 MATLAB 代码的目标是求解 x 的方程并找到平衡点。为此，我使用像这样的

ode45

函数（所有系数都在我的代码中定义，但此处未显示）：

x0 = 0;
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);

function dxdt = odefun(t, x, eps, a, V, B, d0, K)
    dxdt = ((eps * a * V^2) ./ (2 * B * (d0 - x).^2)) + ((K / B) .* x);
end

根据我的老师和几篇论文，这个方程很好，但求解器永远不会收敛，我不明白为什么。所有系数都是根据电容式微机械超声换能器的尺寸定义的。 ode45 的解显示膜的无限位移。这显然是错误的，因为 cmut 的尺寸低于毫米。 你能看出我使用 ode45 求解器的方式有什么错误吗？

完整代码：

% DIMENSIONS DE LA MEMBRANE
e = 500e-9;                 % epaisseur [m]
r = 20e-6;                  % rayon [m]
d0 = 550e-9;                % epaisseur de cavite [m]
a = pi * r^2;               % surface de la membrane [m^2]

% PARAMETRES MECANIQUES
E = 200e9;                  % module d'Young du SiN [Pa]
nu = 0.25;                  % coefficient de poisson du SiN
eta = 18.5e-6;              % viscosité dynamique de l'air [Pa.s]
p0 = 1e5;                   % pression exterieure [Pa]
rhoSiN = 3170;              % densite du SiN [Kg/m^3]
m = rhoSiN * a * e;         % masse membrane [Kg]
K = (16 * E * e^3) / (3 * (1 - nu^2) * r^2);   % raideur [N/m]
B = (eta * pi * r^2) / e;                      % amortissement [N.s/m]

% PARAMETRES ELECTRIQUES
eps0 = 8.85e-12;            % permittivité diélectrique du vide [F/m]
V0 = 10;                    % tension de polarisation [V]

% Conditions initiales et plage de temps
x0 = 0;
ti = 0; tf = 1e-3;
dt = 1e-6;
tspan = linspace(ti, tf, 1/dt);

% Résolution par RK4
[t, x] = ode45(@(t, x) odefun(t, x, eps0, a, V0, B, d0, K), tspan, x0);

plot(t,x,'-')

function dxdt = odefun(t, x, eps, a, V, B, d0, K)
    dxdt = ((eps * a * V^2) ./ (2 * B * (d0 - x).^2)) + ((K / B) .* x);
end

参考资料：

• Y。王L.-M。 He、Z. Li、W. Xu 和 J. Ren，“基于 COMSOL 和 MATLAB/Simulink 的 cMUT 计算高效非线性动态模型”
• T。 Merrien、A. Boulmé 和 D. Certon，“CMUT 阵列元件的集总参数等效电路建模”

我尝试了 MATLAB 中的几种求解器，并根据维基百科示例实现了我自己的 Runge-Kutta 算法。我还根据我的参考资料验证了系数。最后我尝试使用 ode15s 刚性方程求解器和附加选项并得到了相同的结果。

Answer 1

此方程似乎在

x = d0

处不稳定，您需要确保在模拟期间不会经历这种不稳定。

此模拟从

x = 0

开始并向上，当

x = 550e-9 = d0

且解发散时，它会达到不稳定状态。

我不确定正确的分辨率，也许你想开始模拟超越不稳定（例如

x0=1e-8

）或者方程可能有问题。 （你好像在非线性项中从 x 中扣除了

d0

而不是刚度项？）