The coolest ODE function and plot in MATLAB

19 次查看(过去 30 天)
What is the coolest ODE function and plot you came along?

回答(1 个)

Davide Masiello
Davide Masiello 2022-4-19
For me, I think that would be any variation of the Rayleigh-Plesset equation, which describes the volume oscillations of a bubble in a liquid under the effect of a pressure field.
Below, see an example of the Keller-Miksis equation (similar to Rayleigh-Plesset, but accounting for mild liquid compressibility) simulating the radial oscillation of an inertially collpasing bubble.
clear,clc
T0 = 298 ; % Ambient temperature, K
freq = 26500; % Driving frequency, Hz
r0 = 4.5 ; % Initial bubble radius, microns
p0 = 1 ; % Ambient pressure, atm
pA = 1.2 ; % Acoustic pressure, atm
rho0 = 997 ; % liquid density, kg m-3
sigmaL = 0.072 ; % surface tension, N m-1
muL = 1e-3 ; % liquid viscosity, kg m-1 s-1
c = 1485 ; % speed of sound, m s-1
gamma = 1.4 ; % air specific heats ratio, -
tspan = [0 1/freq] ;
IC = [r0*1E-6,0,p0*101325+2*sigmaL/(r0*1E-6),T0] ;
odefun = @(t,w)keller_miksis(t,w,r0*1E-6,T0,rho0,freq,p0*101325,pA*101325,sigmaL,muL,c,gamma) ;
[time,SOL] = ode15s(odefun,tspan,IC) ;
t = time*freq ;
R = SOL(:,1)*1E+6 ;
S = SOL(:,2) ;
p = SOL(:,3)/101325 ;
T = SOL(:,4) ;
set(0,'defaulttextinterpreter','latex')
set(0,'defaultAxesTickLabelInterpreter','latex')
set(0,'defaultLegendInterpreter','latex')
figure(1)
subplot(2,2,1)
plot(t,R)
title('Bubble radius')
xlabel('$\omega t$')
ylabel('$R$, $\mu m$')
subplot(2,2,2)
semilogy(t,p)
title('Gas pressure')
xlabel('$\omega t$')
ylabel('$p$, atm')
subplot(2,2,3)
plot(t,T(:,1))
title('Gas temperature')
xlabel('$\omega t$')
ylabel('$T$, K')
subplot(2,2,4)
plot(t,abs(S/c))
title('Interface Mach number')
xlabel('$\omega t$')
ylabel('$Ma$, -')
function f = keller_miksis(t,w,r0,T0,rho0,freq,p0,pA,sigmaL,muL,c,gamma)
R = w(1) ;
dR = w(2) ;
p = w(3) ;
T = w(4) ;
dpdt = -3*gamma*p*dR/R ;
dTdt = 3*(1-gamma)*T*dR/R ;
pinf = p0-pA*sin(2*pi*freq*t) ;
pL = (p0+2*sigmaL/r0)*(r0/R)^(3*gamma)-2*sigmaL/R-4*muL*dR/R ;
ddR = ((1+dR/c)*(pL-pinf)/(rho0)-1.5*(1-dR/(3*c))*(dR^2)+...
(R/(rho0*c))*(dpdt+(2*sigmaL*dR+4*muL*dR^2)/(R^2)))/((1-dR/c)*R+4*muL/(rho0*c)) ;
f = [ dR ;...
ddR ;...
dpdt ;...
dTdt ];
end
  1 个评论
Shishir Raut
Shishir Raut 2023-2-25
Hello, I had a query with regards to the way you obtained the temperature profile. I assume the pressure profile was obtained considering an adiabatic condition that [PR^(3*gamma) = constant] but what is the equation for temperature based on? I would appreciate your help on this matter.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by