*2.3. Numerical Solution for Inner Nodes*

The convective terms are less important and are omitted from the main analyzed set of equations. This procedure will ensure that the use of interpolation, which greatly affects the numerical solution [2], is excluded. The simplified equations of continuity and momentum are identified as *L*<sup>1</sup> and *L*2, respectively:

$$L\_1 = \frac{\partial p}{\partial t} + c^2(\rho\_l - \rho\_v)\frac{\partial a}{\partial t} + 2\rho\_m c^2 \frac{\partial \varepsilon\_r}{\partial t} + \rho\_m c^2 \frac{\partial}{\partial x} \left(\frac{v}{a}\right) = 0 \tag{25}$$

$$L\_2 = \frac{\partial}{\partial t} \left( \frac{v}{\alpha} \right) + \frac{1}{\rho\_m} \frac{\partial p}{\partial x} + \frac{2}{\rho\_m R} \tau\_m = 0. \tag{26}$$

Combining linearly *L* = *ψL*<sup>1</sup> + *L*2, these equations by the unknown multiplier *ψ* gives:

$$
\psi \left[ \frac{\partial p}{\partial t} + \frac{1}{\rho\_m \psi} \frac{\partial p}{\partial x} \right] + \left[ \frac{\partial}{\partial t} \left( \frac{v}{a} \right) + \psi \rho\_m c^2 \frac{\partial}{\partial x} \left( \frac{v}{a} \right) \right] + \psi c^2 (\rho\_l - \rho\_v) \frac{\partial a}{\partial t} + 2 \psi \rho\_m c^2 \frac{\partial \varepsilon\_r}{\partial t} + \frac{2}{\rho\_m R} \mathfrak{r}\_{ll} = 0. \tag{27}
$$

By examination of the above Equation (27) with the definition of total derivatives, it can be noted that with:

$$\frac{d\chi}{dt} = \frac{1}{\rho\_m \psi} = \psi \rho\_m c^2 \tag{28}$$

it becomes the ordinary differential equation:

$$
\psi \frac{dp}{dt} + \frac{d}{dt} \left(\frac{v}{a}\right) + \psi c^2 (\rho\_l - \rho\_v) \frac{\partial a}{\partial t} + 2\psi \rho\_m c^2 \frac{\partial \varepsilon\_r}{\partial t} + \frac{2}{\rho\_m R} \tau\_m = 0. \tag{29}
$$

The solution of Equation (28) yields two particular values of *ψ*,

$$
\psi = \pm \frac{1}{c\rho\_m}.\tag{30}
$$

By inserting the above values into Equation (28) the particular relation between *x* and *t* is given as follows:

$$\frac{d\mathbf{x}}{dt} = \pm c.\tag{31}$$

This leads to a set of two equations on positive *C*<sup>+</sup> and negative *C*<sup>−</sup> characteristic lines:

$$\mathcal{C}^{+}: \; \frac{1}{c\rho\_{m}}\frac{dp}{dt} + \frac{d}{dt}\left(\frac{v}{a}\right) + \frac{c}{\rho\_{m}}(\rho\_{l} - \rho\_{v})\frac{\partial a}{\partial t} + 2c\frac{\partial \varepsilon\_{r}}{\partial t} + \frac{2}{\rho\_{m}R}\tau\_{m} = 0 \tag{32}$$

$$\mathcal{C}^-: -\frac{1}{c\rho\_m}\frac{dp}{dt} + \frac{d}{dt}\left(\frac{v}{a}\right) - \frac{c}{\rho\_m}(\rho\_l - \rho\_v)\frac{\partial a}{\partial t} - 2c\frac{\partial \varepsilon\_r}{\partial t} + \frac{2}{\rho\_m R}\tau\_m = 0. \tag{33}$$

From Equation (2), it follows that:

$$\alpha = \frac{\rho\_{\mathcal{W}} - \rho\_{\mathcal{v}}}{\rho\_l - \rho\_{\mathcal{v}}}.\tag{34}$$

Considering the above Equation (34), the third term on the left-hand side in Equations (32) and (33) can take the form:

$$\frac{c}{\rho\_m}(\rho\_l - \rho\_\upsilon)\frac{\partial a}{\partial t} = \frac{c}{\rho\_m}\frac{\partial}{\partial t}\left(\frac{\rho\_m - \rho\_\upsilon}{\rho\_l - \rho\_\upsilon}(\rho\_l - \rho\_\upsilon)\right) = \frac{c}{\rho\_m}\frac{\partial}{\partial t}(\rho\_m - \rho\_\upsilon) = \frac{c}{\rho\_m}\frac{\partial}{\partial t}\rho\_m. \tag{35}$$

Shu [8] writes the partial derivative of the density of the mixture in a logarithmic form:

$$\frac{c}{\rho\_m} \frac{\partial}{\partial t} \rho\_m = c \frac{\partial}{\partial t} \ln \left( \frac{\rho\_m}{\rho\_l} \right). \tag{36}$$

Note that under the logarithm, we have the division of *ρ<sup>m</sup>* by *ρl*, but the *ρ<sup>l</sup>* could be exchanged to any other constant, and the above equation would be satisfied anyway.

⎧ ⎪⎪⎨ ⎪⎪⎩ 1 *cρl* <sup>Δ</sup>(*p*−*κpv*)*A*−*<sup>D</sup>* <sup>Δ</sup>*<sup>t</sup>* <sup>+</sup> <sup>Δ</sup> Δ*t v α <sup>A</sup>*−*<sup>D</sup>* <sup>+</sup> *<sup>c</sup>* <sup>Δ</sup> <sup>Δ</sup>*<sup>t</sup> lnρ<sup>m</sup> ρl E* − *D F* − *A* + 2*c ∂ε<sup>r</sup> ∂t E*−*D* + <sup>2</sup> *<sup>ρ</sup>mR* (*τm*)*F*−*<sup>A</sup>* <sup>=</sup> <sup>0</sup> *f or dx dt* = *c* (37)

$$\begin{cases} \begin{array}{c} -\frac{1}{c\rho\_l} \frac{\Delta(p-\mathbf{x}\eta\_r)\_{B-D}}{\Delta l} + \frac{\Delta}{\Delta l} \left(\frac{\mathbf{z}}{a}\right)\_{B-D} - c \frac{\Delta}{\Delta l} \left(\frac{\rho\_m}{\rho\_l}\right)\_{E-D} - 2c \left(\frac{\partial c}{\partial t}\right)\_{E-D} + \frac{2}{\rho\_m \mathbf{z}^R} \left(\mathbf{r}\_{\mathcal{U}}\right)\_{G-B} = 0 \\ \mathbf{G} - \mathbf{B} \end{array} \tag{38}$$
 
$$for \frac{d\mathbf{x}}{dt} = -c$$

where *κ* = 1 + *c*2*ρl*ΞFΔt.

From both characteristics at the inner node, the following explicit system of equations is obtained:

$$\begin{cases} \begin{array}{c} \frac{1}{c\rho\_{l}} \left( (p\_{D} - \mathbf{x}p\_{\mathrm{r}}) - (p\_{A} - \mathbf{x}p\_{\mathrm{r}}) \right) + \frac{\left( \frac{\mathbf{r}\_{D}}{\mathbf{d}\_{D}} - \frac{\mathbf{r}\_{A}}{\mathbf{d}\_{A}} \right)}{\Delta t} + \frac{c}{2\Delta t} \left( \ln \frac{\rho\_{\mathrm{nD}}}{\rho\_{l}} - \ln \frac{\rho\_{\mathrm{nE}} \rho\_{\mathrm{nE}}}{\rho\_{l} \rho\_{\mathrm{nA}}} \right) + 2c \left( \frac{\mathbf{\mathcal{S}}\_{\mathrm{r}}}{\mathbf{d}t} \right)\_{\mathrm{E}-D} + \frac{2}{\rho\_{\mathrm{n}} \mathbb{R}} (\tau\_{\mathrm{nI}})\_{\mathrm{F}-A} = 0 \\\ -\frac{1}{c\rho} \frac{((p\_{D} - \mathbf{x}p\_{\mathrm{r}}) - (p\_{B} - \mathbf{x}p\_{\mathrm{r}}))}{\Delta t} + \frac{\left( \frac{\mathbf{\mathcal{S}}\_{\mathrm{D}}}{\mathbf{d}t} - \frac{\mathbf{r}\_{\mathrm{B}}}{\mathbf{d}t} \right)}{\Delta t} - \frac{c}{2\Delta t} \left( \ln \frac{\rho\_{\mathrm{nD}}}{\rho\_{l}} - \ln \frac{\rho\_{\mathrm{nE}} \rho\_{\mathrm{nE}}}{\rho\_{l} \rho\_{\mathrm{nB}}} \right) - 2c \left( \frac{\mathbf{\mathcal{S}}\_{\mathrm{r}}}{\mathbf{d}t} \right)\_{\mathrm{E}-D} + \frac{2}{\rho\_{\mathrm{n}} \mathbb{R}} (\tau\_{\mathrm{nI}})\_{\mathrm{G}-B} = 0 \end{array} . \end{cases} \tag{39}$$

The system can be further rewritten after introducing:

$$
\mathcal{L}c\left(\frac{\partial \varepsilon\_r}{\partial t}\right)\_{E-D} = p\_D \Sigma cF - \Sigma cG\_E(t) \tag{40}
$$

$$\frac{2}{\rho\_{n}R}(\tau\_{\mathrm{H}})\_{F-A} = \frac{f\_{A}\upsilon\_{A}|\upsilon\_{A}|}{4Ra\_{A}^{2}} + \frac{4\nu\_{mA}}{R^{2}} \sum\_{i=1}^{3} \underbrace{\left[A\_{i}\eta\_{iF} + \eta B\_{i}\left[\frac{\upsilon\_{A}}{a\_{A}} - \frac{\upsilon\_{F}}{a\_{F}}\right] + [1-\eta]\mathbb{C}\_{i}\left[\frac{\upsilon\_{F}}{a\_{F}} - \frac{\upsilon\_{F'}}{a\_{F'}}\right]\right]}\_{\mathrm{Y}\_{iA}} = \frac{\lambda\_{F-A}\upsilon\_{A}|\upsilon\_{A}|}{4Ra\_{A}^{2}}\tag{41}$$

$$\frac{2}{\rho\_{n}R}(\tau\_{\mathrm{n}})\_{G-B} = \frac{f\_{\mathrm{B}}v\_{\mathrm{B}}|v\_{\mathrm{B}}|}{4Ra\_{\mathrm{B}}^{2}} + \frac{4\nu\_{m\mathrm{B}}}{R^{2}} \underbrace{\sum\_{i=1}^{3} \left[A\_{i}y\_{i\mathrm{G}} + \eta B\_{i} \left[\frac{v\_{\mathrm{B}}}{a\_{\mathrm{B}}} - \frac{v\_{\mathrm{G}}}{a\_{\mathrm{G}}}\right] + [1-\eta]\mathbf{C}\_{i} \left[\frac{v\_{\mathrm{G}}}{a\_{\mathrm{G}}} - \frac{v\_{\mathrm{G}'}}{a\_{\mathrm{G}'}}\right] \right]}\_{y\_{\mathrm{B}}} = \frac{\lambda\_{G-B}v\_{\mathrm{B}}|v\_{\mathrm{B}}|}{4Ra\_{\mathrm{B}}^{2}}\tag{42}$$

where:

$$\begin{cases} \lambda\_{F-A} = f\_A + \frac{16\nu\_{nBA}}{R\nu\_A|v\_A|}\sum\_{i=1}^3 \underbrace{\left[A\_i y\_{iF} + \eta B\_i \left[\frac{v\_A}{a\_A} - \frac{v\_F}{a\_F}\right] + [1-\eta]\mathbb{C}\_i \left[\frac{v\_F}{a\_F} - \frac{v\_{F'}}{a\_{F'}}\right]\right]}\_{y\_{iA}} \\\\ \lambda\_{G-B} = f\_B + \frac{16\nu\_{nBA}\frac{a\_B}{R}}{R\nu\_B|v\_B|}\sum\_{i=1}^3 \underbrace{\left[A\_i y\_{iG} + \eta B\_i \left[\frac{v\_B}{a\_B} - \frac{v\_G}{a\_G}\right] + [1-\eta]\mathbb{C}\_i \left[\frac{v\_G}{a\_G} - \frac{v\_{G'}}{a\_{G'}}\right]\right]}\_{y\_{iB}} \end{cases} \tag{43}$$

By transforming the system of Equation (39) in a way that the parameters searched for a given inner node D of the characteristics grid (Figure 1) remain on the left-hand side, one obtains:

$$\begin{cases} \begin{aligned} \frac{p\_D}{c\rho\_l} - \frac{\kappa p\_v}{c\rho\_l} + \frac{v\_D}{a\_D} + \frac{c}{2}In\frac{\rho\_{\rm mD}}{\rho\_l} + p\_D c\Sigma F\Delta t = \mathbb{C}\_A \\\ -\frac{p\_D}{c\rho\_l} + \frac{\kappa p\_v}{c\rho\_l} + \frac{v\_D}{a\_D} - \frac{\xi}{2}In\frac{\rho\_{\rm mD}}{\rho\_l} - p\_D c\Sigma F\Delta t = \mathbb{C}\_B \end{aligned} \end{cases} \tag{44}$$

where *CA* as well *CB* are time-dependent functions that are iteratively calculated using the values known from the previous numerical time step:

$$\begin{cases} \mathbf{C}\_{A} = \frac{\upsilon\_{A}}{a\_{A}} + \frac{p\_{A} - \kappa p\_{v}}{c\rho\_{l}} - \frac{f\_{A}\Delta t \upsilon\_{A}|\upsilon\_{A}|}{4Ra\_{A}^{2}} + \frac{c}{2}ln\frac{\rho\_{m}\varepsilon\rho\_{m}\mathfrak{r}}{\rho\_{l}\rho\_{m\mathcal{A}}} + \mathcal{G}\_{E}(t)\Sigma c\Delta t \\\\ \mathbf{C}\_{B} = \frac{\upsilon\_{B}}{a\_{B}} - \frac{p\_{B} - \kappa p\_{v}}{c\rho\_{l}} - \frac{f\_{B}\Delta t \upsilon\_{B}|\upsilon\_{B}|}{4Ra\_{B}^{2}} - \frac{c}{2}ln\frac{\rho\_{m\mathcal{E}}\rho\_{m\mathcal{G}}}{\rho\_{l}\rho\_{m\mathcal{B}}} - \mathcal{G}\_{E}(t)\Sigma c\Delta t \end{cases} \tag{45}$$

**Figure 1.** Rectangular grid in the method of characteristics.

From the above equations, the final solutions for the inner node D (Figure 1) of the grid of characteristics are obtained. The solution for the mean velocity at cross-section is:

$$v\_D = \frac{\alpha\_D(\mathcal{C}\_A + \mathcal{C}\_B)}{2} \tag{46}$$

and the solution for the pressure is:

$$p\_D = \frac{(\mathbb{C}\_A - \mathbb{C}\_B)c\rho\_l}{2\kappa} + p\_v - \frac{c^2\rho\_l}{2\kappa} \ln \frac{\rho\_{mD}}{\rho\_l}.\tag{47}$$

The analysis of the above formulas shows that in order for *pD* > *pv* (note that then *ρmD* = *ρl*), the condition that *CA* ≥ *CB* must be met. When *CA* ≥ *CB*, there is no cavitation and *<sup>α</sup><sup>D</sup>* <sup>=</sup> 1; *pD* <sup>=</sup> (*CA*−*CB*)*cρ<sup>l</sup>* <sup>2</sup>*<sup>κ</sup>* + *pv*. Otherwise, when *CA* < *CB* cavitation occurs, then *pD* = *pv* and *ρmD* = *ρle*( *CA*−*CB <sup>c</sup>* ) .

Having the instantaneous value of the mixture density *ρmD*, the vapor density *ρv*, and the liquid density *ρl*, the instantaneous value of the liquid phase concentration *α<sup>D</sup>* should be determined from the formula (Equation (2):

$$\alpha\_{D} = \frac{\rho\_{l}e^{(\frac{\mathcal{C}\_{A}-\mathcal{C}\_{R}}{\varepsilon})} - \rho\_{\upsilon}}{\rho\_{l} - \rho\_{\upsilon}}.\tag{48}$$
