Why is my model not converging with ode45 solver ?

Question

ALOYS 2024-8-16，11:52

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2145719-why-is-my-model-not-converging-with-ode45-solver

评论： Sam Chak 2024-8-17，10:37

I am trying to simulate an electronical device that can be modeled by a mass-spring-damper system with an additional non-linear force. The equation at the equilibrium for the system is the following :

The goal of my MATLAB code is to solve this equation for $x$ and find the equilibrium. For this purpose, I use the `ode45` function like this (all coefficient are defined in my code but not shown here) :

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

The equation is good according to my teachers and several papers but the solver never converges and I can't see why. All coefficient are defined according to the dimensions of a capacitive micromachined ultrasound transducer. The solution of ode45 is the following :

This is obviously wrong because the dimensions of a cmut are bellow millimeter.

Can you see any mistakes in the way I use the ode45 solver ?

Full code :

% 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

References :

- Y. Wang, L. -M. He, Z. Li, W. Xu and J. Ren, "A computationally efficient nonlinear dynamic model for cMUT based on COMSOL and MATLAB/Simulink"

- T. Merrien, A. Boulmé and D. Certon, "Lumped-Parameter Equivalent Circuit Modeling of CMUT Array Elements"

I tried several solver from MATLAB and implemented my own Runge-Kutta algorithm based on the Wikipedia example. I also verified the coefficients according to my references.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2024-8-16，12:21

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2145719-why-is-my-model-not-converging-with-ode45-solver#answer_1499799

Hi @ALOYS

My guess is that if the coefficient

is a positive value, then the system will be unstable.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Sam Chak 2024-8-16，12:30

在 MATLAB Online 中打开

Hi @ALOYS

I changed the sign of the coefficient to test if my guess is correct.

% 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;

tf = 2e-9;

dt = 1e-6;

% tspan = linspace(ti, tf, 1/dt);

tspan = [ti, tf];

% Résolution par RK4

opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);

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

plot(t, x), grid on

xlabel('t'), ylabel('x(t)')

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

ALOYS 2024-8-17，9:47

编辑：ALOYS 2024-8-17，9:47

You are right ! The sign on the K/B term was wrong. Now my model gives accurate predictions according to my tests on the real device. This will be the first part of my internship report.

Sam Chak 2024-8-17，10:37

Hi @ALOYS, I'm glad it worked out. If you find the solution helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Your support is greatly appreciated!

Thank you, @William Rose.

请先登录，再进行评论。