Dynamics and long-time behavior of a small bubble in viscous liquids with applications to food rheology

Abstract

In this work, the physical behavior of a small, spherical bubble in a highly viscous, incompressible liquid phase is analyzed under specific pressure impacts. This represents an attractive topic in the food industry, since it is of interest to know under which conditions this two-phase dispersion exhibits a stable state. The specific material law of a second-order liquid, which includes Newtonian and non-Newtonian material constants, provides a nonlinear initial value problem for the radius of the bubble. This system is solved numerically by an efficient version of the classical Runge–Kutta method. By parameter variation, the impact of the dimensionless quantities associated with inertia, non-Newtonian material coefficients, pressure, surface tension and viscosity on the two-phase system is investigated. This particularly yields insights into the stability behavior of the bubble surface. The solution curves show various characteristics such as asymptotic oscillations or monotonically decreasing profiles. These results are transferred to a specific non-Newtonian and Newtonian substance. Finally, by studying stationary solutions, it becomes obvious that only the excitation pressure and the surface tension determine the new equilibrium state of the bubble, which in particular represents its long-time behavior. Furthermore, a sinusoidal driving pressure is used to investigate unstable solutions. The aim of the paper was to bring together these mathematical stability results to practice-oriented considerations.

Appendices

Appendix 1

Remark 1

The complete solution of (19) is obtained from the linearized IVP

$$\begin{aligned} \ddot{y}+2\alpha \dot{y}+\beta y=f(\bar{t}), \end{aligned}$$

(47)

$$\begin{aligned} y(0)=\dot{y}(0)=0, \end{aligned}$$

(48)

$$\begin{aligned} \hbox {with eigenvalues }\lambda _{1,2}=-\alpha \pm \sqrt{\alpha ^{2}-\beta }, \end{aligned}$$

(49)

where for the case $\lambda _{1}\ne \lambda _{2}$ one has

$$\begin{aligned} y(\bar{t})=\mathop {\int }\limits _{0}^{\bar{t}} {\frac{\mathrm{e}^{\lambda _{1}\left( \bar{t}-\bar{s} \right) }-\mathrm{e}^{\lambda _{2}\left( \bar{t}-\bar{s} \right) }}{\lambda _{1}-\lambda _{2}}f(\bar{s})\mathrm{d}\bar{s}} =\mathop {\int }\limits _{0}^{\bar{t}} {\frac{\mathrm{e}^{\lambda _{1}\left( \bar{t}-\bar{s} \right) }-\mathrm{e}^{\lambda _{2}\left( \bar{t}-\bar{s} \right) }}{2\sqrt{\alpha ^{2}-\beta } }f(\bar{s})\mathrm{d}\bar{s}}. \end{aligned}$$

(50)

For $\lambda _{1}=\lambda _{2}$, one has to take the limit $\lambda _{1}\longrightarrow \lambda _{2}=-\alpha$. In particular, it is necessary to show that for a bounded right-hand side f in (47), the limiting case $\bar{t}\longrightarrow \infty$ can be performed. Let $\left| f(\bar{s}) \right| \le K(\bar{t}):=\max \left\{ \left| f(\bar{s}) \right| :0\le \bar{s}\le \bar{t} \right\}$. Then, due to ${0>\lambda }_{1}>\lambda _{2}$:

$$\begin{aligned} \left| y(\bar{t}) \right|\le & {} \mathop {\int }\limits _{0}^{\bar{t}} {\frac{\mathrm{e}^{\lambda _{1}\left( \bar{t}-\bar{s} \right) }-\mathrm{e}^{\lambda _{2}\left( \bar{t}-\bar{s} \right) }}{\lambda _{1}-\lambda _{2}}\left| f(\bar{s}) \right| \mathrm{d}\bar{s}} \le K(\bar{t})\mathop {\int }\limits _{0}^{\bar{t}} {\frac{\mathrm{e}^{\lambda _{1}\left( \bar{t}-\bar{s} \right) }-\mathrm{e}^{\lambda _{2}\left( \bar{t}-\bar{s} \right) }}{\lambda _{1}-\lambda _{2}}\mathrm{d}\bar{s}}\nonumber \\\le & {} \frac{K(\bar{t})}{\lambda _{1}\cdot \lambda _{2}}\left( 1-\frac{\mathrm{e}^{\lambda _{2}\bar{t}}\lambda _{1}-\mathrm{e}^{\lambda _{1}\bar{t}}\lambda _{2}}{\lambda _{1}-\lambda _{2}} \right) \le \frac{K(\bar{t})}{\beta } \quad \hbox {for}\, \lambda _{1}\ne \lambda _{2}. \end{aligned}$$

(51)

This is also valid for $\lambda _{1}=\lambda _{2}$. These estimates can be similarly made for $\dot{y}(\bar{t})$. Hence, these considerations justify the approximation because for small K and not too small $\beta$, the approximation and its derivatives remain small. These assumptions yield

$$\begin{aligned} \left| y(\bar{t}) \right| \le \frac{K(\bar{t})}{4\sigma ^{*}+3}<\frac{K(\bar{t})}{3} \quad \hbox { with }\quad K(\bar{t}):=\max \left\{ \left| 1-\bar{p}_{\infty }(\bar{s}) \right| :0\le \bar{s}\le \bar{t} \right\} . \end{aligned}$$

(52)

Hence, also here the approximation is reasonable. In the stationary case with ${1-\bar{p}}_{\infty }=:K$, one obtains

$$\begin{aligned} y<\frac{K}{4\sigma ^{*}+3}. \end{aligned}$$

(53)

Proposition 2

For the compression of an incompressible liquid between two parallel plates moving toward each other with a constant strain rate (Fig. 20), the pressure function is given by

$$\begin{aligned} \bar{p}_{\infty }(\bar{t})=C_{1}-C_{2}\,{\mathrm{exp}}\left( 2\dot{\bar{\varepsilon }}_\mathrm{c}\bar{t} \right) , \end{aligned}$$

(54)

where $C_{1}\in \mathbb {R}$ and $C_{2}\in [0,\infty [.$

Fig. 20

Compression of a bubble in a liquid medium between two parallel plates which approach each other with a constant strain rate

Proof

For the derivation of the (dimensionless) pressure function used in the numerical analysis, a compression procedure as in Fig. 20 is considered. Let D denote the distance between the bubble and the upper or lower plate. Moreover, let d stand for the bubble diameter. Their ratio is $D{/}d=500\gg 1$. The same holds for the ratio consisting of the distance between the bubble and the side walls and d. This ensures that the elongational uniaxial compression in Fig. 20 leads to a spherical compression of the centrally placed gas bubble. Assuming that the parallel plates with distance $\bar{L}=\bar{L}(\bar{t})$ approach each other with constant strain rate ${\partial \bar{u}} / {\partial \bar{x}}=\dot{\bar{\varepsilon }}_\mathrm{c}$ leads to $\bar{u}( \bar{x})=\dot{\bar{\varepsilon }}_\mathrm{c}\bar{x}+c_{1}$ with $c_{1}\in \mathbb {R}$. The continuity equation for an incompressible liquid then yields

$$\begin{aligned} \frac{\partial \bar{u}}{\partial \bar{x}}+\frac{\partial \bar{v}}{\partial \bar{y}}=0 \Rightarrow \frac{\partial \bar{v}}{\partial \bar{y}}=-\dot{\bar{\varepsilon }}_\mathrm{c} \Rightarrow \bar{v}\left( \bar{y} \right) =-\dot{\bar{\varepsilon }}_\mathrm{c}\bar{y}+c_{2},\quad c_{2}\in \mathbb {R}. \end{aligned}$$

(55)

The $\bar{x}$-component of the stationary Navier–Stokes equations without body forces has the form

$$\begin{aligned} \bar{u}\frac{\partial \bar{u}}{\partial \bar{x}}+\bar{v}\frac{\partial \bar{u}}{\partial \bar{y}}=-\frac{1}{\rho _\mathrm{l}}\frac{\partial \bar{p}}{\partial \bar{x}}+\nu \left( \frac{\partial ^{2}\bar{u}}{\partial \bar{x}^{2}}+\frac{\partial ^{2}\bar{v}}{\partial \bar{y}^{2}} \right) \Rightarrow \frac{\partial \bar{p}}{\partial \bar{x}}=-\rho _\mathrm{l}\dot{\bar{\varepsilon }}_\mathrm{c}\left( \dot{\bar{\varepsilon }}_\mathrm{c}\bar{x}+c_{1} \right) . \end{aligned}$$

(56)

Then, from the expression for $\bar{u}$, it follows that

$$\begin{aligned} \bar{u}=\frac{\partial \bar{L}}{\partial \bar{t}} \Rightarrow \frac{\partial \bar{L}}{\partial \bar{t}}=\dot{\bar{\varepsilon }}_\mathrm{c}\bar{L}+c_{1} \Rightarrow \bar{L}(\bar{t})=c \hbox {exp}\left( \dot{\bar{\varepsilon }}_\mathrm{c}\bar{t} \right) -\frac{c_{1}}{\dot{\bar{\varepsilon }}_\mathrm{c}} \end{aligned}$$

(57)

and finally by (55) that

$$\begin{aligned} \frac{\partial \bar{p}}{\partial \bar{L}}=-\dot{\bar{\varepsilon }}_\mathrm{c}^{2}\bar{L}-\dot{\bar{\varepsilon }}_\mathrm{c}c_{1}\Rightarrow & {} \bar{p}(\bar{t})=-\rho _\mathrm{l}\left( \frac{\left( \dot{\bar{\varepsilon }}_\mathrm{c}c \right) ^{2}}{2}\hbox {exp}\left( 2\dot{\bar{\varepsilon }}_\mathrm{c}\bar{t} \right) -\frac{c_{1}^{2}}{2}+c_{2} \right) \end{aligned}$$

(58)

$$\begin{aligned}\Rightarrow & {} C_{1}:=\rho _\mathrm{l}\left( \frac{c_{1}^{2}}{2}-c_{2} \right) ,\quad C_{2}:=-\frac{\rho _\mathrm{l}}{2}\left( \dot{\bar{\varepsilon }}_\mathrm{c}c \right) ^{2}. \end{aligned}$$

(59)

$\square$

Remark 3

The values of $\rho _\mathrm{l}, \nu$ and $\sigma$ for thickened whole milk in Table 1 were measured in the laboratories of the Institute of Fluid Mechanics (LSTM) at FAU Erlangen-Nuremberg, where the glass pycnometer and the ring method were used to determine $\rho _\mathrm{l}$ and $\sigma$, respectively. For the measurement of the normal stress coefficients $\psi _{1}$ and $\psi _{2}$, a rotational rheometer was used. Since slowly varying motions play a central role, one can assume a hydrostatic pressure for $p_{0}$ and a height of the liquid layer of 1 m, for mathematical convenience (see Proposition 2). Note that two dimensionless numbers in Table 1 on the right-hand side, namely $\eta ^{*}$ and $M^{*}$, are identical because their definitions in (14) and that of the characteristic time lead to

$$\begin{aligned} \eta ^{*}=\frac{\eta }{p_{0}\tau }=\frac{\eta R_{0}^{2}/ \nu }{p_{0}\tau ^{2}}=\frac{R_{0}^{2}\rho _\mathrm{l}}{p_{0}\tau ^{2}}=M^{*}. \end{aligned}$$

(60)

Table 1 Physical quantities and dimensionless parameters for the mixture of sodium alginate with whole milk related to a bubble with initial radius $R_{0}=1\,\mathrm{mm}$

Table 2 Physical quantities from [34] (for a pH value of 6.8 which is approximately the value for milk) and dimensionless parameters for low $\upbeta$-lactoglobulin percentage in distilled water related to a bubble with initial radius $R_{0}=1\,\mathrm{mm}$

Fig. 21

Curve profiles of $Q( k_\mathrm{o})$ for different values of $p_\mathrm{c}$ and $\sigma ^{*}$ related to (32)–(34)

Proposition 4

Properties of the third-order polynomial

$$\begin{aligned} Q( k_\mathrm{o})=p_\mathrm{c}k_\mathrm{o}^{3}+2\sigma ^{*}\left( k_\mathrm{o}^{2}-1 \right) -1. \end{aligned}$$

(61)

• (i) Equation (61) implies

  $$\begin{aligned} Q(0)= & {} -2\sigma ^{*}-1<0 \end{aligned}$$

  (62)
  $$\begin{aligned} Q'( k_\mathrm{o})= & {} 3p_\mathrm{c}k_\mathrm{o}^{2}+4\sigma ^{*}k_\mathrm{o}>0 \end{aligned}$$

  (63)
  $$\begin{aligned} Q( k_\mathrm{o})\longrightarrow & {} \infty \quad \hbox {if}\quad k_\mathrm{o}\longrightarrow \infty \end{aligned}$$

  (64)
  $$\begin{aligned} Q\left( 1 \right)= & {} p_\mathrm{c}-1 \end{aligned}$$

  (65)
  $$\begin{aligned} 0<p_\mathrm{c}<1\Leftrightarrow & {} Q\left( 1 \right) <0 \end{aligned}$$

  (66)
  $$\begin{aligned} p_\mathrm{c}>1\Leftrightarrow & {} Q\left( 1 \right) >0 \end{aligned}$$

  (67)

• (ii) From (62), (63) and (64), it follows for $Q( k_\mathrm{o})=0$:

  There is one unique positive zero $k_\mathrm{o}>0$ of Q.

• (iii) The relations (66) and (67) yield

  $$\begin{aligned} 0<p_\mathrm{c}<1\Leftrightarrow & {} k_\mathrm{o}>1 \end{aligned}$$

  (68)
  $$\begin{aligned} p_\mathrm{c}>1\Leftrightarrow & {} 0<k_\mathrm{o}<1 \end{aligned}$$

  (69)
  $$\begin{aligned} p_\mathrm{c}=1\Leftrightarrow & {} k_\mathrm{o}=1. \end{aligned}$$

  (70)

Proof

• (i) The verification of (62)–(67) requires simple computations.

• (ii) Due to (62) and (64), there is a change of sign from minus to plus. The uniqueness follows from the fact that Q is strictly monotonically increasing, see (63).

• (iii) To see this, let $T\,{:}\,k_\mathrm{o}\longrightarrow p_\mathrm{c}$ be the map

  $$\begin{aligned} p_\mathrm{c}=T( k_\mathrm{o})=\left( 1+2\sigma ^{*} \right) \frac{1}{k_\mathrm{o}^{3}}-2\sigma ^{*}\frac{1}{k_\mathrm{o}}=\frac{1}{k_\mathrm{o}^{3}}\left( \left( 1+2\sigma ^{*} \right) -2\sigma ^{*}k_\mathrm{o}^{2} \right) ({\ge }0). \end{aligned}$$

  (71)

  One directly concludes the equivalences (68)–(70) by rewriting (71) as

  $$\begin{aligned} p_\mathrm{c}=T( k_\mathrm{o})=\frac{1}{k_\mathrm{o}^{3}}\left( 1+2\sigma ^{*}\left( 1-k_\mathrm{o}^{2} \right) \right) . \end{aligned}$$

  (72)

Hence,

$$\begin{aligned}&p_\mathrm{c}>1 \quad \mathrm{if} \quad 0<k_\mathrm{o}<1 \end{aligned}$$

(73)

$$\begin{aligned}&p_\mathrm{c}=1 \quad \mathrm{if} \quad k_\mathrm{o}=1 \end{aligned}$$

(74)

$$\begin{aligned}&p_\mathrm{c}<1 \quad \mathrm{if} \quad k_\mathrm{o}>1 \end{aligned}$$

(75)

which completes the proof. $\square$

Remark 5

Because of (31), all curves must begin at a negative value of Q(0), then increase continuously to infinity, which is obviously the case in Fig. 21. If $p_\mathrm{c}=1$, then a variation of the value of $\sigma ^{*}$ effectuates a change in value of Q(0), but the zero remains unchanged at $k_\mathrm{o}=0.$ Apart from this, the zero varies in the ranges indicated by (32)–(34).

Appendix 2: Symbols and notations

1.1 Coordinates and functions

$\left( r,\vartheta ,\varphi \right) ^\mathrm{T}$ :

Spherical coordinates

$s, t, \left( t_{\mathrm{e}_{R_{0}}}, t_{\mathrm{j}_{R_{0}}} \right)$ :

(Specific) time variables

${\mathbf {v}}=\left( v(r,t),0,0 \right) ^\mathrm{T}$ :

Velocity vector field in spherical coordinates

$R=R(t), \left( R_{0}=R(0) \right)$ :

Bubble radius (at rest)

d :

Bubble diameter

D :

Distance between bubble and plate

$\tau$ :

Characteristic excitation time

$c_{1}, c_{2},c_{3},c_{4},C,K_{1},K_{2}, \chi$ :

Constant functions

$\mathrm{e}^{t}, \mathrm{exp}(t)$ :

Exponential function

${\mathcal {O}}(.)$ :

Big-O notation

1.2 Physical and rheological quantities

$p, \tilde{p}$ :

Isotropic pressure

$p_\mathrm{g}, p$ :

Gas and liquid pressure

$p_{\mathrm{g}_{0}}=p_\mathrm{g}(0), p_{0}=p(0)$ :

Gas and liquid pressure at rest

$p_{\infty }=p_{\infty }(t), p_\mathrm{c}$ :

Far-field pressure, constant pressure

$V, \left( V_{0} \right)$ :

Volume (at rest)

$\dot{\varepsilon }_\mathrm{c}, \left( \dot{\varepsilon }_{0} \right)$ :

(Reference) rate of compression, strain rate

$f, (\omega )$ :

(Angular) frequency

A :

Amplitude of sinusoidal pressure excitation

$g=9.81\,\mathrm{m/s}^{2}$ :

Gravitational acceleration

$N_{i}$ :

Normal stress differences

$\dot{\gamma }$ :

Shear rate

$\rho _\mathrm{g}, \rho _\mathrm{l}, \varDelta \rho :=\rho _\mathrm{g}-\rho _\mathrm{l}$ :

Gas and liquid density, density difference

$\delta _{i,j}$ :

Kronecker delta

$\lambda$ :

Bulk viscosity

$\lambda _{1}, \lambda _{2}$ :

Eigenvalues

$\eta , \left( \eta _\mathrm{l} \right)$ :

(Liquid) dynamic viscosity

$\nu , \left( \nu _\mathrm{l} \right)$ :

(Liquid) kinematic viscosity

$\psi _{1}, \psi _{2}, \left( \psi _{10}, \psi _{20} \right)$ :

(Limits of) normal stress coefficients

$\alpha , \beta$ :

Non-Newtonian transport coefficients

$\sigma$ :

Surface tension

1.3 Dimensionless numbers and variables

$Bo$ :

Bond number

$Ga$ :

Galileo number

$Re$ :

Reynolds number

$A^{*}, B^{*}$ :

Non-Newtonian transport coefficients

$M^{*}$ :

Inertia constant

$\sigma ^{*}$ :

Surface tension constant

$\eta ^{*}$ :

Viscosity constant

$p_{\mathrm{g}_{0}}^{*}$ :

Gas pressure at rest

N :

Combination of $B^{*}, k_\mathrm{o}$ and $M^{*}$

${\mathfrak {R}}_\mathrm{o}$ :

Combination of $N, k_\mathrm{o}, \eta ^{*}, p_\mathrm{c}$ and $\sigma ^{*}$

$\bar{s}, \bar{t}$ :

Time variables

$\bar{R}=\bar{R}(\bar{t}), a=a(\bar{t})$ :

Bubble radius, linearized bubble radius

$k_\mathrm{o}$ :

Stationary solution

$k_{\mathrm{o},\mathrm{e},\mathrm{s}}, k_{\mathrm{o},\mathrm{j},\mathrm{s}}, k_{\mathrm{o},\mathrm{twm}}, k_{\mathrm{o},\mathrm{blg}}$ :

Specific stationary solutions

$Q=Q( k_\mathrm{o})$ :

Cubic polynomial

$f=f(\bar{t}), K=K(\bar{t}),y=y(\bar{t})$ :

Scalar functions

$\bar{p}=\bar{p}(\bar{t}), \bar{p}_{\mathrm{c},\infty }, \bar{p}_{0}$ :

Time-dependent pressure, constant pressures

$\dot{\bar{\varepsilon }}_\mathrm{c}$ :

Rate of compression

$\bar{\omega }$ :

Angular frequency

1.4 Differential operators and tensors

$\dot{R}:=\frac{\partial R}{\partial t}, \ddot{R}:=\frac{\partial ^{2}R}{\partial t^{2}}$ :

Single and double time derivative of bubble radius

$\mathrm{div} \, {\mathbf {f}}=\nabla \cdot \left( f_{1}, f_{2}, f_{3} \right) ^\mathrm{T} = \frac{\partial f_{1}}{\partial x_{1}}+\frac{\partial f_{2}}{\partial x_{2}}+\frac{\partial f_{3}}{\partial x_{3}}$ :

Divergence

${\mathbf {I}}:=\varvec{\delta }_{ij}=\mathrm{diag}\left( 1, 1, 1 \right)$ :

Unit tensor

${\mathbf {A}}_{1}, {\mathbf {A}}_{2}$ :

Rivlin–Ericksen tensors

$\varvec{\sigma }, \sigma _{rr}, \sigma _{\vartheta \vartheta }, \sigma _{\varphi \varphi }$ :

Stress tensor, diagonal components in spherical coordinates

$\varvec{\tau }$ :

Viscous stress tensor

Abbreviations

BLG:

$\upbeta$-Lactoglobulin

IVP:

Initial value problem

ODE:

Ordinary differential equation

RPE:

Rayleigh–Plesset equation

RKM:

Runge–Kutta method

RKS:

Runge–Kutta solution