The coolest ODE function and plot in MATLAB
9 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
What is the coolest ODE function and plot you came along?
0 Kommentare
Antworten (1)
Davide Masiello
am 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 Kommentar
Shishir Raut
am 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.
Siehe auch
Kategorien
Mehr zu General Applications finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
