Hello dear all,
Continuing my work on the Rayleigh-Plesset equation, I wrote this code:
using OrdinaryDiffEq, Plots
#SI units
#external parameters
Pₕ = 1.01325e5 #hydrostatic pressure (Pa)
#Parameters (liquid = water)
Pᵥ = 2.5e3 #vapour pressure (Pa)
σ = 73.0e-3 #surface tension (N/m)
ρ = 1000.0 #volumic mass (Kg.m-3)
η = 1.0e-3 #dynamic viscosity (Pa.s)
#Parameter (gas = Ar)
γ = 1.666 #polytropic ratio (-)
#Acoustic Parameters
Pₐ = 1.3e5 #acoustic pressure (Pa)
fₐ = 2.0e4 #acoustic frequency (s-1)
#Initial conditions
R₀ = [4.0e-6] #initial radius (m)
dR₀ = [1.0e-2] #initial radial velocity (ms-1)
tspan = (0.0,150.0e-6) #(s) -- 1 acoustic cycle = 50 μs
#Num
function rpnnp(ddR, dR, R, p, t)
@. ddR = (1/ρ*R)*((Pₕ-Pᵥ-(2*σ/R₀))*(R₀/R)^(3*γ)-(4*η*dR/R)-(2*σ/R)-Pₕ+Pᵥ-p(t)) - (1.5*dR*dR)/R
end
P = t->Pₐ*sin(2*pi*fₐ*t + pi/2) #external acoustic wave
probleme = SecondOrderODEProblem(rpnnp, dR₀, R₀, tspan, P)
solution = solve(probleme, DPRKN6())
@show solution.t
@show solution.u
#graph
plot(solution, vars=(2), linewidth=2,
xaxis = "t ", yaxis = "R",
framestyle = :box)
I am expecting a result like in the 4th graph concerning the “radius versus the time” curve:
However with the code as above the bubble radius evolves hardly.
Do you have a suggestion since I do not see where my error is.
Thank you in advance,
Thierry