# 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: $$$$\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}$$$$

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?

• The reason you aren't seeing more oscillations is because you are missing the radiation term (energy loss to acoustics). Your "*(R_0/R)**(3*k)" term just needs to be modified to "*(R_0/R)**(3*k)*(1-u/c)" where c is the speed of sound in your liquid. Apr 1, 2021 at 15:29

I see at least one important problem. On the right hand side you have a term that looks like

$$P_o \left( \frac{\dot{R}}{R} \right)^{3 \kappa}$$

This term is dimensionally inconsistent with the other terms in the brackets, which have dimensions of a pressure. This term should actually be

$$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.

• Thank you. I have updated the code, and it (finally) gives a "normal" output. Unfortunately however, the code only gives the first oscillation, and all the subsequent ones are missing... Maybe I should try with MATLAB (ODE45 package might work maybe)? Dec 22, 2019 at 9:21
• I was unable to reproduce your result using your code. However, you'll want to use adaptive time stepping near the collapse times. Surface tension and viscosity are not required to observe this behavior (aside from the decaying amplitude of the growth phase).
– user20857
Dec 22, 2019 at 11:49