The coolest ODE function and plot in MATLAB

19 views (last 30 days)
What is the coolest ODE function and plot you came along?

Answers (1)

Davide Masiello
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
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.

Sign in to comment.

Tags

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!