Unexplained gain in my code when solving coupled differential equations using solve_ivp() in python

Question

I am trying to solve two coupled differential equations in python using the solve:ivp() function from scipy.integrate but I get an unexplained gain on one of my terms that breaks the equation. It is for the movement of bubbles in an acoustic field and i am trying to solve for the radius, rate of change of radius, position (z) and velocity (u_b); however, I get that the u_b value increases exponentially to unrealistic values crashing the solver. The code that i am using is the following:

import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from scipy.integrate import OdeSolver
import matplotlib.pyplot as plt

def solver_ode(t, y):
    rho_l = 997.4  # Liquid density (kg/m^3)
    mu_l = 8.9e-4  # Dynamic viscosity (Pa·s)
    sigma = float(0.072)  # Surface tension (N/m)
    P_l0 = 1e5  # Far-field pressure (Pa)
    gamma = 1.4  # Polytropic index
    R0 = 0.00025  # Initial bubble radius (m)
    f = 600  # Acoustic frequency (Hz)
    omega = 2 * np.pi * f  # Angular frequency (rad/s)
    c = 1500  # Speed of sound in medium (m/s)
    lamda = c / f  # Acoustic wavelength (m)
    k = 2 * np.pi / lamda  # Wavenumber (1/m)
    g = 9.81  # Gravity (m/s^2)
    dPa = 60000  # Acoustic amplitude (Pa)
    rho_b = 1.2  # Bubble density (kg/m^3)

    R, R_dot, z, u_b = y  # Unpacking condition


    # Gas pressure

    P_v = 0
    Pg0 = ((2 * sigma) / R0) + P_l0 - P_v
    Pg = Pg0 * (R0 / R) ** (3 * gamma)

    # Acoustic pressure
    P_ac = dPa * np.cos(omega * t) * np.sin(k*z)
    P_inf = P_ac + P_l0
    P_b = Pg + P_v

    # Radial acceleration
    R_ddot = (1 / (rho_l * R)) * (P_b - P_inf) - \
             (3 / (2 * R)) * (R_dot ** 2) - \
             ((4 * mu_l) / (rho_l * (R ** 2))) * R_dot - \
             (2 * sigma / (rho_l * (R ** 2)))
    
    V = (4 / 3) * np.pi * R ** 3
    V_dot = 4 * np.pi * R**2 * R_dot
    m_b = rho_b * ((4 / 3) * np.pi * R0 ** 3)    
    A = 4 * np.pi * R ** 2
    

    # Forces
    u_l = ((-k * dPa) / (omega * rho_l)) * np.sin(omega * t) * np.cos(2 * np.pi * k*z)
    u_l_dot = ((-k * dPa) / (rho_l)) * np.cos(omega * t) * np.cos(2 * np.pi * k*z)

    F_bj = -V * (dPa * k * np.cos(omega * t) * np.cos(2 * np.pi * k*z))
    V_dot = 0

    u_b_dot = (1/(m_b + 0.5*rho_l*V))*(F_bj - 0.5*rho_l*V_dot*(u_b - u_l) //
                + 0.5*rho_l*V*(-k*(dPa/(omega*rho_l)*(omega*np.cos(omega*t)*np.cos(k*z) -                k*np.sin(omega*t)*np.sin(k*z)*u_b))) - 0.5*rho_l*(u_b - u_l)**2*A*(27*(abs(mu_l/(2*rho_l*(u_b - u_l)*R0)))**0.78) + V*(rho_l - rho_b)*g)


    return [R_dot, R_ddot, u_b, u_b_dot]


# Initial conditions: [R, R_dot, pos, u_b]
y0 = [0.00025, 0, 0.001, 0.00001]  # Start at a quarter wavelength

t_span = (0, 1)  # Time span (s)
t_eval_1 = np.linspace(0, 1, 1000)  # Time steps

# Solve the system

solution = solve_ivp(
    solver_ode, t_span, y0,  method='RK45', max_step=0.0001,
    rtol=1e-6, atol=1e-9
)

print(solution)

Thank you in advance for any help.

I have tried altering the parameters and simplying the eqautions to just solve for the postion and velocity of the equation by keeping the radius the same. This has not work and the value of the velocity u_b still grows exponentially eventually crashing the function. I have tried the other solvers available in solve_ivp() but the same happens.

Answer 1

OK, I'm guessing that the source of the equations is https://doi.org/10.1016/j.ijmultiphaseflow.2021.103733 , so that's what I'm coding to.

Your expression for the bubble acceleration, u_b_dot, is far too complicated and over-bracketed; it is impossible to decipher. You need to break it down into the original forces. (These are given in the reference above, though in that reference: the added-mass force has the wrong sign; the Reynolds number is positive, so needs an abs(); the drag force, whilst proportional in magnitude to the square of the relative velocity, must always have the opposite sign). Area A is used in the drag force, so should be the projected area normal to the flow, not the surface area of a sphere as you have coded.

The following version of your code aims for one of the bubble diameters in Figure 2a of the reference above. It is VERY sensitive to bubble diameter and acoustic frequency. The maximum timestep must be set small enough to resolve the high acoustic frequency.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

rho_l = 997.4                                                   # Liquid density (kg/m^3)
mu_l = 8.9e-4                                                   # Dynamic viscosity (Pa s)
sigma = 0.072                                                   # Surface tension (N/m)
P_l0 = 1e5                                                      # Far-field pressure (Pa)
gamma = 1.4                                                     # Polytropic index
R0 = 0.000565                                                   # Initial bubble radius (m)
f = 14580                                                       # Acoustic frequency (Hz)
omega = 2 * np.pi * f                                           # Angular frequency (rad/s)
c = 1500                                                        # Speed of sound in medium (m/s)
lamda = c / f                                                   # Acoustic wavelength (m)
k = 2 * np.pi / lamda                                           # Wavenumber (1/m)
g = 9.81                                                        # Gravity (m/s^2)
dPa = 20000                                                     # Acoustic amplitude (Pa)
rho_b = 1.2                                                     # Bubble density (kg/m^3)

m_b = rho_b * (4 / 3) * np.pi * R0 ** 3                         # Bubble mass (kg)
P_v = 0                                                         # Vapour pressure (Pa)
Pg0 = 2 * sigma / R0 + P_l0 - P_v                               # Initial gas pressure (Pa)


def solver_ode(t, y):
    R, R_dot, z, u_b = y  # Unpacking condition

    # Gas pressure
    Pg = Pg0 * (R0 / R) ** (3 * gamma)
    P_b = Pg + P_v                                              # Bubble pressure

    # Acoustic pressure
    P_ac = dPa * np.cos(omega * t) * np.sin(k * z)
    P_inf = P_ac + P_l0

    # Radial acceleration                                       # Rayleigh-Plesset equation
    R_ddot = (  (P_b-P_inf)/rho_l - 1.5*R_dot**2 - 4*mu_l*R_dot/rho_l/R - 2*sigma/rho_l/R   ) / R
    
    # Geometric quantities
    V = (4 / 3) * np.pi * R ** 3
    V_dot = 4 * np.pi * R**2 * R_dot
    A = np.pi * R ** 2                                          # Projected area (used in drag definition)
    
    # Liquid motion
    u_l     = -(k / omega) * ( dPa / rho_l ) * np.sin(omega * t) * np.cos(k * z)
    u_l_dot = - k          * ( dPa / rho_l ) * np.cos(omega * t) * np.cos(k * z)
    u_r = u_b - u_l                                             # Relative velocity

    # Forces
    # Bjerknes force
    dPldz = k * dPa * np.cos(omega * t) * np.cos(k * z)         # Pressure gradient
    F_bj = -V * dPldz

    # Drag force
    Re = 2 * rho_l * abs(u_r) * R / mu_l                        # Reynolds number
    cD = 27 / Re ** 0.78                                        # Drag coefficient
    F_D = -0.5 * rho_l * u_r * abs( u_r ) * A * cD

    # Buoyancy force
    F_buoy = ( rho_l - rho_b ) * V * g

    # Added-mass force (part treated implicitly via m_a)
    m_a = 0.5 * rho_l * V
    F_a = -0.5 * rho_l * u_r * V_dot + 0.5 * rho_l * V * u_l_dot

    u_b_dot = ( F_bj + F_D + F_buoy + F_a ) / ( m_b + m_a )     # Acceleration = force / (effective mass)

    return [R_dot, R_ddot, u_b, u_b_dot]



y0 = [R0, 0, 0.3 * lamda, 0.00001]                              # Initial conditions [R, R_dot, pos, u_b]
t_span = (0, 0.6)                                               # time span (s)
solution = solve_ivp( solver_ode, t_span, y0, max_step=0.00001 )    # solve
print(solution)

plt.plot( solution.t, solution.y[2] / lamda, label=f'R0={1000*R0} mm' )
plt.xlim( 0, 0.6 )
plt.ylim( 0.3, 0.8 )
plt.xlabel( 't (s)' )
plt.ylabel( r'z/$\lambda$' )
plt.legend()
plt.show()