I have been trying to numerically solve the Rayleigh-Plesset equation for a sonoluminescence bubble in Python. You can read about this phenomenon here: https://iopscience.iop.org/article/10.1088/0143-0807/34/3/679.
The Rayleigh-Plesset equation is a non-linear ODE, which can be solved to find the Radius of a bubble subject to non-linear oscillations due to an external driving sound wave (Sonoluminescence). Here is the form of the equation I used: \begin{equation}\ddot{R} = \frac{1}{R\rho}\Big(P_g-P_0-P(t)-4\mu\frac{\dot{R}}{R}-\frac{2\sigma}{R}+\big(\frac{2\sigma}{R_0}+P_0\big)\big(\frac{R_0}{R}\big)^{3\kappa}\Big)-\frac{3\dot{R}^2}{2R}\end{equation}
I rewrote this as a system of differential equations (so that ODEINT would process it): \begin{cases} \dot{R}=u\\ \dot{u} = \frac{1}{R\rho}\Big(P_g-P_0-P(t)-4\mu\frac{u}{R}-\frac{2\sigma}{R}+\big(\frac{2\sigma}{R_0}+P_0\big)\big(\frac{u}{R}\big)^{3\kappa}\Big)-\frac{3u^2}{2R} \end{cases}
I used the following parameters and initial conditions: \begin{cases}R(0) = 2.0\times10^{-6} \text{ m}\\u(0) = 0 \text{ ms}^{-1}\\ \rho = 10^3 \text{ kg m}^3\\ \sigma = 7.25\times10^2 \text{ Nm}^{-1}\\ \mu = 8.9\times10^{-4}\text{ Pa s}\\ P_g = 2330 \text{ Pa}\\ P_0 = 10^5 \text{ Pa}\\ \kappa = 1.33 \end{cases}
The driving pressure from the sound waves is a sine function: $P(t) = 1367900\sin{(2\pi(25600t))}$Here is the python code I have written. It plots both the radius of the bubble and its radial velocity.
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint
# define equations
def equation(y0, t):
R, u = y0
return u, (P_g-P_0-70000*np.sin(2*np.pi*31700*t)-2*sigma/R-
4*miu*u/R+(2*sigma/R_0+P_0-P_g)*(R_0/R)**(3*k))/(R*rho)-3*u**2/(2*R)
# initial conditions
R_0 = 0.000005
u_0 = 0
# parameters
rho = 1000
sigma = 0.0728
miu = 8.9*10**(-4)
P_g = 2330
P_0 = 10000
k = 1.33
time = np.arange(0, 0.0003, 0.00000000025)
R_1 = odeint(equation, [R_0, u_0], time)
V = R_1[:,1]
R = R_1[:,0]*10**6
mtimes = time*10**6
#plot results
fig, ax1 = plt.subplots()
ax1.set_xlabel("T/$\mu$s")
ax1.set_ylabel("R/$\mu$m", color = "red")
ax1.plot(mtimes, R, linewidth = 0.7, label = "Bubble Radius", color =
"red")
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-
axis
ax2.set_ylabel("$\dot{R}$/$ms^{-1}$") # we already handled the x-label
with ax1
ax2.plot(mtimes, V, linestyle = "dashed", color = "black", linewidth =
0.7, label = "Radial Velocity Bubble")
ax1.legend()
ax2.legend(loc = "lower right")
fig.tight_layout() # otherwise the right y-label is slightly clipped
plt.show()
Here is the output I got:
This is clearly wrong. Firstly, the maximum bubble radius is around 600$\mu$m, which is way too high. Moreover, if I were to extend the time range, the periodic behaviour would not be exhibited, something that is present in all Sonoluminescence papers I have read so far. Not only, increasing the time range (say $0<t<0.003$) gives a weird output:
As you can see the radius bubble graph remains the same, but the radial velocity changes completely. All information from the previous graph is lost, and all I see is this spike.
I think this problem may have to do with the equation in itself, hence why I am posting this on Computational Sciences. I think the driving pressure should be a cosine function, and not a sine one (although I don't really see what difference it would make). Maybe my pressure parameters are wrong?
EDIT: I have updated the code with the correct ODE, and here is the output:
. This coincides with literature, for example:
So i think I've got the first oscillation down, all I need to do now is add all the other oscillations that follow. How could the other oscillations be added?
