Bubble models¶

A bubble model in jbubble has four parts:

1. An equation of motion (EoM): the ordinary differential equation for the radius \(R(t)\).
2. A gas model: the pressure inside the bubble.
3. A shell model: surface tension and the stresses of a coating.
4. A medium model: the viscous and elastic stresses of the surrounding liquid or tissue.

The equation of motion asks the other three for the liquid pressure at the bubble wall,

and relates it to the wall acceleration \(\ddot R\). Where an equation needs a time derivative, such as \(\mathrm{d}p_L/\mathrm{d}t\) in Keller-Miksis, JAX computes it by automatic differentiation. That's why any gas works with any shell, medium, and equation of motion, including models that you write yourself.

Every code block on this page runs as written, in order. The first block defines the parts that the later ones reuse:

import jax
import jax.numpy as jnp

from jbubble import run_simulation
from jbubble.bubble.eom import KellerMiksis, ModifiedRayleighPlesset, RayleighPlesset
from jbubble.bubble.gas import PolytropicGas
from jbubble.bubble.medium import NewtonianMedium
from jbubble.bubble.shell import NoShell
from jbubble.pulse import ToneBurst
from jbubble.pulse.shapes import Sine

water = NewtonianMedium(mu=1e-3, rho_L=998.0, c_L=1500.0)  # [Pa s], [kg/m³], [m/s]
parts = {
    "gas": PolytropicGas(gamma=1.4),
    "shell": NoShell(sigma=0.072),  # [N/m]
    "medium": water,
    "R0": 2e-6,  # [m]
    "P_amb": 101325.0,  # [Pa]
}

Equations of motion¶

Every equation of motion takes the same parameters: gas, shell, medium, R0, and P_amb. The medium holds the liquid's density rho_L and speed of sound c_L, which default to water at 20 °C; see Medium models.

Rayleigh-Plesset¶

The classic equation for a spherical bubble in an incompressible liquid, where \(p_\text{ac}(t)\) is the driving pressure:

Modified Rayleigh-Plesset¶

Adds the radiation damping of the gas pressure only, a first-order compressibility correction:

Keller-Miksis¶

Accounts for liquid compressibility to first order in the wall Mach number \(M = \dot R / c_L\):

Compare the equations¶

The following code drives the same free bubble at 400 kPa, hard enough for an inertial collapse, with three equations of motion:

pulse = ToneBurst(freq=1e6, pressure=400e3, shape=Sine(), cycle_num=3)
eoms = {
    "Rayleigh-Plesset": RayleighPlesset(**parts),
    "modified Rayleigh-Plesset": ModifiedRayleighPlesset(**parts),
    "Keller-Miksis": KellerMiksis(**parts),
}
for name, eom in eoms.items():
    result = run_simulation(eom, pulse)
    R = result.radius / eom.R0
    print(f"{name:25s} R_max/R0 = {R.max():.2f}, R_min/R0 = {R.min():.3f}")

Rayleigh-Plesset has no radiation damping, so it overestimates the growth and the collapse. For the full comparison, with plots, see the example Equations of motion.

Gas models¶

PolytropicGas is the standard model:

The polytropic exponent \(\gamma\) is 1 for an isothermal process, where heat flows faster than the bubble oscillates, and the adiabatic index for an adiabatic one: 1.4 for air and about 1.095 for the SF6 of lipid-coated contrast agents.

VanDerWaalsGas adds a hard core of radius \(h = h_\text{frac} R_0\), the radius of the gas compressed to its excluded volume, which stops a violent collapse from reaching \(R = 0\):

from jbubble.bubble.gas import VanDerWaalsGas

air = PolytropicGas(gamma=1.4)
air_with_core = VanDerWaalsGas(gamma=1.4, h_frac=0.1)  # h = 0.1 R0

The hard-core fraction h_frac depends on the gas species and its density at equilibrium.

P_gas0 isn't a parameter: the equation of motion computes it from the equilibrium at \(R_0\), where the gas pressure balances the ambient pressure and the Laplace pressure of the shell.

Shell models¶

Uncoated bubbles¶

NoShell contributes only the Laplace pressure \(2\sigma/R\) of the liquid's surface tension, 0.072 N/m for water.

Lipid shells¶

LipidShell models a lipid monolayer with a radius-dependent surface tension \(\sigma(R)\) and a surface-dilatational viscosity \(\kappa_s\):

The surface tension law carries the shell's elasticity. Three laws share the parameters R_buckle_ratio (the buckling radius as a fraction of \(R_0\)), chi (the elasticity \(\chi\)), and sigma_rupture:

In the elastic regime, all three follow \(\sigma = \chi\left[(R/R_b)^2 - 1\right]\) with \(R_b\) the buckling radius. The smoothed law stays within smoothing * ln(2) * sigma_rupture, about 0.5 mN/m at the defaults, of the piecewise law, and converges to it as smoothing goes to zero. The Gompertz law doesn't converge to the Marmottant law, so a \(\chi\) that you fit with it differs from a Marmottant \(\chi\); its Notes give the size of the difference.

A surface tension law is a Property: a function of the bubble state. The following code evaluates the three laws below, at, and above the buckling radius:

from jbubble.bubble.shell import (
    GompertzSurfaceTension,
    LipidShell,
    MarmottantSurfaceTension,
    SmoothMarmottantSurfaceTension,
)
from jbubble.bubble.state import BubbleState

shell_params = {"R_buckle_ratio": 0.98058, "chi": 0.5, "sigma_rupture": 0.072}
laws = {
    "Marmottant": MarmottantSurfaceTension(**shell_params),
    "smoothed": SmoothMarmottantSurfaceTension(**shell_params),
    "Gompertz": GompertzSurfaceTension(**shell_params),
}
for name, law in laws.items():
    values = [law(BubbleState(R=x * 2e-6, R0=2e-6)) for x in (0.97, 1.0, 1.05)]
    print(f"{name:10s}", " ".join(f"{1e3 * v:6.2f}" for v in values), "mN/m")

lipid = LipidShell(sigma=laws["smoothed"], kappa_s=7.5e-9)  # kappa_s [N s/m]

Thick shells¶

ThickShell models a polymer or protein shell of thickness \(d_s\), shear modulus \(G_s\), and shear viscosity \(\mu_s\) as an incompressible viscoelastic layer (Church 1995):

where \(V_s = R_0^3 - (R_0 - d_s)^3\). For a thin shell, these reduce to the model of Hoff et al. (2000), which adds a stiffness of \(12 G_s d_s / R_0\) per unit radial strain.

from jbubble.bubble.shell import ThickShell

polymer = ThickShell(sigma=0.04, d_s=20e-9, G_s=11.7e6, mu_s=0.45)

Medium models¶

The two elastic laws agree for small strains, \(4G(R - R_0)/R_0\), but only the neo-Hookean law holds at finite strain. For a power-law liquid, mu is the consistency index \(K\) in Pa sⁿ, and eps regularises the shear rate where \(\dot R\) changes sign.

Every medium also holds the liquid's density rho_L [kg/m³] and speed of sound c_L [m/s], which the equations of motion and the emission models read. They're keyword-only and default to water at 20 °C: rho_L=998.0 and c_L=1500.0. They're plain constants, not Property fields, because the equations of motion assume a liquid of constant density and sound speed.

from jbubble.bubble.medium import NeoHookeanMedium, PowerLawMedium

tissue = NeoHookeanMedium(mu=5e-3, G=10e3, rho_L=1060.0, c_L=1540.0)  # [Pa s], [Pa], [kg/m³], [m/s]
shear_thinning = PowerLawMedium(mu=1.5e-2, n_exp=0.7)  # K [Pa s^n]

Properties: constant, state-dependent, or learned¶

Every physical coefficient of a gas, shell, or medium model, such as sigma, kappa_s, G, or gamma, except the liquid's rho_L and c_L, is a Property: a function from the bubble state to a scalar. A plain number becomes a ConstantProperty. Pass any other Property to make the coefficient depend on the state, for example a shear modulus that stiffens with strain:

from jbubble.bubble.medium import KelvinVoigtMedium
from jbubble.bubble.property import Property


class StrainStiffening(Property):
    """Shear modulus G0 (1 + alpha ε²) at the radial strain ε = R/R0 - 1."""

    G0: float  # [Pa]
    alpha: float

    def __call__(self, state):
        strain = state.R / state.R0 - 1.0
        return self.G0 * (1.0 + self.alpha * strain**2)


gel = KelvinVoigtMedium(mu=1e-3, G=StrainStiffening(G0=5e3, alpha=20.0))

NeuralProperty wraps a neural network that maps \(R/R_0\) to the value. Bound its output with the network's final activation, for example to keep a surface tension between 0 and the rupture value:

import equinox as eqx

from jbubble.bubble.property import NeuralProperty

net = eqx.nn.MLP(
    in_size=1,
    out_size=1,
    width_size=16,
    depth=2,
    final_activation=lambda x: 0.072 * jax.nn.sigmoid(x),  # [N/m]
    key=jax.random.key(0),
)
learned_shell = LipidShell(sigma=NeuralProperty(net=net), kappa_s=7.5e-9)

To fit such a network to data, see Learn a constitutive law with a neural network. To write a new medium or shell model rather than a coefficient, subclass MediumModel or ShellModel; the example Custom physics shows how.

Set the initial state¶

By default, a simulation starts at rest at the equilibrium radius, eom.initial_state(). To start elsewhere, pass state0 to run_simulation. The following code releases a bubble at rest from 1.5 times its equilibrium radius, with no drive:

eom = KellerMiksis(**parts)
silence = ToneBurst(freq=1e6, pressure=0.0, shape=Sine(), cycle_num=5)
state0 = eom.initial_state(R=1.5 * eom.R0)  # R_dot defaults to 0
result = run_simulation(eom, silence, state0=state0)
print(f"R_min/R0 after release: {result.radius.min() / eom.R0:.3f}")

state0 holds the radius R, the wall velocity R_dot, and the equilibrium values R0 and P_gas0, which stay constant during the solve. BubbleState(R=1.5 * R0) gives the same start: a zero R0 or P_gas0 means "unset", and the solver fills it from the equation of motion. For an empty cavity, set P_gas0 to a tiny positive value, such as 1e-12.

Choose models for your application¶

The presets free_bubble, lipid_bubble, and thick_shell_bubble assemble the first, third, and fourth rows with cited parameters. Small or viscous bubbles can make the equations stiff; see Solvers and stiffness.