*2.2. Wall Shear Stress and Retarded Strain*

The wall shear stress in transient pipe flow can be calculated with the help of convolutional theory. Zielke [41] for laminar flow and later Vardy-Brown [42] for turbulent flow presented an initial version of this equation. For homogeneous bubble flow, the mixture density should be taken into account:

$$\tau\_{\mathfrak{m}} = \left(\frac{fv|v|\rho\_{\mathfrak{m}}}{8a^2} + \frac{2\mu\_m}{R} \int\_0^t \frac{\partial}{\partial t} \left(\frac{v}{a}(u)\right) \cdot w\_{\text{LF}}(t-u) du\right). \tag{17}$$

The function *wUF*(*t* − *u*) [-] is a so-called *weighting function*. In our work, this function is identical to the functions used in other cavitational and single-phase flow models. For detailed form and more information about weighting functions, please refer to paper [43].

The numerical modified effective solution of the above convolution integral, which is used in this work, is based on the improved solution for single-phase flows [44,45]:

$$\tau(t+\Delta t) \approx \frac{\rho\_{\rm mf}\upsilon\_{(t)}\Big|\upsilon\_{(t)}\Big|}{8a\_{(t)}^{2}} + \frac{2\mu\_{\rm mf}}{R} \sum\_{i=1}^{3} \underbrace{\left[A\_{i}y\_{i}(t) + \eta B\_{i}\left[\frac{\upsilon\_{(t)}}{a\_{(t)}} - \frac{\upsilon\_{(t-\Delta t)}}{a\_{(t-\Delta t)}}\right] + [1-\eta]C\_{i}\left[\frac{\upsilon\_{(t-\Delta t)}}{a\_{(t-\Delta t)}} - \frac{\upsilon\_{(t-2\Delta t)}}{a\_{(t-2\Delta t)}}\right]\right]}\_{y\_{i}(t+\Delta t)}.\tag{18}$$

where Δt [s] is a constant time step in the method of characteristics; *μ<sup>m</sup>* is Dukler's [46] two-phase mixture dynamic viscosity [Pa·s], and

$$\eta = \frac{\int\_0^{\Delta \hat{t}} w\_{class.}(u) du}{\int\_0^{\Delta \hat{t}} w\_{eff.}(u) du};\ A\_i = e^{-\eta\_i \Delta \hat{t}};\ B\_i = \frac{m\_i}{\Delta \hat{t} n\_i} [1 - A\_i];\ \mathcal{C}\_i = A\_i B\_i. \tag{19}$$

The *mi* and *ni* coefficients are determined from the analytical formulas presented in a recent paper [47], and Δˆ*t* = *<sup>μ</sup><sup>m</sup> <sup>ρ</sup>mR*<sup>2</sup> <sup>Δ</sup>*<sup>t</sup>* is dimensionless time. In the case of turbulent flow, these coefficients must be scaled in accordance with the guidelines presented in paper [48].

The next step is to evaluate the partial derivative of retarded strain using a convolution integral. According to [49,50], the retarded strain can be written in a simpler form than the original one [38]:

$$\frac{\partial \varepsilon\_r(t)}{\partial t} = \frac{D}{2\varepsilon} \tilde{\zeta} \int\_0^t \frac{\partial}{\partial t} (p(u) - p(0)) \cdot \left( \sum\_{i=1}^k \frac{I\_i}{T\_i} e^{-\frac{t-u}{T\_i}} \right) du = \frac{\Xi}{2} \int\_0^t \frac{\partial p(u)}{\partial t} \cdot w\_I(t-u) du \tag{20}$$

where *wJ*(*<sup>t</sup>* − *<sup>u</sup>*) is the creep weighting function [Pa−1·s−1]; *<sup>J</sup>*<sup>i</sup> is the creep compliance of the i spring of the Kelvin–Voigt element [Pa<sup>−</sup>1]; *T*<sup>i</sup> is the retardation time of the dashpot of i-element [s].

Worth noting is the analogy of the above convolution integral with the convolutional integral representing the wall shear stress ∑*<sup>k</sup> i*=1 *Ji Ti e* <sup>−</sup> *<sup>t</sup>*−*<sup>u</sup> Ti* = *wJ*(*t* − *u*). This analogy made it possible to solve numerically the above convolution integral in an effective manner using Schohl's effective scheme [51]:

$$\frac{\partial \varepsilon\_{r}}{\partial t}(t + \Delta t) \approx \frac{\Xi}{2} \sum\_{i=1}^{k} \underbrace{\left(z\_{i}(t) \cdot e^{-\frac{\Delta t}{I\_{i}}} + \frac{I\_{i}}{\Delta t} \left[1 - e^{-\frac{\Delta t}{I\_{i}}}\right] \left(p\_{(t + \Delta t)} - p\_{(t)}\right)\right)}\_{z\_{i}(t + \Delta t)}.\tag{21}$$

The above equation (Equation (21) may be written in a simpler form:

$$\frac{\partial \varepsilon\_r}{\partial t}(t + \Delta t) = \left(p\_{(t + \Delta t)}F - G(t)\right)\frac{\Sigma}{2} \tag{22}$$

where:

$$F = \sum\_{i=1}^{k} M\_{\bar{i}} ; \; G(t) = \sum\_{i=1}^{k} \left( M\_{\bar{i}} p\_{(t)} - z\_{\bar{i}}(t) \cdot N\_{\bar{i}} \right) ; \; M\_{\bar{i}} = \frac{f\_{\bar{i}}}{\Delta t} \left[ 1 - e^{-\frac{\Delta t}{T\_{\bar{i}}}} \right] ; \; N\_{\bar{i}} = e^{-\frac{\Delta t}{T\_{\bar{i}}}} . \tag{23}$$

The use of convolution integrals (Equations (17) and (20)) in the basic system of Equation (16) results in:

$$\begin{cases} \begin{array}{c} \rho\_{m}\frac{d}{dt}\left(\frac{\upsilon}{a}\right) + \frac{\partial p}{\partial x} + \frac{2}{\overline{\kappa}} \left(\frac{fv|v|\rho\_{m}}{8a^{2}} + \frac{2\mu\_{m}}{\overline{\kappa}} \int\_{0}^{t} \frac{\partial}{\partial t} \left(\frac{\upsilon}{a}(u)\right) \cdot w\_{\text{LF}}(t-u) du\right) = 0\\ \frac{1}{\rho\_{m}c^{2}} \frac{dp}{dt} + \frac{(\rho\_{l} - \rho\_{v})}{\rho\_{m}} \frac{d\kappa}{dt} + \frac{\partial}{\partial x} \left(\frac{\upsilon}{a}\right) + \Sigma \int\_{0}^{t} \frac{\partial p(u)}{\partial t} \cdot w\_{I}(t-u) du = 0 \end{array} . \end{cases} . \tag{24}$$
