Rayleigh–Plesset equation (problem with ode45)

조회 수: 17 (최근 30일)
Md Ishak Khan
Md Ishak Khan 2016년 11월 8일
편집: Sam Chak 2025년 11월 21일
Hello, I'm trying to solve Rayleigh–Plesset equation to compute bubble radius at certain time instants.
The image above shows the equation. I've chosen fluid density, kg/m³, Pascal. I have to vary the Pressure at infinity () as a function of time. As the paper shows (from which I have taken the values of P at infinity and time range), the bubble radius must increase first and then it should collapse. But I have been getting the radius increment out of bound. What would be the solution to this problem? Am I missing something obvious?
The pulse code (varying P at infinity) is here:
function z = pulse(t)
z = (-.0002*t^3) + (0.0851*t^2) - (8.7073*t) + 52.853; % code
end
The differential equation code is here:
function f = RP(t, y)
p_ref = 3200; % reference pressure (Pa)
rho = 1e3; % density of the surrounding fluid (kg/m3)
f = [y(2);
(p_ref - pulse(t))/(rho*y(1)) - 1.5*(y(2)^2)/y(1)];
end
I'm taking the initial radius as 50 micrometer. The plotting code is here:
[t, y] = ode45(@RP, [0, 0.3], [50e-6, 0]);
plot(t, y(:,1)), grid
xlabel('time'),
ylabel('$R, \;(\rm{m})$', 'interpreter', 'latex')
title('Radius of the spherical bubble')
for i = 1:length(t)
z(i) = pulse(t(i));
end
plot(t, z), grid
xlabel('time'),
ylabel('$P_{\infty}, \;(\rm{Pa})$', 'interpreter', 'latex')
title('External pressure infinitely far from the bubble')
Any help with the code will be appreciated. Thanks.

답변 (2개)

Sulekha
Sulekha 2018년 2월 16일
Were you able to find a solution to this issue? I am also facing the same problem.

Jhonath wensenber
Jhonath wensenber 2025년 11월 21일
편집: Jhonath wensenber 2025년 11월 21일

Lo que sucede es que el término de (PR-PL)/rho debe tener unidades de m2/s2 para ello modifica la expresión a: g(PR-PL)/(g*rho), de tal modo que, [g]= m/s2, [PR-PL]=N/m2 y [g*rho]=9810 N/m3.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by