The coolest ODE function and plot in MATLAB
19 views (last 30 days)
Show older comments
What is the coolest ODE function and plot you came along?
0 Comments
Answers (1)
Davide Masiello
on 19 Apr 2022
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 Comment
Shishir Raut
on 25 Feb 2023
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.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!