Numerically Solving the Rayleigh Plesset Equation

  • Python
  • Thread starter KDPhysics
  • Start date
  • Tags
    Rayleigh
In summary: 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-P_g\big)\big(\frac{u}{R}\big)^{3\kappa}\Big)-\frac{3u^2}{2R}##
  • #1
KDPhysics
74
23
I have been trying to numerically solve the Rayleigh-Plesset equation for a sonoluminescing 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-P_g\big)\big(\frac{\dot{R}}{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)=70000\sin(2\pi(31700t))$$​

Here is the python code I have written. It plots both the radius of the bubble and its radial velocity.
Python:
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.00001
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.00008, 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()

The results from this code however doesn't coincide with literature. Indeed, my code only gives the first collapse, whereas it should show various non-linear oscillations. Not only, the velocity plot shows that the bubble expands at supersonic speeds (which also doesn't coincide with literature given the present parameters).

Here is what the output should look like:
ejp456493f1_online.jpg


Here is what the code gives:
Schermata 2019-12-21 alle 14.02.19.png


Please CompSci geniuses help me ...

EDIT: here is the stack exchange thread I created regarding the same question. Hopefully, someone will answer on one of the two forums.
 
Last edited:
Technology news on Phys.org
  • #2
Check your work again, you have lost a ## P_g ## between equation (1) and its restatement as a first order system, although you seem to have found it again in the code. However you have an R_0 where it looks like you should have a u.

Once you have got this right I expect the plot should take similar values to the original plot up to just before what appears to be a singular point at t ≈ 17.5μs, but I don't think you should be surprised if it doesn't get any further. You probably need to transform your system into something better behaved. You might also have more success with scipy.integrate.solve_ivp which is replacing the deprecated scipy.integrate.odeint.
 
  • #3
Equation 1 and your python code don't agree.
Equation 1
\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-P_g\big)\big(\frac{\dot{R}}{R}\big)^{3\kappa}\Big)-\frac{3\dot{R}^2}{2R}\end{equation}

Your python code with difference noted in bold red:
(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)

It seems to me that your python code should have ##\frac u R##, not ##\frac {R_0} R##.

Also, but not relevant to the problem you're seeing, There's a difference between Equation 1 (above) and what you wrote following "I rewrote this as a system of differential equations ..."

##
\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}##
Equation 1 has ##(\frac{2\sigma}{R_0} - P_0 - P_g)##, but the equation just above is missing the ##-P_g## term.
 
  • #4
o_O
Turns out the ODE i wrote in the code was the correct one, and that the one i wrote in latex is incorrect...
I might try out scipy.integrate.solve_ivp
 

Related to Numerically Solving the Rayleigh Plesset Equation

1. What is the Rayleigh Plesset Equation?

The Rayleigh Plesset equation is a mathematical model used to describe the behavior of a gas bubble in a liquid under the influence of an external pressure, such as sound waves or changes in depth. It takes into account factors such as surface tension, viscosity, and compressibility of the gas and liquid.

2. Why is it important to numerically solve the Rayleigh Plesset Equation?

Numerically solving the Rayleigh Plesset equation allows us to accurately predict the behavior of gas bubbles in various situations, such as in medical procedures or in oceanography. It also helps us understand the physical mechanisms behind the behavior of gas bubbles, which has implications in fields such as fluid dynamics and acoustics.

3. What methods are typically used to numerically solve the Rayleigh Plesset Equation?

The most commonly used methods for numerically solving the Rayleigh Plesset equation are finite difference methods, finite element methods, boundary element methods, and spectral methods. Each method has its own advantages and disadvantages, and the choice of method depends on the specific problem being solved.

4. What are some challenges in numerically solving the Rayleigh Plesset Equation?

One of the main challenges in numerically solving the Rayleigh Plesset equation is accurately modeling the complex dynamics of gas bubbles, such as the formation and collapse of bubbles, and the effects of non-linearities and instabilities. Another challenge is dealing with the large range of length and time scales involved in such problems.

5. What are some applications of numerically solving the Rayleigh Plesset Equation?

Numerically solving the Rayleigh Plesset equation has a wide range of applications in fields such as biomedical engineering, oceanography, and acoustics. It is used to study the behavior of gas bubbles in ultrasound imaging, drug delivery, and decompression sickness. It is also used to understand the behavior of gas bubbles in the ocean, such as in sonar applications and the formation of acoustic noise.

Similar threads

  • Other Physics Topics
Replies
1
Views
2K
Replies
1
Views
717
Replies
1
Views
754
  • Advanced Physics Homework Help
Replies
7
Views
1K
  • Special and General Relativity
Replies
11
Views
447
Replies
13
Views
1K
  • Advanced Physics Homework Help
Replies
11
Views
2K
Replies
2
Views
585
  • Programming and Computer Science
Replies
13
Views
3K
  • Special and General Relativity
Replies
15
Views
2K
Back
Top