New functionality: Power-law liquids

.github/workflows/test_run.yml

@@ -108,6 +108,8 @@ jobs:
       run: ${{github.workspace}}/examples/ultrasound/build/ultrasound_apecss -options ${{github.workspace}}/examples/ultrasound/zener.apecss -tend 5e-6 -freq 1e6 -amp 1e6
     - name: Run ultrasound (Oldroyd-B)
       run: ${{github.workspace}}/examples/ultrasound/build/ultrasound_apecss -options ${{github.workspace}}/examples/ultrasound/oldroydb.apecss -tend 3e-6 -freq 3e6 -amp 400e3
+    - name: Run ultrasound (power-law)
+      run: ${{github.workspace}}/examples/ultrasound/build/ultrasound_apecss -options ${{github.workspace}}/examples/ultrasound/powerlaw.apecss -freq 636619 -amp 25331 -tend 25e-6
 
   linux:
     runs-on: ubuntu-latest
@@ -206,3 +208,5 @@ jobs:
       run: ${{github.workspace}}/examples/ultrasound/build/ultrasound_apecss -options ${{github.workspace}}/examples/ultrasound/zener.apecss -tend 5e-6 -freq 1e6 -amp 1e6
     - name: Run ultrasound (Oldroyd-B)
       run: ${{github.workspace}}/examples/ultrasound/build/ultrasound_apecss -options ${{github.workspace}}/examples/ultrasound/oldroydb.apecss -tend 3e-6 -freq 3e6 -amp 400e3
+    - name: Run ultrasound (power-law)
+      run: ${{github.workspace}}/examples/ultrasound/build/ultrasound_apecss -options ${{github.workspace}}/examples/ultrasound/powerlaw.apecss -freq 636619 -amp 25331 -tend 25e-6

.vscode/settings.json

@@ -10,6 +10,8 @@
         "stdio.h": "c",
         "*.tcc": "c",
         "random": "c",
-        "ratio": "c"
-    }
+        "ratio": "c",
+        "apecss.h": "c"
+    },
+    "cmake.sourceDirectory": "${workspaceFolder}/lib"
 }

README.md

@@ -23,7 +23,7 @@ Key features of APECSS are:
 - Acoustic interactions of multiple bubbles under different assumptions (incompressible, quasi-acoustic).
 - Interbubble interactions through their acoustic emissions
 - Prediction of the formation and attenuation of shock fronts emitted by the bubble.
-- Viscoelastic media (Kelvin-Voigt, Zener, Oldroyd-B).
+- Non-Newtonian (power-law) and viscoelastic (Kelvin-Voigt, Zener, Oldroyd-B) media.
 - Lipid monolayer coating of the bubble as used for ultrasound contrast agents.
 - APECSS has, aside from the C standard library, no external dependencies.
 
@@ -91,7 +91,7 @@ If you use APECSS for your scientific work, please consider citing the [paper](h
 
     F. Denner and S. Schenke, APECSS: A software library for cavitation bubble dynamics and acoustic emissions. Journal of Open Source Software 8 (2023), 5435. https://doi.org/10.21105/joss.05435 
 
-and, in the interest of reproducibility, the version of APECSS you've used for your work, e.g.
+and, in the interest of reproducibility, the version of APECSS you've used for your work, for instance
 
     F. Denner, S. Schenke and P. Coulombel, APECSS (v1.7), Zenodo (2024). https://doi.org/10.5281/zenodo.13850831
 

documentation/chapters/bubble.tex

@@ -145,20 +145,21 @@ \section{The gas}
 \section{The liquid}
 \label{sec:liquid}
 
-In the same way that every bubble contains a gas, in APECSS every bubble is surrounded by a liquid, which requires to select an appropriate equation of state and fluid type, as well as define meaningful properties.
+In the same way that every bubble contains a gas, in APECSS every bubble is surrounded by a liquid, which requires to select an appropriate equation of state and liquid type, as well as define meaningful properties.
 
 \vspace{0.8em}
 
 \noindent
-\begin{tabular}{p{0.1\textwidth} p{0.35\textwidth} p{0.5\textwidth}}
+\begin{tabular}{p{0.08\textwidth} p{0.385\textwidth} p{0.485\textwidth}}
     \textbf{Section} &\textbf{Command} & \textbf{Description} 
 \vspace{1mm} \\ \hline
 {\tt LIQUID} & {\tt EoS Tait} & The Tait EoS is applied to the liquid. Only relevant for the Gilmore model and acoustic emissions based on the Kirkwood-Bethe hypothesis.\\ 
 & {\tt EoS NASG} & The Noble-Abel-stiffened-gas EoS is applied to the liquid. Only relevant for the Gilmore model and acoustic emissions based on the Kirkwood-Bethe hypothesis.\\ 
-& {\tt LiquidType Newtonian} & Newtonian fluid. This is the default.\\ 
+& {\tt LiquidType Newtonian} & Newtonian liquid. This is the default.\\ 
 & {\tt LiquidType KelvinVoigt} & Kelvin-Voigt solid.\\ 
 & {\tt LiquidType Zener} & Zener solid.\\ 
-& {\tt LiquidType OldroydB} & Oldroyd-B (or upper-convected Maxwell) fluid.\\ 
+& {\tt LiquidType OldroydB} & Oldroyd-B (or upper-convected Maxwell) liquid.\\
+& {\tt LiquidType PowerLaw} & Power-law liquid.\\ 
 & {\tt PolytropicExponent <float>} & Polytropic exponent $\Gamma_{\ell}$.\\
 & {\tt ReferencePressure <float>} & Reference pressure $p_{\ell,\text{ref}}$.\\
 & {\tt ReferenceDensity <float>} & Reference density $\rho_{\ell,\text{ref}}$.\\
@@ -169,6 +170,8 @@ \section{The liquid}
 & {\tt PolymerViscosity <float>} & Polymer viscosity $\eta_\ell$ associated with viscoelasticity.\\
 & {\tt ShearModulus <float>} & Shear modulus $G_\ell$ associated with viscoelasticity.\\
 & {\tt RelaxationTime <float>} & Relaxation time $\lambda_\ell$ associated with viscoelasticity.\\
+& {\tt PowerLawExponent <float>} & Exponent $n$ of a power-law liquid.\\
+& {\tt PowerLawConsistencyCoeff <float>} & Consistency coefficient $k$ of a power-law liquid.\\
  \hline
 \end{tabular} \vspace{1em}
 
@@ -178,10 +181,21 @@ \section{The liquid}
 \end{equation}
 where $p_\text{G}$ is the gas pressure, see Section \ref{sec:gas}, $\sigma$ is the surface tension coefficient of the interface, see Section \ref{sec:interface}, and $\mu_\ell$ is the liquid (Newtonian) viscosity. The derivative of Eq.~\eqref{eq:pL} follows as
 \begin{equation}
-  \dot{p}_\mathrm{L} = \dot{p}_\mathrm{G} + \frac{2 \, \sigma \, \dot{R}}{R^2} + 4 \, \mu_\ell \left(\frac{\dot{R}^2}{R^2} - \frac{\ddot{R}}{R}\right).
+  \dot{p}_\mathrm{L} = \dot{p}_\mathrm{G} + \frac{2 \sigma \dot{R}}{R^2} + 4 \, \mu_\ell \left(\frac{\dot{R}^2}{R^2} - \frac{\ddot{R}}{R}\right).
   \label{eq:dotpL}
 \end{equation}
 
+The pressure of a power-law liquid at the wall of a spherical bubble is given as \citep{Kaykanat2024}
+\begin{equation}
+  p_\text{L} = p_\text{G} - \frac{2 \sigma}{R} - \frac{4 \, k^\star}{n} \, |\dot{R}|^{n-1} \frac{\dot{R}}{R^n}, \label{eq:pL_powerlaw}
+\end{equation}
+where $p_\text{G}$ is the gas pressure, see Section \ref{sec:gas}, $\sigma$ is the surface tension coefficient of the interface, see Section \ref{sec:interface}, $k^\star = k(2\sqrt{3})^{n-1}$ is an extended consistency coefficient \citep{Kaykanat2024}, and $n$ is the power-law exponent. The derivative of Eq.~\eqref{eq:pL_powerlaw} follows as
+\begin{equation}
+  \dot{p}_\mathrm{L} = \dot{p}_\mathrm{G} + \frac{2 \sigma \dot{R}}{R^2} + 4 \, k^\star \, |\dot{R}|^{n-3} \dot{R}^2 \, \left(\frac{\dot{R}^2}{R^2} - \frac{\ddot{R}}{R}\right).
+  \label{eq:dotpL_powerlaw}
+\end{equation}
+The power-law liquid, described by Eqs.~\eqref{eq:pL_powerlaw} and \eqref{eq:dotpL_powerlaw}, reduces to a Newtonian liquid, described by Eqs.~\eqref{eq:pL} and \eqref{eq:dotpL}, for $n=1$ and $k = \mu_\ell$.
+
 \subsection{Equation of state}
 
 For the Gilmore model \eqref{eq:gilmore} and the acoustic emissions based on the Kirkwood-Bethe hypothesis (see Section \ref{sec:emissionskb}), an equation of state (EoS) for the liquid has to be defined. Two EoS are currently available in APECSS: the Tait EoS and the NASG EoS. 
@@ -285,7 +299,7 @@ \subsubsection{Zener model}
 
 \subsubsection{Oldroyd-B model}
 
-The Oldroyd-B model is a widely used constitutive model for viscoelastic fluids. 
+The Oldroyd-B model is a widely used constitutive model for viscoelastic liquids. 
 Following the work of \citet{Jimenez-Fernandez2005}, the liquid pressure at the bubble wall including the Oldroyd-B model is given as
 \begin{equation}
      p_\text{L} = p_\text{G} - \frac{2 \sigma}{R} - 4 \mu_\ell \frac{\dot{R}}{R} + \mathcal{S} \label{eq:pL_OldroydB}.
@@ -301,7 +315,7 @@ \subsubsection{Oldroyd-B model}
     \mathcal{S}_{1,n+1} &= \mathcal{S}_{1,n} + \Updelta t \frac{-\left(4 \lambda_\ell \dfrac{\dot{R}}{R}+1\right) \mathcal{S}_{1,n}- 2 \eta_\ell \dfrac{\dot{R}}{R}}{\lambda_\ell+\Updelta t} \label{eq:ode_oldroydB1disc} \\
     \mathcal{S}_{2,n+1} &= \mathcal{S}_{2,n} + \Updelta t \frac{-\left(\lambda_\ell \dfrac{\dot{R}}{R}+1\right) \mathcal{S}_{2,n}- 2 \eta_\ell \dfrac{\dot{R}}{R}}{\lambda_\ell+\Updelta t}.\label{eq:ode_oldroydB2disc}
 \end{align}
-For $\lambda_\ell = 0$ Eqs.~\eqref{eq:ode_oldroydB1disc} and \eqref{eq:ode_oldroydB2disc} still give a meaningful result and reduce to a Newtonian fluid with $\mathcal{S} = - 4 \eta_\ell \dot{R}/R$.
+For $\lambda_\ell = 0$ Eqs.~\eqref{eq:ode_oldroydB1disc} and \eqref{eq:ode_oldroydB2disc} still give a meaningful result and reduce to a Newtonian liquid with $\mathcal{S} = - 4 \eta_\ell \dot{R}/R$.
 
 \section{The interface}
 \label{sec:interface}

documentation/references.bib

@@ -161,6 +161,20 @@ @article{Jimenez-Fernandez2005
   doi     = {10.1016/j.ultras.2005.03.010}
 }
 
+@article{Kaykanat2024,
+  title   = {Shape Stability of a Microbubble in a Power--Law Liquid},
+  author  = {Kaykanat, S. Ilke and Uguz, Kerem},
+  year    = {2024},
+  month   = sep,
+  journal = {The European Physical Journal Special Topics},
+  volume  = {233},
+  number  = {8-9},
+  pages   = {1625--1635},
+  issn    = {1951-6355, 1951-6401},
+  doi     = {10.1140/epjs/s11734-024-01174-7},
+  urldate = {2025-07-01}
+}
+
 @article{Keller1980,
   title   = {Bubble Oscillations of Large Amplitude},
   author  = {Keller, Joseph B. and Miksis, Michael},

examples/run_all.sh

@@ -65,4 +65,7 @@ rm Gilmore_R1.000e-06_fa1.000e+06_pa1.000e+06.txt
 ./build/ultrasound_apecss -options oldroydb.apecss -tend 3e-6 -freq 3e6 -amp 400e3
 python3 plot_result_oldroydb.py
 rm RP_R1.000e-06_fa3.000e+06_pa4.000e+05.txt
+./build/ultrasound_apecss -options powerlaw.apecss -freq 636619 -amp 25331 -tend 25e-6
+python3 plot_result_powerlaw.py
+rm RP_R5.000e-06_fa6.366e+05_pa2.533e+04.txt
 cd ../

examples/ultrasound/README.md

@@ -9,6 +9,9 @@ This example simulated a lipid-coated microbubble using the Gilmore model and re
 #### Sonoluminescence (with emissions)
 Using [sonolumin_emissions.apess](sonolum_emissions.apecss) with ````./build/ultrasound_apecss -options sonolum_emissions.apecss -tend 40e-6 -freq 23.5e3 -amp 145e3````, this example reproduces the argon bubble studied by [Holzfuss, _Proc. R. Soc. A: Math. Phys. Eng. Sci._ 466 (2010), 1829](https://doi.org/10.1098/rspa.2009.0594) in the context of sonoluminesence. The acoustic emissions at _pLmax_ may be compared to Figure 5 of Holzfuss' work.
 
+#### Power-law liquid
+Using [powerlaw.apess](powerlaw.apecss) with ````./build/ultrasound_apecss -options powerlaw.apecss -freq 636619 -amp 25331 -tend 25e-6````, this example reproduces Figure 2c of [Kaykanat & Uguz, _The European Physical Journal Special Topics_ 233 (2024), 1625](https://doi.org/10.1140/epjs/s11734-024-01174-7).
+
 #### Kelvin-Voigt
 Using [kelvinvoigt.apess](kelvinvoigt.apecss) with ````./build/ultrasound_apecss -options kelvinvoigt.apecss -tend 6e-6 -freq 1e6 -amp 3e6````, this example reproduces Figure 5b of [Yang & Church, _Journal of the Acoustical Society of America_ 118 (2005), 3595](https://doi.org/10.1121/1.2118307).
 

examples/ultrasound/plot_result_powerlaw.py

@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+plt.rcParams['font.family']='serif'
+plt.rcParams['font.serif']=['Times New Roman'] + plt.rcParams['font.serif']
+plt.rcParams['mathtext.fontset']='stix'
+plt.rcParams['font.size']=10
+
+cm = 1/2.54
+
+Bubble = np.genfromtxt("RP_R5.000e-06_fa6.366e+05_pa2.533e+04.txt", delimiter=" ")
+
+fig1 = plt.figure(figsize=(17*cm,5*cm))
+ax1 = plt.subplot2grid((1,3),(0,0),colspan=1)
+ax2 = plt.subplot2grid((1,3),(0,1),colspan=1)
+ax3 = plt.subplot2grid((1,3),(0,2),colspan=1)
+plt.subplots_adjust(wspace=1.2*cm,hspace=1.2*cm)
+
+ax1.set(xlabel=r'$t\, f$',ylabel=r'$R(t)/R_0$')
+ax1.set_xlim(xmin=0,xmax=15)
+ax1.grid(color='gainsboro', linestyle='-', linewidth=0.5)
+ax1.plot(Bubble[:, 1]/1.570796e-6, Bubble[:, 3]/5e-6, linestyle='solid', linewidth=1,color='steelblue')
+
+ax2.set(xlabel=r'$t\, f$',ylabel=r'$\dot{R}(t)$[m/s]')
+ax2.set_xlim(xmin=0,xmax=15)
+ax2.grid(color='gainsboro', linestyle='-', linewidth=0.5)
+ax2.plot(Bubble[:, 1]/1.570796e-6, Bubble[:, 4], linestyle='solid', linewidth=1,color='steelblue')
+
+ax3.set_yscale('log')
+ax3.set(xlabel=r'$t\, f$',ylabel=r'$p_\mathrm{G}(t)$ [Pa]')
+ax3.set_xlim(xmin=0,xmax=15)
+ax3.grid(color='gainsboro', linestyle='-', linewidth=0.5)
+ax3.plot(Bubble[:, 1]/1.570796e-6, Bubble[:, 5]/1e3, linestyle='solid', linewidth=1.0,color='steelblue')
+
+ax1.xaxis.set_label_coords(0.5,-0.24)
+ax2.xaxis.set_label_coords(0.5,-0.24)
+ax3.xaxis.set_label_coords(0.5,-0.24)
+
+ax1.yaxis.set_label_coords(-0.28, 0.5)
+ax2.yaxis.set_label_coords(-0.25, 0.5)
+ax3.yaxis.set_label_coords(-0.25, 0.5)
+
+fig1.savefig('ultrasound_powerlaw.pdf', bbox_inches='tight',pad_inches=0.035)

examples/ultrasound/powerlaw.apecss

@@ -0,0 +1,34 @@
+#########################################################
+#                                                       #
+#  APECSS Options File                                  #
+#                                                       #
+#########################################################
+
+BUBBLE
+InitialRadius 5.0e-6
+RPModel RP
+PressureAmbient 1.01325e5
+END
+
+GAS
+EoS IG
+ReferencePressure 1.0e5
+ReferenceDensity 1.2
+PolytropicExponent 1.4
+END
+
+LIQUID
+LiquidType PowerLaw
+ReferencePressure 1.0e5
+ReferenceDensity 1000.0
+PowerLawExponent 0.75
+PowerLawConsistencyCoeff 0.002
+END
+
+INTERFACE
+SurfaceTensionCoeff 0.0729
+END
+
+RESULTS
+Bubble
+END

include/apecss.h

@@ -2,7 +2,7 @@
 // for the computation of pressure-driven bubble dynamics and acoustic
 // emissions in spherical symmetry.
 //
-// Copyright (C) 2022-2024 The APECSS Developers
+// Copyright (C) 2022-2025 The APECSS Developers
 //
 // The APECSS Developers are listed in the README.md file available in
 // the GitHub repository at https://github.com/polycfd/apecss.
@@ -102,6 +102,7 @@ typedef double APECSS_FLOAT;
 #define APECSS_LIQUID_KELVINVOIGT (2)  // Kelvin-Voigt solid
 #define APECSS_LIQUID_ZENER (3)  // Zener solid, requires two additional ODEs
 #define APECSS_LIQUID_OLDROYDB (4)  // Oldroyd-B liquid, requires two additional ODEs
+#define APECSS_LIQUID_POWERLAW (5)  // Power-law liquid
 
 // Lipid coating model (bit-wise)
 #define APECSS_LIPIDCOATING_NONE (1)  // Clean interface
@@ -191,6 +192,9 @@ struct APECSS_Liquid
   APECSS_FLOAT G;  // Shear modulus [Pa]
   APECSS_FLOAT eta;  // Polymer viscosity [Pa s]
   APECSS_FLOAT lambda;  // Relaxation time [s]
+  APECSS_FLOAT n;  // Power-law exponent
+  APECSS_FLOAT k;  // Power-law consistency index [Pa s^n]
+  APECSS_FLOAT k_ext;  // Power-law extended consistency index [Pa s^n]
 
   // Pointers to the functions describing the properties of the liquid
   APECSS_FLOAT (*get_density)(APECSS_FLOAT p, struct APECSS_Liquid *Liquid);
@@ -622,8 +626,11 @@ APECSS_FLOAT apecss_liquid_pressurederivative_bubblewall_explkelvinvoigt(APECSS_
 APECSS_FLOAT apecss_liquid_pressurederivative_bubblewall_zener(APECSS_FLOAT *Sol, APECSS_FLOAT t, struct APECSS_Bubble *Bubble);
 APECSS_FLOAT apecss_liquid_pressurederivative_bubblewall_exploldroydb(APECSS_FLOAT *Sol, APECSS_FLOAT t, struct APECSS_Bubble *Bubble);
 APECSS_FLOAT apecss_liquid_pressure_viscous(APECSS_FLOAT R, APECSS_FLOAT U, struct APECSS_Bubble *Bubble);
+APECSS_FLOAT apecss_liquid_pressure_viscous_powerlaw(APECSS_FLOAT R, APECSS_FLOAT U, struct APECSS_Bubble *Bubble);
 APECSS_FLOAT apecss_liquid_pressurederivative_viscous_expl(APECSS_FLOAT R, APECSS_FLOAT U, struct APECSS_Bubble *Bubble);
 APECSS_FLOAT apecss_liquid_pressurederivative_viscous_impl(APECSS_FLOAT R, struct APECSS_Bubble *Bubble);
+APECSS_FLOAT apecss_liquid_pressurederivative_viscous_powerlaw_expl(APECSS_FLOAT R, APECSS_FLOAT U, struct APECSS_Bubble *Bubble);
+APECSS_FLOAT apecss_liquid_pressurederivative_viscous_powerlaw_impl(APECSS_FLOAT R, struct APECSS_Bubble *Bubble);
 APECSS_FLOAT apecss_liquid_pressurederivative_viscous_nonimpl(APECSS_FLOAT R, struct APECSS_Bubble *Bubble);
 
 // ---------------------