The coolest ODE function and plot in MATLAB

What is the coolest ODE function and plot you came along?

Answer 1

[Davide Masiello]

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

[Shishir Raut]

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.