# Numerical Solution to Rayleigh Plesset Equation in Python

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) 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?

$$P_o \left( \frac{\dot{R}}{R} \right)^{3 \kappa}$$
$$P_o \left( \frac{R_o}{R} \right)^{3 \kappa}$$
The paper you link at the beginning has this correct in equation (6), assuming $$h$$ is zero. I also that advise you mess around with the RPE with $$P(t) = const.$$, $$\mu = \sigma = 0$$ for a bit at first.