Non-Newtonian coefficient condition for a stable long-time behavior of a single bubble: existence and characteristics of stable solutions

The present paper takes up the underlying nonlinear initial value problem from a preceding author’s work about the dynamics of a single bubble in a highly viscous liquid medium under different pressure impacts. The arising ordinary differential equation is mainly based on the constitutive relation of a second-order liquid that in particular includes two non-Newtonian material constants. In this article, the significance of these coefficients is mathematically analyzed in detail by proving the existence of stable solutions of the named initial value problem. This is achieved by special transformations of the differential equation at hand and the introduction of appropriate Lyapunov functions. It particularly turns out that a combined condition of the non-Newtonian coefficients and diverse restrictions to the external pressure impact are decisive for the validity of the existence results. Furthermore, the convergence speed of solutions is investigated by considering the linearized equation associated with the present initial value problem and by applying a special variant of Gronwall’s lemma. The main theoretical result, being the prementioned strong condition for the non-Newtonian coefficients, is finally compared to real data sets.

of unit cells of foams during storage and transport are necessary to maintain specific positive organoleptic characteristics and a certain shelf life of a food product [1][2][3]. But bubbles and foams can also be obstacles for production cycles which is currently considered in [4]. The task of the mathematical analysis is now to determine ranges of relevant physical and rheological parameters in which stability of bubbles can be expected.
In [5], a two-phase system consisting of a small, spherical gas bubble in a highly viscous, incompressible liquid was theoretically examined by the conservation equations for mass and momentum at changing pressure conditions. The important point to note here is that the material law of a second-order liquid was applied for the stress tensor to take into account slow and slowly varying bubble motions [6][7][8][9]. The resulting nonlinear initial value problem (IVP) was numerically solved by a Runge-Kutta method of order four. Performing (dimensionless) parameter studies of the occurring physical and rheological quantities showed that the pressure excitation and the surface tension represent the strongest influencing factors. They determine the new equilibrium state of the bubble which simultaneously represents its long-time behavior. Furthermore, an increasing liquid viscosity effects considerably damped bubble surface oscillations. The constitutive equation of the second-order fluid includes two non-Newtonian material coefficients which may also have bubble stabilizing effects. The one connected to the nonlinear character of the fluid causes at increasing values a damping of bubble surface vibrations, while that one connected to the memory effects leads to an increment in surface oscillations almost similar to the effects of inertia. Deriving an analytical solution by linearization of the named IVP from [5] yields an explicit expression for the bubble radius taking into account only one of the non-Newtonian material constants. For stationary solutions of the original ordinary differential equation (ODE), an upper bound could be derived only depending on surface tension.
The main objective of the present paper is to provide a solid mathematical basis for a certain class of practical applications for which stability of single bubbles in high-viscosity liquids can be expected under specific assumptions. In this context, the aim is to extract more information from the above-named nonlinear IVP by transforming the underlying differential equation into a well-structured, mathematically handable ODE. This is necessary for a conclusive, mathematical argumentation to find solutions under more general conditions, in particular involving the aforementioned non-Newtonian material coefficients. Regarding the scientific literature, it becomes apparent that there is very rare work related to the mathematical significance of the non-Newtonian contribution of a second-order liquid in conjunction with bubble dynamics. A theoretical model for a corresponding two-phase system in connection with the rheology of incompressible and compressible non-Newtonian fluids is examined in [10,11], where slow fluid motions play a central role. Related bubble analyses are carried out by the mentioned second-order approximation of the stress tensor, where the main attention is devoted to the bulk viscosity. A phenomenological description of the fluid-mechanical transport coefficients for the compressible case is thoroughly performed in [10]. In the incompressible case, the author obtains a special variant of the Rayleigh-Plesset equation (RPE) (comparable to that in [5]), but only the linearization including just one of the non-Newtonian coefficients is taken into consideration. In [12], the second-order constitutive relation is used for incompressible fluids in plane creeping flows. A commonly used assumption is that the sum of both coefficients is zero [13,14]. Moreover, experimental measurements reveal that the coefficients differ from each other only by a negative sign [9]. The authors from [14] investigate steady motions of second-order and even third-grade fluids in aperture domains by analyzing existence and regularity of (weak) solutions. Also here, in the case of third-order fluids, standard assumptions are applied, see also [15].
In the following section, the derivation of the prementioned IVP is briefly taken up to initiate the transformation process. For this purpose, appropriate functions are introduced to define, among others, integral expressions. The latter help to derive a nonlinear first-order ODE system so that an appropriate Lyapunov function can be introduced. Its properties are crucial for the validity of the convergence of stable solutions. Section 3 considers this aspect over long periods of time. Next, in Sect. 4, the ODE system is modified to formulate the main result of this work, namely an existence theorem including a strongly restricting condition for the two non-Newtonian coefficients of the second-order constitutive relation as well as specific boundedness conditions for the far-field pressure. Furthermore, it is mathematically shown by linearization of the original differential equation and by application of a special variant of Gronwall's lemma that the solutions and therefore also the physical process come to a stable rest state. The last section refers to the connection between the gained strong coefficient condition for the non-Newtonian contributions and application-oriented data sets.

Preparations for transformation of IVP
The following mathematical analysis is based on the consideration of a nonlinear ODE derived from the model of a two-phase system consisting of a small, single gas-filled bubble immersed in an incompressible, highly viscous liquid under different outer pressure impacts. The bubble is assumed to be spherically symmetric and the origin of the coordinate system is placed in the center of the bubble. Combination of the corresponding continuity and momentum equation connected with the condition of incompressibility, i.e., ∂ ∂t (r 2 v(r, t)) = 0, and the kinematic boundary condition on the bubble surface, i.e., v(R, t) =Ṙ, yields a reduced equation of motion for the bubble radius R = R(t); see [5] for more details. The derived dimensionless ODE is then given byR where R, p and t represent the bubble radius, the pressure and the time variable, respectively. The index 0 denotes the corresponding values at rest, that are R 0 = R(0) and p 0 = p(0), and τ stands for the excitation time τ := R 2 0 /ν l with the liquid kinematic viscosity ν l . Moreover, the non-dimensional coefficients are defined by Here α and β are non-Newtonian material coefficients arising from the constitutive relation of a second-order liquid [6][7][8][9], namely with the Rivlin-Ericksen tensors A 1 and A 2 . The remaining quantities are the liquid density l , the surface tension σ , the dynamic viscosity η, and the gas pressure at rest p g 0 . For steady solutions assumeR and letp Then, (1) becomes a third-degree polynomial Equations (1) and (7) are rewritten as To qualitatively analyze the ODE (8), one has to apply appropriate substitutions. Using For the sake of simplification, one defines for Thus (9) becomesẍ As preparation for the subsequent theorems, one considers some help functions. For an appropriate constant c 0 > 0, the expression is computed by using partial fraction decomposition and integration: The term (16) has the following behavior: where O(.) is the Big-O notation. Moreover, By defining one deduces For F, it is easy to check that In what follows, the subsequent behavior will be useful:

Convergence behavior for long times
On the basis of these preparations, the first existence theorem for stable solutions will be formulated [16]: Then for the solution x of (14) holds and then for the solution of (8) Proof Consider the following equivalent system to (14) x The respective calculations can be found in Appendix, Remark 2. As positive definite help function (referring to Lyapunov functions [17]) one considers is strictly monotonically decreasing. Therefore one finds x(t), z(t) for allt ≥ 0 and This implies the boundedness of x(t) and z(t). Furthermore, If now lim supt →∞ |x(t)| > 0 would hold, then one could find, sinceẋ is bounded, an infinite sequence of t−intervals of identical length above which |x(t)| > ε > 0. But the positivity of the integrand in (26) furnishes a contradiction. Hence, the convergence of the integral (26), with the boundedness of the derivative, brings along: limt →∞ x(t) = 0.

Preparatory work
Our next objective is to perform some preliminary calculations to ODE (23) to facilitate the derivation of the main result of the article given by Theorem 2. Using (22) leads to and, due to (7), the numerator becomes A complete calculation of (28) can be found in Appendix, Remark 3.
Set e(t) := −( p c −p ∞ (t)) and get Concerning this, it will be shown that P(0) = e(t)k 3 o and P − e(t) A more detailed calculation of the chain of inequalities (30) is given in Appendix, Remark 4. With (27) and an alternative form of (23) is given bẏ

Stability condition
The preceding mathematical preparations form the basis for the main result of this paper, namely the convergence of the solution of Eqs. (14) and (8), respectively, under special assumptions to the non-Newtonian material coefficients and the external pressure excitation. On the one hand, a strong relation between the dimensionless non-Newtonian coefficients A * and B * , which may be considered as physically reasonable, is required. Furthermore, boundedness conditions for the pressure excitationp ∞ are necessary for the validity of the subsequent theorem. Step 1 If x(t) ≥ K 0 with a constant K , then, due to (32), it follows thatż(t) < 0, hence z(t) strictly monotonically decreases. This impliesẋ(t) < 0 and x(t) also decreases. Thus, x(t) → ∞ can be excluded. (29)] holdsż(t) > 0 and z(t) strictly monotonically increases, which impliesẋ(t) > 0 so that x(t) also increases. Consequently, x(t) + k o → 0 is impossible. This ensures the existence of the functions (x(t), z(t)) for all timest > 0.
Step 2 For large timest with e(t) := p c −p ∞ (t) the Lyapunov function The solution (x(t), z(t)) is inserted. Then, differentiation with respect tot yields According to (15), (21), (29), and (34) holds and, withṗ ∞ (t) ≤ 0, one has fort ≥ T o The last relation results from the boundedness for − k o < x < 0 with A * 4B * < 1, and from On the basis of (37) and (38), by integration from T o tot, one derives for W and fort → ∞ the estimate which shows the importance of assumption (33) and immediately entails the boundedness of x(t), z(t) and then also ofẋ(t).
Step 3 The same argumentation as in the proof of Theorem 1 is used. The solution (x(t), z(t)) of (23) with and the total derivative is calculated as Integration gives

s)) h (x(s)) z(s) e(s)ds. (41)
The boundedness of the functions and the assumption (34) ensure the convergence of the second integral and also of

Convergence behavior
Theorem 2 above states that the solution of (14) quantitatively tends to zero. Under this assumption the final subsection gives a result for the convergence speed of the solution of Theorem 2 in proximity to zero. The solution converges exponentially and in accordance with |( p c −p ∞ (t))|. For this purpose, the linearization of Eq. (1) by Taylor expansion, see [5], is used, that is The reason why one may also consider this linearized ODE system is due to a stability theorem in [18] which allows the transfer of the knowledge about the stability behavior of the solution of (42) to that of the original nonlinear problem (1). Equation (42) is transformed into an ODE system of first order by letting y = (y 1 , y 2 ) T with y 1 =R − k o and y 2 =Ṙ so thatẏ 1 =Ṙ andẏ 2 =R. It followṡ with the 2 × 2 matrix A and its components From (44) and one obtains negative eigenvalues of A which fulfills the assumption of the next theorem.
In the subsequent result, for the sake of convenience, the parametert is replaced by t and the bold letters from (43) except A are rewritten by normal letters.

Theorem 3 Let
A be a matrix with eigenvalues λ j with negative real parts, and consider If lim t→∞g (t) = 0, then for the solution of (47) holds A variant of Gronwall's lemma [18] is needed for the proof of Theorem 3.
Then the following inequality holds The proof of this lemma can be found in Appendix, Remark 5.
Proof of Theorem 3 At first, ODE (47) endowed with y(0) = y 0 is transformed into an integral equation:
Due to (52) the limiting process lim t→∞ y p (t) = 0 holds. This can also be applied on the above integral because Hence, where the convergence speed for large times t approximatively behaves as ce − λ 2 t +g(t). Table 1 The experimental data for the non-Newtonian coefficients α and β for thickened whole milk are taken from [5], those for polyisobutylene and osteoarthritic fluid come from [20] and for silicone oil from [21] Material The corresponding dimensionless parameters A * and B * are given as well

Connection between theory and real parameters
In this section, the strong dimensionless condition from the assumptions of Theorem 2 in Sect. 4, given by is compared with concrete experimental data of diverse viscoelastic materials from food processing [5], bioengineering [19] and technical applications [20,21]. For this purpose, real measurement values of the non-Newtonian contributions α and β from (4) are necessary to obtain results for the dimensionless ratio of the left-hand side of inequality (54). More precisely, the definitions of A * and B * from (3) associated with the nonlinear ODE (8) [or (9)] are applied for the computation.
Remark 1 Transforming (54) by means of (3) into its dimensional form yields the equivalent condition: Hence, the non-dimensional relation only depends on the non-Newtonian contributions α and β because liquid pressure at rest p 0 and characteristic excitation time τ vanish during the calculation process.
Regarding the domain of food industries, there is very few literature available about experimental measurements of normal stress coefficients. Therefore, only tabular values from [5] specifically for a medium composed of 2% sodium alginate in whole milk (3.5% fat) are used. Inserting the associated dimensionless values A * and B * from Table 1 into the left-hand side of (54) yields This value obviously fulfills condition (54) and is consistent with the mathematical statement from Theorem 2. For further application-oriented comparisons one can consider specific polymers from chemical engineering. The authors from [20] used corresponding non-Newtonian material parameters for the organic polymer polyisobutylene (in cetane) and for a biological osteoarthritic synovial fluid. According to Table 1, one obtains for polyisobutylene and for the osteoarthritic synovial fluid A * 4B * ≈ −1.5 < 1.
From [21] one derives for silicone oil: Clearly, all computed values for the considered non-Newtonian liquids fulfill condition (54). This is mainly due to the negativity of the coefficient β which leads to a negative calculation value in (56)-(59). The present practical-oriented result of (56) confirms the gained knowledge from [5] that there is negligible contribution of the non-Newtonian coefficients α and β to the bubble behavior. For the sake of completeness, it should be noted that the assumptions (33) and (34) from Theorem 2 always hold in the long term for the applied pressure functions, see (6) and [5]. From the above, it can be concluded, that for the different types of high-viscosity liquids used here, the developed theory forecasts a stable bubble behavior.

Summary
The governing differential equation of a two-phase system consisting of a gas bubble immersed in a highly viscous, incompressible liquid medium is theoretically analyzed in detail with special focus to the mathematical significance of the non-Newtonian contributions of the second-order fluid on the bubble behavior over the whole timescale. Transforming the underlying problem into a well-structured nonlinear ODE makes possible the application of convergence arguments to show the existence of stable solutions under restrictions to the external pressure excitation and the named non-Newtonian coefficients. A first existence proof needs the convergence of the outer pressure excitation and is based on an equivalent ODE system of first order to enable the application of an appropriately defined Lyapunov function with required properties. These assumptions are also relevant for the main result which entails a strong correlated condition for the two non-Newtonian coefficients of the material law for a second-order liquid. Necessary to this end are boundedness conditions for the pressure comprising an upper bound for the time integral over the absolute value of the associated pressure difference, where the latter itself must be bounded and the pressure function has to be monotonically decreasing. Finally, by considering the linearized problem and using a variant of Gronwall's lemma, it is proved that the convergence speed of the solutions has exponential character. According to the linear theory of ordinary differential equations, this result may also be transferred to the original nonlinear system whose solutions achieve a stable rest state. Finally, the theoretical condition for the non-Newtonian contributions is compared to experimental values of different viscoelastic materials. It turns out that this strict coefficient condition is fulfilled for all considered non-Newtonian substances. Hence, after multiplication of the first integral term, and again by (50), the inequality (51) is obtained: