<figure><img src="image_1.png" /></figure>

HAL
open science

# A 3D numerical strategy for the computations of shock-induced bubble collapse near a wall

Ksenia Kozhanova, Yannick Hoarau, Eric Goncalves da Silva

<figure><img src="image_2.png" /></figure>

## To cite this version:

Ksenia Kozhanova, Yannick Hoarau, Eric Goncalves da Silva. A 3D numerical strategy for the computations of shock-induced bubble collapse near a wall. Computers and Fluids, 2025, 293, pp.106609. 10.1016/j.compfluid.2025.106609. hal-05017998

HAL Id: hal-05017998

https://hal.science/hal-05017998v1

Submitted on 2 Apr 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A 3D numerical strategy for the computations of
shock-induced bubble collapse near a wall.

Ksenia Kozhanova$^{a,b,*}$, Yannick Hoarau$^{c}$, Eric Goncalves da Silva$^{b}$

$^{a}$Aix Marseille University, CNRS, Centrale Marseille, M2P2, 13451 Marseille, France

$^{b}$ISAE-ENSMA, Institut Pprime, UPR 3346 CNRS, 86961 Chasseneuil Cedex, France

$^{c}$University of Strasbourg, ICUBE, UMR 7357 CNRS, 67000 Strasbourg, France

**Abstract**

The importance of modelling two-phase flows involving shock waves arises from
many engineering and medical applications. The presence of strong shock waves
and their interactions with bubble interfaces, the high density ratio between
phases and the large variation of material properties makes the resolution of
such problems a complicated task for the numerical methods. While the variety
of numerical techniques to solve these problems exist, e.g. the sharp interface or
the diffuse interface methods, these strategies can lead to spurious oscillations of
the solution near the interface. It is well known that it is difficult to achieve both
a high order accuracy of the scheme and the monotonicity of the solution. In this
paper a four-equation two-phase model is employed and integrated in an explicit
fully parallelised finite-volume solver with HLLC numerical scheme coupled with
WENO reconstruction methods and Hancock predictor-corrector scheme and
non-uniform mesh based on stretching function in order to compute a 3D shock-
induced bubble collapse near a wall. The novelty of our work is improved
accuracy of computations of such a problem with optimised computational cost
thanks to the non-uniform mesh introduction in 3D computations.

Keywords: Diffuse interface methods, Riemann solver, WENO scheme,
bubble collapse, shock waves

<footnote>*Corresponding author.
Email address: ksenia.kozhanova@univ-amu.fr (Ksenia Kozhanova)</footnote>
# 1. Introduction

Cavitation erosion is an important problem which has to be tackled in many domains, such as hydraulic and marine applications. This process occurs as a result of bubble collapse near the solid walls causing the material damage. The physical process of such phenomena has been studied in many experimental and numerical investigations by examining the single bubble collapse under shock wave loading (for instance, [1, 2, 3, 4]) in order to understand in detail the bubble collapse mechanism and consequential effect on the material. It has been understood that generated intense pressure waves and their interactions between each other and high-speed jets generated as a result of bubble collapse are main reasons for material damages and even micro-deformations.

From numerical computations point of view, the modelling of bubble collapse is challenging since it involves a simulation of a stiff multi-phase flows kind. Even though the industrial demand of such a modelling is high, its complexity is yet a challenge for the applied mathematics and numerical methods due to high variations of thermodynamic properties, strong pressure waves and their interaction with interfaces and the involved fast dynamics.

There are several widely used strategies for such two-phase problems. One of the often used is a method based on an average approach. This method led to the development of the variety of different perspectives based on the physical assumptions for the slip condition between phases and mechanical and thermo-dynamic equilibrium. The most complete modelling approach is known to be the two-fluid model [5], which has a high computational complexity since it consists of seven equations. An alternative method, which is more suitable for practical applications, is a reduced five-equation model with the underlying assumption of pressure and velocity equilibrium between phases [6, 7]. Furthermore, this model can be reduced by using the assumption of thermal equilibrium. Such a model consists of three equations for conservation laws and a fourth one for the void ratio [8].

One of the important and difficult domains of studying the two-phase com-
pressible flows is the suitable numerical methods. One challenge comes from the characteristic waves structure of the hyperbolic system describing this physical phenomena, which yields to the unstable simulations. Moreover, the problem of discontinuities of thermodynamic variables and equations of state occurring due to the material interfaces causes difficulties in the derivation of Riemann solvers. The speed of sound component is also a challenge. Due to the large differences between its values in mixture and liquid phases, some large variations are present and cause the non smooth behaviour and inaccuracies in the solution. Thus, non-physical oscillations can appear in the solution of the whole wave structure [9]. The numerical dissipation of the numerical schemes which is required to follow the thermodynamic law is also a result of these inaccuracies of the modelling.

In order to improve the order of accuracy of finite-volume solvers, high-order reconstruction methods have to be used. However, the choice of such methods is not straightforward. The numerical instabilities due to the material interface, for instance, can be improved by using the primitive variables reconstruction, as shown in [10], where authors successfully applied Weighted Essentially Non-Oscillatory (WENO) for gas-gas two-phase problems. However, this strategy is not as suitable when a large density ratio between phases is present leading to high gradients in the flow. Indeed, it has been shown that high-order schemes can lead to spurious oscillations when applied to stiff compressible two-phase flows. Two methods have been suggested to tackle this issue in [11]. One of them, the hybrid construction of the scheme, where the fifth order WENO is locally reduced to the second order reconstruction near the interface. Another method is a smoothing the material interface as a part of initial condition. However, there are disadvantages in both strategies. The first method is problem-dependent, since the interface condition is different based on the physical phenomena involved. On the other hand, the second method has a problem of thickening the interface artificially, which can cause further numerical instabilities. In order to avoid these problems, an alternative method have been suggested by Wang et al. [12], where an incremental 2- and 3-point
stencil is used to obtain a fifth order WENO scheme. This choice between the smaller and larger stencil is conditioned based on the discontinuity criteria.

We consider a gas bubble immersed in water and collapsing under the inci-
dence of shock wave in the vicinity of a wall [13]. The small spatial-temporal
scales associated to physical phenomena lead to the requirement of very high
number of points in computational domain and very small $\Delta t$. The CPU cost
becomes quickly prohibitive. Thus, in addition to the high-order numerical
methods we require an introduction of non-uniform mesh. Several mesh reduc-
tion methods exist. The adaptive mesh refinement (AMR), one of the most often
used in recent literature, is a method where the solution accuracy is adapted
based on the location of certain critical solution components. The main un-
derlying idea of this method is to apply higher numerical precision by means
of finer mesh in the zones of computational domain where the high gradients
are detected and reduction of this precision otherwise. The original AMR algo-
rithm had been presented in [14] and [15] where authors presented a dynamic
gridding or local adaptive mesh refinement. An alternative and simpler non-
uniform mesh method is a mesh stretching. The mesh stretching can be achieved
by introduction of some function to distribute the points non-uniformly. The
easiest function can be based on the geometrical progression and more com-
plex one is on smoother functions, such as sinus, tangent etc. The analytical
studies of these functions suitable for the grid generation have been studied in
[16, 17, 18, 19]. While both mesh reduction techniques have been successfully
used to introduce non-uniform mesh, some difficulties may arise. Namely, the
suitability of the numerical schemes for such meshes can be questioned, since
majority of the schemes are derived for the uniform meshes. The modified
derivation of the high-order reconstruction is required should it be used with
not equal grid spacing. The WENO-type methods have been addressed in the
framework of non-uniform grid in [20]. Present study discusses the class of
stretching functions based on a one-sided stretching function by using hyper-
bolic tangent distribution by [21]. The strategy of appropriate mesh stretching
for the problem of shock-induced bubble collapse near a wall is proposed.
The present paper aims to perform a 3D computations of the shock bubble
collapse by using the four-equation model presented in [22] and by means of 3D
finite volume based solver with high-order spacial discretisation. Successful 3D
computations of similar problem have been already performed in [8] and [23]
on uniform mesh with second order accuracy by using MUSCL discretisation.
The novelty of present work is to extend these studies in terms of well-resolved
3D shock-induced bubble collapse by using high-order discretisation and non-
uniform mesh for specific test case from [23], i.e. 3D shock-induced bubble
collapse next to a wall with a stand-off distance equals to two. This study could
be used as a benchmark test case for potential future numerical validation and a
starting point for the future pressure loads on solid material studies with higher
accuracy.

The paper is organised as follows: The governing equations and equations
of state are provided in Section 2, followed by their discretisation in Section 3.
Section 4 provides an overview of the optimised parallel 3D solver SCB with
numerical results of a highly-resolved shock-induced bubble collapse near a wall
by using third order discretisation and stretched mesh. We conclude our work
in Section 5.

## 2. Equations and models

The present study is based on the one-fluid mixture approach with following
assumptions:

* strong coupling of phases with the same velocity

• thermal and mechanical equilibrium between the phases

• saturation state of liquid

• viscous effects and surface tension are neglected

2.1. *Four-equation model*

The model is based on a four-equation system where three equations consist
of conservation laws for mixture quantities and the fourth one of the void ratio
transport equation [4]. Using the conservative variables representation, i.e. $w = (\rho, \rho\vec{V}, \rho E)$ and the volume fraction of gas $\alpha$, the 2D inviscid system can be written as

$$ 
\begin{aligned}
& \frac{\partial \rho}{\partial t} + \text{div}(\rho \vec{V}) = 0 \\
& \frac{\partial (\rho \vec{V})}{\partial t} + \text{div}(\rho \vec{V} \otimes \vec{V} + P I_d) = \vec{0} \\
& \frac{\partial (\rho E)}{\partial t} + \text{div}(\rho \vec{V} H) = 0 \\
& \frac{\partial \alpha}{\partial t} + \vec{V} \cdot \text{grad}(\alpha) = K \text{div}(\vec{V})
\end{aligned}
\qquad (1) $$

where $\vec{V} = (u, v)$ denotes the velocity, $E = e + V^2/2$ is the total energy of the mixture and $H = h + V^2/2$ its total enthalpy, while $e$ and $h$ are internal energy and enthalpy, respectively. The reflection of the change in each phase volume and speed of sound of pure phases $c_k$ are included into the term $K$, which is derived according to the following formulation,

$$ K = \frac{\rho_l c_l^2 - \rho_g c_g^2}{\frac{\rho_l c_l^2}{1-\alpha} + \frac{\rho_g c_g^2}{\alpha}} $$

where $c_g$ and $c_l$ are constants and index $g$ stands for the gas phase and $l$ for the liquid one.

An equation of state (EOS) for the mixture is required to close the system. The formulation is deduced from convex stiffened gas EOS [24] for both phases and obtained from the thermal and mechanical equilibrium assumptions. The equations of mixture for the pressure and temperature are,

$$ P(\rho, e, \alpha, Y) = (\gamma(\alpha) - 1)\rho(e - q(Y)) - \gamma(\alpha)P_{\infty}(\alpha), \qquad (2) $$

$$ T(\rho, h, Y) = \frac{h - q(Y)}{C_p(Y)}, \qquad (3) $$

where $\gamma(\alpha)$ and $P_{\infty}(\alpha)$ are mixture parameters deduced from the mass frac-
tion of gas $Y = \alpha \frac{\rho_g}{\rho}$ and

$$ \frac{1}{\gamma(\alpha) - 1} = \frac{\alpha}{\gamma_g - 1} + \frac{1 - \alpha}{\gamma_l - 1} \\ q(Y) = Yq_g + (1 - Y)q_l \qquad (4) $$

$$ C_p(Y) = YC_{pg} + (1 - Y)C_{pl} $$

$$ P_{\infty}(\alpha) = \frac{\gamma(\alpha) - 1}{\gamma(\alpha)} \left[ \alpha \frac{\gamma_g}{\gamma_g - 1} P_{\infty}^{g} + (1 - \alpha) \frac{\gamma_l}{\gamma_l - 1} P_{\infty}^{l} \right] \qquad (5) $$

where $C_p$ and $C_v$ are thermal capacities, $q$ is the energy of formation, $\gamma = C_p/C_v$ is the heat capacity ratio and $Y = \alpha \rho_v / \rho$ is the mass fraction of gas. The quantity $P_{\infty}$ denotes a constant reference pressure.

The present system is hyperbolic and the eigenvalues are $\lambda_1 = u - c_{wallis}$, $\lambda_{2,3,4} = u$, $\lambda_5 = u + c_{wallis}$, where $c_{wallis}$ corresponds to the propagation of the acoustic waves without mass and heat transfer and has the following value,

$$ \frac{1}{\rho c_{wallis}^2} = \frac{\alpha}{\rho_g c_g^2} + \frac{1 - \alpha}{\rho_l c_l^2} \qquad (6) $$

where $c_g$, $c_l$ is the speed of sound of vapour and liquid phases, respectively.

The presented system could be derived from the two-component model (involving phasic mass equations) assuming the mechanical equilibrium between phases (see [25]). However, mathematically, the equivalence between these systems have not been demonstrated. The similarity of the void transport equation in (1) could be noted in comparison with the five-equation model (so-called Kapila model, [26]). However, the latter is a two-temperature formulation in which a heat transfer term appears in the void ratio equation (see [7], for example). The present four-equation model assumes the temperature equilibrium between phases. This four-equation model showed its robustness and its ability to compute accurately various cases such as cavitating flows [22, 27], shock/expansion tube problems [4] and shock-bubble interactions [8, 23]. Moreover, compared to the five-equation model, similar results were obtained on shock-bubble interactions with a lower computational cost [8].

One of the ongoing debates about such a system is, however, its non-conservative formulation. Indeed, in the presence of shock waves the convergence of the nu-
merical scheme to the exact solution is not guaranteed. Yet, practically, this model has been validated on various shock tube problems against existing reference results or analytical solutions. Reader is referred to the discussion in [28], for example.

### 3. Numerical discretisation

The system (1) is discretised by using the finite volume method with explicit time integration. This section discusses the basic discretisation in space and time followed by high-order spacial numerical method, WENO3-Z. The later is used in order to obtain the accurate converged results in shorter computational time.

The choice of both, basic and advanced discretisation, is based on the existing literature focused on the compressible two-phase flows and the ability of the methods to capture the shocks robustly while resolving smooth solutions accurately. Furthermore, the derivation of the schemes to account for the non-uniform mesh should be readily accessible.

#### 3.1. Basic discretisation

We start by rewriting the system (1) in 1D form which is used for the directional splitting. The resulting formulation is,

$$ \frac{\partial U}{\partial t} + \frac{\partial G(U)}{\partial x} + B(U) \frac{\partial u}{\partial x} = 0, \qquad (7) $$

where $U = \begin{pmatrix} \mathbf{w} \\ \alpha \end{pmatrix}$, $G(U) = \begin{pmatrix} F(\mathbf{w}) \\ \alpha u \end{pmatrix}$, $B(U) = \begin{pmatrix} 0 \\ -K - \alpha \end{pmatrix}$, and $F$ denotes the convective flux.

We note that present non-conservative formulation poses several known problems, e.g. the difficulties associated with sharp gradients reconstruction without non-physical oscillations, accurate reconstruction of smooth areas of the solution etc.

Presently, the computational spatial and temporal domains are divided into regular meshes of the uniform length $\Delta x$ and uniform intervals $\Delta t$, respectively.
We denote a cell by index $i$ and its interfaces on the left and right by $i-1/2$ and $i+1/2$, respectively. The finite volume method is used to obtain the discretised form of (1),

$$ \Delta x \frac{U_i^{n+1} - U_i^n}{\Delta t} + G_{i+1/2}^n - G_{i-1/2}^n + \int_{x_{i-1/2}}^{x_{i+1/2}} B(U^n) \frac{\partial u}{\partial x} dx = 0. \quad (8) $$

The discretisation of the non-conservative term is discussed below. While many formulations are available for the estimation of the numerical flux $G$, such as Roe-type [29], AUSM-type [30], Jameson-Schmidt-Turkel [31], the HLLC [32] scheme is used in this work due to its known good performance in presence of shocks and discontinuities.

A brief presentation of above mentioned HLLC scheme, the approximate Riemann solver, is proposed here. It is based on two intermediate states, $w_L^*$ and $w_R^*$, which are separated by the speed contact wave $S_M$. The conservative part of the equation (7) is discretised according to the classical HLLC flux formulation [32]. The discretisation of the non-conservative part, on the other hand, is based on the approach proposed in [33]. The $B(U)$ term is approximated as

$$ \int_{x_{i-1/2}}^{x_{i+1/2}} B(U) \frac{\partial u}{\partial x} dx \approx B(U_i)(u_{i+1/2} - u_{i-1/2}). \quad (9) $$

The interface cell $u_{i+1/2}$ is then computed as,

$$ u_{i+1/2} = \begin{cases} u_L & \text{if } S_L > 0, \\ \frac{S_L - u_L}{S_L - S_M} S_M & \text{if } S_L \le 0 < S_M, \\ \frac{S_R - u_R}{S_R - S_M} S_M & \text{if } S_M \le 0 \le S_R, \\ u_R & \text{if } S_R < 0. \end{cases} \quad (10) $$

The HLLC formulation standalone leads to the first order scheme in space. However, the spatial order can be straight forwardly increased by using MUSCL extrapolation technique with TVD slopes limiters to avoid non-physical oscillations or more sophisticated high order numerical methods. While the focus in this study is on the latter, the results obtained by using MUSCL formula-
tions coupled with the van Albada slope limiter are presented for comparison purposes.

In terms of the time integration, the solver uses a predictor-corrector scheme, referred as the Hancock's method [34], which is discussed in detail in the next subsection. The numerical treatment of the boundary conditions is based on the use of the characteristic relations of the four-equation system (see [8] for more details).

### 3.2. Hancock's predictor-corrector

The main interest of this method lies in its remarkable robustness and the ability to carry out stable computations of complex problems involving discontinuities and shocks [35]. The method can be sketched by using the following steps:

1. denoting the set of variables $\mathbf{W}$ and a reconstructed variable $W$, perform the reconstruction on linear subcell distribution over the whole chosen set. In general, the choice of variables does not alter structure of algorithm. We define

$$ \delta W_j = \text{limiter} (W_j - W_{j-1}, W_{j+1} - W_j), \quad (11) $$

2. perform the advancement of solution by half time step, using the following form of equations

$$ \mathbf{W}_t + (\mathbf{A}_W)_j \mathbf{W}_x = 0, \quad (12) $$

where $\mathbf{A}_W$ is a matrix which depends on the variables choice for the reconstruction. For the primitive variables, this matrix $\mathbf{A}_W$ in direction $x$ is
$$A_x = \begin{pmatrix} u & \rho & 0 & 0 & 0 \\ 0 & u & 0 & \frac{1}{\rho} & 0 \\ 0 & 0 & u & 0 & 0 \\ 0 & \rho c^2 & 0 & u & 0 \\ 0 & -K & 0 & 0 & u \end{pmatrix}.$$

Note, that the source term is included to this matrix A. The predictor step is then computed as,

$$\tilde{\mathbf{W}}_j = \mathbf{W}_j - \frac{\Delta t}{2\Delta x} \mathbf{A}_W \delta \mathbf{W}_j \quad (13)$$

3. using previously computed gradients, solve for the time-centered interface values,

$$\begin{aligned} \tilde{\mathbf{W}}_{j-\frac{1}{2}R} &= \tilde{\mathbf{W}}_j - \frac{1}{2} \delta \mathbf{W}_j, \\ \tilde{\mathbf{W}}_{j+\frac{1}{2}R} &= \tilde{\mathbf{W}}_j + \frac{1}{2} \delta \mathbf{W}_j, \end{aligned} \quad (14)$$

4. perform the computation for time-centered interface fluxes using the chosen numerical approximation flux scheme (e.g. HLLC), i.e.

$$\tilde{\mathbf{F}}_{j+\frac{1}{2}} = \mathbf{F}(\tilde{\mathbf{W}}_{j+\frac{1}{2}L}, \tilde{\mathbf{W}}_{j+\frac{1}{2}R}), \quad (15)$$

5. advance the solution over the full time-step

$$\mathbf{U}_j^{n+1} = \tilde{\mathbf{U}}_j - \frac{\Delta t}{\Delta x}(\tilde{\mathbf{F}}_{j+\frac{1}{2}} - \tilde{\mathbf{F}}_{j-\frac{1}{2}}) - \frac{\Delta t}{\Delta x} B(U_j)(\tilde{u}_{j+1/2} - \tilde{u}_{j-1/2}). \quad (16)$$

The reconstruction step can be applied to either the primitive or characteristic variables. Some authors note the latter choice leads to more robust solutions. Should the characteristic variables be chosen, the limiting step (11) is to be replaced by the following transformation step

$$\delta W_j = R(\mathbf{W}_j)\text{limiter} (\Delta_j - \Delta_{j-1}, \Delta_{j+1} - \Delta_j) \quad (17)$$

where $\Delta_i^j = L(W_{j+1} - W_j)$, L is the matrix of left eigenvectors arranged as rows, and R is the matrix of right eigenvectors arranged as columns. Since both
Jacobian matrix and primitive coefficient matrix relates to the same transformation, the left and right eigenvectors associated with primitive variables can be straightforwardly derived. For instance, the following associated matrices have been derived for the direction $x$,

$$R_x = \begin{bmatrix} 1 & 0 & 0 & -\frac{\rho}{c} & \frac{\rho}{c} \\ 0 & 0 & 0 & 1 & 1 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\rho c & \rho c \\ 0 & 0 & 1 & \frac{K}{c} & -\frac{K}{c} \end{bmatrix} \quad L_x = \begin{bmatrix} 1 & 0 & 0 & 0 & \frac{\rho}{K} \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & \frac{\rho c^2}{K} \\ 0 & 1 & 0 & -\frac{1}{\rho c} & 0 \\ 0 & 1 & 0 & \frac{1}{\rho c} & 0 \end{bmatrix} \quad (18)$$

The limiting formulation (11) can be performed by either simple MUSCL TVD techniques [36] or high-order numerical reconstructions methods.

### 3.3. High-order numerical reconstruction

The strategy to improve the spatial order is to use a higher order reconstruction method, i.e. by changing the calculation of the slope limiter $\delta W_j$. The present work concerns an improved third and fifth order WENO formulations, i.e. WENO3-Z and WENO5-IS. The detailed formulations can be found in [37] and [12] for WENO3-Z and WENO5-IS, respectively. The derivations of classical WENO scheme has been presented in [38]. The modified formulations to account for the non-uniform meshes are given below and have been adapted from [20].

**Spatial stencil** In order to reformulate the key formulations of WENO class of schemes, the formal formulation for the quadratic polynomial is recalled from [38],

$$P_r(x) = \sum_{j=0}^{2} C_{rj}(x) \bar{u}_{i-r+j}, \text{ with } r = 0, 1, 2, \quad (19)$$

where $\bar{u}$ denotes the cell average, superscript $i$ is dropped and

$$C_{rj}(x) = B_{rj}(x)h_{3-r+j}, \text{ where } h_m = \Delta x_{i-3+m}, m = 1, \dots, 5 \quad (20)$$
$$B_{rj}(x) = \sum_{m=j+1}^{3} \frac{\sum_{l=0, l \neq m}^{3} (\prod_{q=0, q \neq m, l}^{3} (x - x_{i-r+q-\frac{1}{2}}))}{\prod_{l=0, l \neq m}^{3} (x_{i-r+m-\frac{1}{2}} - x_{i-r+l-\frac{1}{2}})} \quad (21)$$

where $r, j = 0, 1, 2$.

Now, by defining $\hat{c}_{rj} = C_{rj}(x_{j+\frac{1}{2}})$ and $\hat{b}_{rj} = B_{rj}(x_{j+\frac{1}{2}})$ and recalling that WENO reconstruction is

$$w_{i+1/2}^{+} = \sum_{r=0}^{2} \omega_{r}^{i+1} P_{r}^{i+1}(x_{i+1/2}) \quad \text{and} \quad w_{i+1/2}^{-} = \sum_{r=0}^{2} \omega_{r}^{i} P_{r}^{i}(x_{i+1/2}) \quad (22)$$

we can compute $P_r(x_{i+\frac{1}{2}})$,

$$P_r(x_{i+\frac{1}{2}}) = \sum_{j=0}^{2} \hat{c}_{rj} \bar{u}_{i-r+j}, \quad \text{where} \quad \hat{c}_{rj} = \hat{b}_{rj} h_{3-r+j} \quad (23)$$

and

$$\hat{b}_{rj} = \sum_{m=j+1}^{3} \frac{\sum_{l=0, l \neq m}^{3} (\prod_{q=0, q \neq m, l}^{3} (x_{i+\frac{1}{2}} - x_{i-r+q-\frac{1}{2}}))}{\prod_{l=0, l \neq m}^{3} (x_{i-r+m-\frac{1}{2}} - x_{i-r+l-\frac{1}{2}})} \quad (24)$$

We provide the explicit formulations for the spacing $h_m = \Delta x_{i-3+m}$ with $m = 1, \dots, 5$, for $\hat{b}_{rj}$ in Appendix A.

Similarly, we can define $P_r(x_{i-\frac{1}{2}}) = \sum_{j=0}^{2} c_{rj} \bar{u}_{i-r+j}$ with $\hat{c}_{rj} = C_{rj}(x_{j-\frac{1}{2}}) = \hat{b}_{rj} h_{3-r+j}$ and

$$\hat{b}_{rj} = B_{rj}(x - j - \frac{1}{2}) = \sum_{m=j+1}^{3} \frac{\sum_{l=0, l \neq m}^{3} (\prod_{q=0, q \neq m, l}^{3} (x_{i-\frac{1}{2}} - x_{i-r+q-\frac{1}{2}}))}{\prod_{l=0, l \neq m}^{3} (x_{i-r+m-\frac{1}{2}} - x_{i-r+l-\frac{1}{2}})} \quad (25)$$

We can demonstrate that $b_{rj} = \hat{b}_{r-1,j}$ while $j = 0, 1, 2$ and $r = 1, 2$. Then, the only coefficients left to be computed are when $r = 0$ and $j = 0, 1, 2$. See Appendix A for coefficients, $b_{0,j}$.

Thus, we have explicit formulations for the approximate solution at the cell boundaries $x_{i-1/2}$ and $x_{i+1/2}$, $P_r(x_{i-\frac{1}{2}})$ and $P_r(x_{i+\frac{1}{2}})$, respectively with spacial
stencils formulations needed to account for the non-uniform mesh. The coefficients $c_{rj}$ and $\hat{c}_{rj}$ have to be computed only once since they depend only on the spatial spacing $h_m$.

**Smoothness measure** We recall that the smoothness indicator function are introduced in WENO class of the methods in order to improve the solution on the border between the areas with discontinuities and area where the solution is smooth. This function is formulated in classical WENO as,

$$ \beta_r = \int_{x-1/2}^{x+1/2} h_3(P'_r(x))^2 dx + \int_{i-1/2}^{i+1/2} (h_3)^3 (P''_r(x))^2 dx \quad (26) $$

Note, that when the above expression satisfies the condition that the approximate solution $P_r$ is smooth in $I_r$, that is, the computational stencil $S_r$ is smooth, we have,

$$ \beta_r = (P'_r(x_i)h_3)^2(1 + \mathcal{O}((h_3)^2)) \quad (27) $$

If, however, the stencil $S_r$ is not smooth, the following is satisfied,

$$ \beta_r = \mathcal{O}(1) \quad (28) $$

By using the fact that $P_r(x)$ is a quadratic polynomial, it follows that $P'_r(x)$ and $P''_r$ are linear and constant, respectively. The latter one is derived as,

$$ P''_r = \sum_{j=0}^{2} C''_{rj} \bar{u}_{i-r+j} = \sum_{j=0}^{2} B''_{rj} h_{3-r+j} \bar{u}_{i-r+j} \quad (29) $$

The argument $x$ has been dropped since $C''_{rj}(x)$ and $B''_{rj}$ are constants. We now can have the expressions for the $B''_{rj}$ (see Appendix A).

By substituting these coefficients into (26), the second integral of the equation can be rewritten,

$$ \int_{x_{i-1/2}}^{x_{i+1/2}} h_3(P''_r(x))^2 dx = (h_3)^4 \left( \sum_{j=0}^{2} B''_{rj}(x) h_{3-r+j} \bar{u}_{i-r+j} \right)^2 \quad (30) $$
The first integral can be reformulated by using the property that if the first derivative $P'_r(x)$ is linear in $x$, then $(P''_r(x))^2$ is quadratic in $x$. That leads to the way to compute the exact value for the first integral in (26) by applying the Simpson's quadrature rule,

$$ \int_{x_{i-1/2}}^{x_{i+1/2}} h_3(P'_r(x))^2 dx = (h_3)^2((P'_r(x-1/2))^2+4(P'_r(x_i))^2+(P'_r(x_{i+1/2}))^2) \quad (31) $$

It follows from 19, that the values for $P'_r(x_{i-1/2})$, $P'_r(x_i)$, $P'_r(x_{i+1/2})$ are,

$$ P'_{r}(x) = \frac{1}{6} \sum_{j=0}^{2} C'_{rj}(x) \bar{u}_{i-r+j} = \frac{1}{6} \sum_{j=0}^{2} B'_{rj}(x) h_{3-r+j} \bar{u}_{i-r+j} \quad (32) $$

The formulations for the coefficients $B'_{rj}$ are obtained from (21), (see Appendix A.

**Weights coefficients.** We recall that the fifth order approximation to the solution $u(x, t_n)$ is unique and, thus the weights can be derived from following relations,

$$ \tilde{P}(x_{i-\frac{1}{2}}) = \sum_{r=0}^{2} d_r P_r(x_{i-\frac{1}{2}}) \quad \text{and} \quad \tilde{P}(x_{i+\frac{1}{2}}) = \sum_{r=0}^{2} \hat{d}_r P_r(x_{i+\frac{1}{2}}) \quad (33) $$

The explicit formulations for the constants $d_r$ and $\hat{d}_r$ can be found in Appendix Appendix A.

Lastly, the stability condition for the above described discretization follows standard CFL calculations, i.e. $CFL = \Delta t(\sum_{i=1}^{n} \frac{u_i}{\Delta x_i}) \le C_{max}$, where $C_{max} = 1$. Generally, CFL number should not be above 0.3 or even 0.1 for the test cases presented in this paper. Both, $\Delta t$ and/or CFL number are specified for each test case presented below.

# 4. Validation and preliminary results

The solver and the numerical schemes are firstly validated on 1D and 2D problems involving compressible two-phase flows [39], [4]. Particularly, a good
<table><caption>Table 1: Air-helium EOS parameters and post-shock conditions.</caption><thead><tr><th></th><th>γ</th><th>P<sub>∞</sub></th><th>C<sub>p</sub></th></tr></thead><tbody><tr><td>air</td><td>1.4</td><td>0 Pa</td><td>1004.5 J/K.kg</td></tr><tr><td>helium</td><td>1.648</td><td>0 Pa</td><td>3955.2 J/K.kg</td></tr></tbody><thead><tr><th></th><th>P</th><th>ρ</th><th>u</th></tr></thead><tbody><tr><td>post-shock</td><td>1.444 × 10<sup>5</sup> Pa</td><td>1.51 kg/m<sup>3</sup></td><td>93.65 m/s</td></tr></tbody></table>

performance of HLLC present variant has been demonstrated in [4], while validation test cases in [39] have shown good agreement with reference solution in results obtained by using MUSCL-Hancock scheme coupled with higher order reconstruction methods, WENO methods, for example. These test cases included the water-gas mixture shock tube and vacuum-like problems, e.g. water-gas mixture expansion shock tube with stronger discontinuity in velocity. In order to evaluate the performance of considered numerical scheme to compute efficiently and accurately a bubble collapse in a liquid medium we propose two validation test cases compared with experimental results (2D air-helium shock-bubble interaction) and previous publications (2D bubble collapse in a free field and next to a wall).

### 4.1. 2D air-helium shock-bubble interaction

The validation of the solver is performed by studying the 2D problem of an air-helium shock-bubble interaction. This problem is chosen due to the existence of experimental studies by [1] including visualizations at different times. The helium bubble has an initial diameter of 4 cm and is impacted by a normal shock wave moving at the Mach number $M_{sh}=1.175$. The EOS parameters and post-shock conditions are presented in Table 1. The volume fraction $\alpha$ represents here the ratio of the lighter gas to the carrier gas.

The computations are performed on a half-domain due to the symmetry of the problem. Unless stated otherwise, the results have been obtained by using
<figure><img src="image_4.png" /><figcaption>Figure 1: Experimental visualization by Haas and Sturtevant, extracted from [1] (top), numerical results computed on 4000×400 nodes, $\Delta t=2.5 \times 10^{-9}$ s, HLLC Hancock WENO5-IS (bottom).</figcaption></figure>

an uniform mesh discretization of 4000 × 400 cells, a time step of $2.5 \times 10^{-9}$ s and by using the HLLC flux formulation.

First, the density gradient evolution compared to the experimental visualizations is shown in Figure 1 for simulations performed using the WENO5-IS reconstruction. The bubble deformation starts by its flattening in the direction of the shock wave and eventually takes the shape of the kidney as a result of the high speed air jet formation located at the upstream interface. On the other hand, the downstream interface jet leads to the appearance of the counter-rotating vortices which cause the bubble elongation. The computations obtained with all schemes (MUSCL, WENO3-Z and WENO5-IS) are able to reproduce the main elements of the bubble shape evolution compared to the experimental data.

The comparison of three numerical methods performance is now proposed on Figure 2. Certainly, the qualitative evaluation demonstrates more diffused result computed by using MUSCL reconstruction. However, the evaluation of WENO3-Z and WENO5-IS shows that the reconstruction is fairly similar between two schemes with WENO5-IS having more structure on the bubble surface and interanl vortices.

Lastly, the effect of the characteristic variables reconstruction has been analysed in [39]. No discrepancy on the results by using characteristic or primitive
<figure><img src="image_5.png" /><figcaption>Figure 2: Air-helium shock-bubble interaction, comparison of the numerical schemes. Mesh: 4000×400, $\Delta t=2.5 \times 10^{-9}$ s, $t_{final}=0.05$ ms, left to right: MUSCL, WENO3-Z, WENO5-IS.</figcaption></figure>

variables have been observed.

The study of the mesh convergence of the schemes has been performed on the maximum pressure component computed at the final time $t_{final}=0.05$ ms with CFL=0.3 by using MUSCL, WENO3-Z and WENO5-IS reconstruction methods. The evolution of the maximum pressure reached is plotted in Figure 3 (left part). The solution of the maximum pressure has three main peaks. The first one of the approximate value 1.8 bar occurs at time 0.003 ms. WENO schemes led to the similar pressure values while MUSCL has a slightly lower pressure intensity. The second and highest peak of approximately 2.2 bar happens at time 0.012 ms, which has been reconstructed similarly by all schemes. Finally, the third peak of 1.9 bar with two parts occurring at times between 0.029-0.035 ms has been once again computed similarly by WENO numerical schemes and lower peak intensity obtained by MUSCL reconstruction.

The normalized pressure value of the highest peak (estimated by using the solution computed on the finest mesh 4000×400 nodes with the WENO5-IS formulation) is shown in Figure 3 (right part) against the gradually increased mesh resolution indicated by the number of points per bubble diameter (ppbd). The numerical rate of convergence of this value (similar to [40]) of approximately second order is verified by computing the relative error compared with the finest
<figure><img src="image_6.png" /><figcaption>Figure 3: Air-helium shock-bubble interaction, maximum pressure comparison, mesh: 4000×400 (left), convergence of the maximum pressure (right). CFL=0.3</figcaption></figure>

mesh pressure value. The second order is consistent with basic reconstruction
schemes employed in present work which is the MUSCL-Hancock reconstruction
with HLLC solver. It is important to note that the problems involving shocks
cannot demonstrate the order of convergence higher than one. This is the case
for third and even fifth order numerical schemes. In order to demonstrate a
proper numerical order of convergence, a test case with smooth solution is re-
quired, which is out of scope of this paper. That is a reason why the rate of
convergence is analysed instead of order of convergence.

4.2. 2D Shock-induced bubble collapse in a free field

The present test is based on one presented in [41] and computed by various
authors (for example [42]). An air bubble with diameter 6 mm, immersed in
water at rest, is impacted by a normal shock wave for which the Mach number
$M_{sh}=1.7$. The physical conditions are initialised as following: $P = 10^5$ Pa, $\rho_{air} = 1$ kg/m^3 and $\rho_{water}$ 1000 kg/m^3. Only a half of the bubble is considered due to the symmetry of the problem. The domain has top and bottom boundaries set as slip walls, left and right boundaries have imposed non-reflective conditions.
The computations are made with spatial discretization of 2000×1000 cells using either a fixed $\Delta t = 10^{-9}$ s or a fixed CFL number. The EOS parameters and post-shock condition are stated in Table 2.

Firstly, a description of the main physical phenomena is proposed in Figure
<table><caption>Table 2: Water-air EOS parameters and post-shock condition</caption><thead><tr><th></th><th>γ</th><th>P<sub>∞</sub></th><th>C<sub>p</sub></th></tr></thead><tbody><tr><td>air</td><td>1.4</td><td>0 Pa</td><td>4180 J/K.kg</td></tr><tr><td>water</td><td>4.4</td><td>6 × 10<sup>8</sup> Pa</td><td>1487 J/K.kg</td></tr><tr><th></th><th>P</th><th>ρ</th><th>u</th></tr><tr><td>post-shock</td><td>1.9 × 10<sup>9</sup> Pa</td><td>1323.65 kg/m<sup>3</sup></td><td>681.58 m/s</td></tr></tbody></table>

where the density gradient modulus (on the top) and the pressure field (on the bottom) obtained from the WENO5-IS solution are plotted. The obtained solution is in agreement with previous similar computations in [4], [23]. The collision of the shock with the bubble causes a strong rarefaction wave to be reflected in the backward direction and a weak shock wave to penetrate the bubble (at t=2 µs). The pressure difference leads to the bubble deformation which becomes kidney shaped (at time t=3 µs). This particular shape is a result of the vorticity occurred at the bubble edge due to the water jet caused by the passing shock wave. The high pressure zone is generated once this water jet hits the standing still water in front of the bubble causing the emission of a strong blast wave (time t=4 µs). Moreover, further propagation of the water jet and its high speed makes the blast front expand asymmetrically. Finally, a high pressure area is generated at t=4.8 µs due to the recollapse of the bubble pieces by the leftward front of the blast wave. While the blast wave continues to expand, the cavity shrinks and the low pressure fields are observed at the vortices core.

The time and the intensity of the maximum pressure peaks reached during the collapse are analysed. The mesh convergence is attested through the evolution of the highest pressure peak. Figure 5 (left) presents the evolution of the maximum pressure. This curve presents three main peaks (see Table 3). The first peak occurs at approximately 3.6 µs due to the generation of the blast wave. Its value is around 50,000 bar, the highest being predicted by the WENO5-IS,
<figure><img src="image_7.png" /><figcaption>Figure 4: Bubble collapse time evolution, (a) 2 µs, (b) 3 µs, (c) 4 µs, (d) 4.8 µs, HLLC Hancock WENO5-IS, mesh 2000 × 1000, CFL=0.3.</figcaption></figure>

while the lowest intensity is reproduced by the MUSCL reconstruction. The second peak happens at approximate time 4.4 µs due to the generation of another strong shockwave in the liquid when the expelled-gas shock focuses on the bubble's most upstream point. Finally, the third and strongest peak is observed at approximate time 4.7-4.8 µs due to the recollapse of the bubble fragments. This peak has the most deviations in results

The convergence of the schemes is studied by refining the mesh and examining the third peak pressure intensity. Figure 5 (right) presents the convergence by using the second order MUSCL reconstruction and higher order methods, WENO3-Z and WENO5-IS. The computations have been performed by setting CFL=0.1. The oscillating nature of the solution is noted across all schemes which makes it problematic to perform the numerical convergence study. Thus,
<table><caption>Table 3: Shock-induced bubble collapse, comparison of the maximum pressure peaks P<sub>max</sub>.</caption><thead><tr><th colspan="7">Water-air shock-bubble, Pmax (bar)</th></tr><tr><th>Scheme</th><th>t1</th><th>P1</th><th>t2</th><th>P2</th><th>t3</th><th>P3</th></tr></thead><tbody><tr><td>MUSCL</td><td>3.6µs</td><td>49,910</td><td>4.4µs</td><td>58,100</td><td>4.8µs</td><td>86,630</td></tr><tr><td>WENO3-Z</td><td>3.6µs</td><td>50,590</td><td>4.4µs</td><td>59,620</td><td>4.8µs</td><td>86,280</td></tr><tr><td>WENO5-IS</td><td>3.6µs</td><td>50,730</td><td>4.4µs</td><td>67,290</td><td>4.8µs</td><td>89,330</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 5: Shock-induced bubble collapse, maximum pressure, mesh 2000 × 1000 (left), convergence of the third pressure peak (right), CFL=0.1.</figcaption></figure>

in order to perform approximate convergence study, the Savitsky-Golay filter of fourth order with 100 points is used to smooth the data [43]. The better convergence of high-order WENO3-Z and WENO5-IS is noted compared with second-order MUSCL scheme where slight discrepancy is noted. The mesh independent solution is achieved at around 500 points per bubble diameter, which corresponds to the mesh 2000×1000 nodes.

The density gradient along with the pressure profile evolution at the time of the third pressure peak $t_3$ is proposed in Figure 6. The presented results have been computed on the mesh 2000×1000 with CFL=0.3. The variation in the solution obtained by using different order schemes is noticed. The most detailed solution is obtained by using WENO5-IS and WENO3-Z schemes where all the details of the bubble core structure can be distinguished. However, the least amount of details is produced by the second order MUSCL approach, which is a consistent result based on the theoretical order of accuracy of these schemes.
<figure><img src="image_9.png" /><figcaption>Figure 6: Shock-induced bubble collapse, contour of density gradient (top) and pressure field (bottom), MUSCL (left), WENO3-Z (middle) and WENO5-IS (right). Mesh 2000 × 1000, CFL=0.3, tpeak3.</figcaption></figure>

The pressure profiles are illustrated at the bottom of the plots. The low and high pressure areas are present on the plots of all three schemes and are similar to each other.

### 4.3. 2D Bubble collapse near a wall

The present test has been proposed by Paquette et al. [13] and then computed in 2D and in 3D in [23] by using MUSCL-Hancok. A wall is placed behind the bubble. A normal shock wave initially located at 0.1 mm is moving with a Mach number $M_{sh}=1.09$ causing the bubble collapse. The EOS parameters stay the same and the post shock conditions are given below:

$$ \begin{pmatrix} P \\ \rho \\ u \end{pmatrix}_{\text{post-shock}} = \begin{pmatrix} 1.2 \times 10^8 \text{ Pa} \\ 1049 \text{ kg/m}^3 \\ 75 \text{ m/s} \end{pmatrix} \qquad (34) $$

The mesh convergence is analysed on the evolution of the maximum pressure peaks. Figure 7 presents the maximum pressure evolution next to a wall (top left part) and mesh convergence of the MUSCL, WENO3-Z and WENO5-IS schemes. The intensity of the peak depending on the employed numerical scheme is presented in Table 4. The mesh is ranged between 750×1250 and 1908×3180 (250 to 650 ppbd). Three main peaks in the pressure solution are observed. The first one, with the lowest intensity, occurs at around 0.13 µs and is due to the impact of the incident shock wave on the wall. The second one occurs at
<table><caption>Table 4: Bubble collapse near a wall, comparison of the maximum pressure peaks $P_{max}$.</caption><thead><tr><th colspan="7">Water-air shock-bubble with a wall, maximum pressure $P_{max}$ (bar)</th></tr><tr><th>Scheme</th><th>$t_1$</th><th>$P_1$</th><th>$t_2$</th><th>$P_2$</th><th>$t_3$</th><th>$P_3$</th></tr></thead><tbody><tr><td>MUSCL</td><td>0.13µs</td><td>2,573</td><td>0.28µs</td><td>9,444</td><td>0.36µs</td><td>8,829</td></tr><tr><td>weno3-Z</td><td>0.13µs</td><td>2,565</td><td>0.28µs</td><td>9,453</td><td>0.36µs</td><td>8,938</td></tr><tr><td>weno5-IS</td><td>0.13µs</td><td>2,563</td><td>0.28µs</td><td>9,456</td><td>0.35µs</td><td>9,138</td></tr></tbody></table>

0.28 µs and is due to the impact of the generated blast wave on the wall. The third one at 0.36 µs is due to the impact of a secondary wave emitted by the recollapse of the bubble. All three numerical schemes lead to approximately similar time for these peaks with, however, the highest third pressure peak being obtained by computing the solution with the schemes of higher order, i.e. WENO5-IS and WENO3-Z, while the first and second pressure peaks having very small difference caused by the numerical scheme. The tests have shown the convergence of all schemes with a mesh independent solution achieved using 500 points per bubble diameter.

Next we propose the comparison of the solution based on the density gradients and pressure fields. The computations are made with fixed CFL=0.1 and 650 points per bubble diameter, which corresponds to a spatial discretization of 1908×3180 nodes. All obtained solutions are in close agreement with computations presented in [23]. The solution improvement based on the enlargement of the density gradient is illustrated. The gradient of the density with corresponded pressure profile at time of the third pressure peak is proposed in Figure 8. For comparison purposes five main elements are distinguished in the bubble vortices curvature, i.e. c1, c2, c3, c4 and c5. The numerical solution obtained by MUSCL method is similar to the computations in [23], however, does not reconstruct the vortices as accurately as compared with WENO5-IS and WENO3-Z. All the elements of curvature have been smoothed out and, moreover, the element c4 has a different shape and almost horizontal angle which might correspond to the delayed occurrence of the physical phenomena.
<figure><img src="image_10.png" /><figcaption>Figure 7: Bubble collapse next to a wall, maximum wall pressure obtained with the mesh 1908 × 3180, CFL=0.1 (top left), convergence of the maximum pressure peaks Pmax obtained with MUSCL (top right), WENO3-Z (bottom left), WENO5-IS (bottom right), CFL=0.1 (right).</figcaption></figure>

The pressure profile has also considerably less area of high pressure values. We believe that higher order schemes WENO have led to the improvement of the solution with, however, not big difference between WENO3-Z and WENO5-IS.

The numerical solution assessment is concluded by analysing the CPU cost. Indeed, an increase of the CPU cost is present as the order of the method improves, in particular due to the increase of the computational stencil size. The computational time of the solutions obtained by using the key numerical schemes discussed in (i.e. MUSCL, WENO3-Z, WENO5-IS) are proposed in Table 5. Although the solution obtained by high-order numerical methods has been improved, the cost has been increased by factor 1.3 between MUSCL and WENO3-Z and 1.4 between MUSCL and WENO5-IS schemes. This cost of computations is considerably more critical for higher resolution computations and, particularly, for 3D problems. Thus, given only a small differences in the qual-
<figure><img src="image_11.png" /><figcaption>Figure 8: Bubble collapse next to a wall, comparison of the reconstruction method, contour of density gradient (top), pressure field (bottom), mesh: 1908×3180, CFL=0.1, WENO5-IS (left), WENO3-Z (middle) and MUSCL (right).</figcaption></figure>

<table><caption>Table 5: CPU cost, bubble collapse near a wall. Mesh: 1908×3180. CFL=0.1.</caption><thead><tr><td>Scheme</td><td>CPU (h)</td><td>ratio (%)</td></tr></thead><tbody><tr><td>MUSCL</td><td>94.7</td><td>1</td></tr><tr><td>WENO3-Z</td><td>121.6</td><td>1.28</td></tr><tr><td>WENO5-IS</td><td>133.7</td><td>1.41</td></tr></tbody></table>

itative and quantitative examination of results by WENO3-Z and WENO5-IS, the choice of third order WENO3-Z for 3D computations is motivated. However, given the estimated number of points per bubble diameter necessary to obtain the converged results, the stretching mesh strategy is required to reduce the CPU cost in 3D. The stretching mesh function details are described in Appendix B.

## 5. Computations of the spherical bubble collapse

This section presents the results of the shock-induced bubble collapse near a wall in fully parallelised solver, SCB with introduction of stretched mesh. The numerical methods have been validated by using academic 1D and 2D problems where an appropriate spatial-temporal resolution has been established. Moreover, the variety of stretching techniques have been analysed in [39]. The
novelty of this study is to extend 3D computations in [23] to non-uniform mesh
and improve the spacial resolution by using WENO3-Z with a goal to make
computations considerably more affordable and the results improved thanks to
the third order spacial accuracy.

The section starts with a brief presentation of the SCB solver. The following
subsections are dedicated to 3D results of shock bubble collapse near a wall with
an introduction of the mesh stretching strategy.

5.1. SCB solver

SCB is an efficient and simple parallel multi-phase solver developed to sim-
ulate various compressible multiphase flows, in particular bubble collapse. The
parallelisation of the solver is based on hybrid approach of using OpenMP and
MPI libraries. Recently the parallelisation strategy of combination of MPI and
OpenACC has been performed given large developments of CPU and GPU-
based supercomputers, see [44].

**The distributed memory parallelisation** is based on the MPI library.
The computational stencil of SCB relies on five points per direction, which
implies that each unknown is computed by using 13 neighbours. The HLLC
Riemann solver is parallelised by employing global arrays distribution between
the processes. A 3D block partitioning of the matrix is introduced in order
to perform a decomposition. This technique is based on the introduction of
"ghost" cells at each subdomain and data communication between neighbouring
cells of one subdomain. A Cartesian processor topology is used to organise the
subdomains. The computations are coupled by using two layers of auxiliary cells
which are defined on the boundaries of each subdomain. The better performance
can be achieved if the subdomains has equal size and square topology.

**Performance and scalability** A problem of shock-induced bubble collapse requires intensive computing power due to the high-speed dynamics and small spatio-temporal scales. The physical phenomena of this problem can be resolved only with very fine mesh, e.g. $10^9$ points in 3D. Additionally, the time period of
the high solution variation is very short and measured in microseconds. Moreover, the complexity of the mathematical model has to be taken into account. This means that explicit time integration has to be used.

**The non-uniform mesh** is generated by separately developed program, which takes as inputs initial location and size of space steps $dx$, $dy$ and $dz$. Since the mesh stretching is a problem-dependent, every strategy has to be implemented based on the test case in consideration. The non-uniform mesh is generated and written into the data file, which is then used by SCB. The distribution of the grid coordinates values is required for the parallel computations.

**An introduction of the high order schemes in SCB** requires a communication of additional neighbours per unknown. This is particularly the case for the schemes of fifth order where the high communication cost is questionable in terms of its proportionality to the solution improvement. At this stage WENO3-Z scheme has been introduced in the framework of OpenMP and MPI. The implementation of the method on the uniform mesh does not require an enlargement of the existing stencil. However, the new communication is required for the non-uniform mesh, where we have to share the values for the non-uniform mesh coefficients needed for the scheme. While this does not enlarge a number of neighbours for MUSCL method, it certainly does for WENO3 reformulated for the non-uniform grid, where two additional nodes have to be defined at the first time iteration to compute the coefficients. This can increase a communication cost. In order to avoid this, we note that due to the smooth mesh stretching function used in this study, the variation of coefficients values needed for the non-uniform mesh is very small. Thus, the coefficients computed for WENO3 at the boundaries of the domain are set equal to its neighbour.

## 5.2. Spherical bubble collapse near a wall

This problem is essentially the same as the previous one but extended to 3D.

The computational domain is of a size 0.3×0.5×0.5 mm and the initial location of the bubble is at (0.2,0) mm. Exploiting the symmetry of the problem,
only a quarter of the bubble is computed. The configuration at initial time is
presented on Figure 9. The computations are done by setting $\Delta t = 10^{-10}$ s.

<figure><img src="image_12.png" /><figcaption>Figure 9: Configuration at initial time, $P_{sh}=1200$ bar</figcaption></figure>

The previous 2D validation tests have shown the convergence of all schemes
with a mesh independent solution achieved by using at least 500 points per
bubble diameter. Additionally, a demonstration of fine mesh requirement in 3D
by using stretched mesh is proposed in Figure 12. Two key components of the
solution are examined, i.e. maximum pressure inside the fluid and on the wall.
The obtained maximum pressure curve next to a wall does not have a critical
variation at highest intensity peak. However, the pressure inside the fluid has
a gain of more than 30% in highest peak. This increase of intensity which is
a result of increased number of nodes in computational domain confirms the
requirement of high number of points per bubble diameter in order to obtain an
accurate reconstruction of physical phenomena.
<figure><img src="image_13.png" /><figcaption>Figure 10: Time evolution of pressure. Comparison between stronger and weaker stretching.
$P_{sh} = 1200$ bar</figcaption></figure>

The time evolution of maximum pressure inside the fluid and along the wall is presented on Figure 11. We introduce a 2D result computed on the corresponding uniform mesh for comparison. As noted in [23] the collapse of spherical bubble develops faster and its intensity is stronger. Indeed, the same observation is confirmed in our computations. The intensity of maximum pressure peak is of factor 5 compared to 2D solution and it occurs approximately 0.5µs faster. However, the wall pressure (the right part of the plot) has very similar intensity at maximum peak, as indeed has been observed in [23].

<figure><img src="image_14.png" /><figcaption>Figure 11: Time evolution of the maximum pressure inside the fluid (left) and near the wall (right). WENO3-Z. 2D, uniform mesh: 1500×2500; 3D, stretched mesh: 1050×850×850.
$P_{sh} = 1200$ bar</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 12: Time evolution of pressure. Comparison between coarse and fine mesh. $P_{sh} = 1200$ bar</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13: Visualisation of 3D bubble collapse at different times for a stand-off distance $L/R_0 = 2$. Dimensionless longitudinal velocity component $u/u_{sh}$, dimensionless wall pressure $P/P_{sh}$, Schlieren-like representation and isosurface of void ratio. MUSCL, $t=0.28$ µs. $P_{sh} = 1200$ bar</figcaption></figure>

It is important to establish the tolerance level for the stretching factor which affects the solution accuracy away from the uniform mesh area. This degree of stretching is examined by using MUSCL-Hancok 2nd order in space and time scheme due to its smaller computational stencil and the initial spacial interval corresponding to the mesh of 750×1250×1250 nodes. The uniform mesh box is
set similarly in all computations and is located around initial bubble diameter.
The 3D visualisations presented on Figure 13 illustrate the deviation of the
solution with increased stretching factor. This is especially the case for the wall
pressure component where the distribution values are not smoothly propagated
and have a square pattern. However, the solution around the bubble, where the
uniform mesh is set, remains accurate and similar in all tests. This emphasises
the need of proper definition of how much stretching can be applied to obtain the
solution of good quality. Indeed, the solution accuracy for the problems where
non-uniform mesh is used is dependent on the ratio between $dx_{max}/dx_{min}$,
$dy_{max}/dy_{min}$ and $dz_{max}/dz_{min}$. We limit our consideration on the directions
Y and Z since the uniform mesh area is relatively small for Y and Z axes. Saying
that, the ratio for the first stretched mesh configuration (550×900×900) is 1.7
and the obtained solution is visually of good quality, while the second stretched
mesh configuration (550×600×600) with approximate ratio equal to 3.2 yield to
the deformation of the pressure values distribution. Thus, the initialisation for
the finer stretched mesh computations are set according to the ratio between
maximum and minimum spacing in directions of y and z equal to 1.7.
<figure><img src="image_17.png" /><figcaption>Figure 14: Visualisation of 3D bubble collapse at different times for a stand-off distance $L/R_0 = 2$. Dimensionless longitudinal velocity component $u/u_{sh}$, dimensionless wall pressure $P/P_{sh}$, Schlieren-like representation and isosurface of void ratio. WENO3-Z. Comparison of stretched mesh factor. $P_{sh} = 1200$ bar</figcaption></figure>

An illustration of the bubble collapse evolution is proposed on Figure 15. The third order WENO3-Z scheme is employed in parallel framework in its form for the non-uniform mesh. The non-uniform mesh corresponding to 500 points per bubble diameter is generated and stretched in a way that the ratio between $dy_{max}/dy_{min}$ and $dz_{max}/dz_{min}$ does not exceed value of 1.7. The resulted stretched grid consists of 1050×1780×1780 points, with initialised spacing corresponding to 1500×2500×2500 points. We present the normalised longitudinal velocity by $u_{sh} = 75.07$ m/s on the vertical symmetry plane, normalised pressure by $P_{sh} = 1.2 \times 10^8$ Pa on the wall, the density gradient on the horizontal symmetry plane and void ratio inside the volume. The incident shock wave reaches the wall and reflected wave is generated. At time $t = 0.15\mu s$, the reflection wave is located near the bubble. We note the bubble deformation compared with its initial spherical shape at $t = 0.12\mu s$. This happens due to the difference in pressure between two sides. The water jet penetrating the bubble causes its toroid-like shape at $t = 0.18\mu s$ and the reflected shock wave has reached the bubble interface. The spherical intense blast wave is generated at
time $t = 0.20\mu s$ as the water-jet impacts the opposite side of bubble interface.
The high velocity values are generated and reach 1500 m/s. The occurrence of
strong pressure peak caused by the blast wave reaching the wall is observed at
$t = 0.23\mu s$. The reflected wave propagates in direction of toroid-like bubble.
Finally, this wave impacting the bubble leads to its recollapse at $t = 0.25\mu s$.
Consequently, another pressure peak is generated as can be observed on Figure
15.

Finally, once again we revisit the question about the stronger stretching
factor. Figure 14 presents side-by-side comparison of 3D visualisation at the
time of bubble recollapse under the incidence of reflections wave. A square-
shaped reconstruction of pressure profile is noted away from the bubble. It is
certainly due to the stretching factor since the computations on uniform mesh
and on stretched mesh with lower factor do not exhibit it. While the effect of
such strong stretching has been also seen in 2D computations [39], it did not
affect the pressure region around the bubble.

The comparison of maximum pressure inside the fluid and next to a wall is
proposed by using similar configuration of stretching mesh on Figure 10. The
highest and most critical peaks of pressure have not been affected by stronger
degree of stretching.

The CPU cost for present problem (Table 6) is recorded to be 6 hours for
the stretched mesh of 1050×850×850 and 34 hours for the stretched mesh of
1050×1780×1780 by using 4,000 processors. In comparison, the computational
time of similar problem in 3D with uniform coarser spacing of 750×1250×1250
and Δ*t*-value twice higher, is only 0.7 h less. We believe that given computa-
tional cost and solution accuracy inside the critical area of physical phenomena,
stronger mesh stretching (i.e. stretched mesh of 1050×850×850) can be used to
recover key characteristics of the bubble collapse process. Moreover, the stretch-
ing factor and the initialisation of non-uniform mesh can be adjusted depending
on the CPU cost appetite.
<table><caption>Table 6: CPU cost, 3D bubble collapse near a wall, WENO3-Z, stretched mesh. JeanZay
4,000 processors</caption><thead><tr><td>Problem</td><td>Total nb points</td><td>nppbd</td><td>CPU (h)</td></tr></thead><tbody><tr><td>P<sub>sh</sub> = 1200 bar</td><td>1050×1780×1780</td><td>500</td><td>34.9</td></tr><tr><td>P<sub>sh</sub> = 1200 bar</td><td>1050×850×850</td><td>500</td><td>5.9</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Figure 15: Visualisation of 3D bubble collapse at different times for a stand-off distance $L/R_0 = 2$. Dimensionless longitudinal velocity component $u/u_{sh}$, dimensionless wall pressure $P/P_{sh}$, Schlieren-like representation and isosurface of void ratio. WENO3-Z, stretched mesh: 1050×1780×1780</figcaption></figure>
## 6. Conclusion and perspectives

We have presented the 3D computations of the shock-bubble collapse near a wall. The non-uniform mesh based on hyperbolic tangent one-sided stretching function has been introduced for improvement of CPU time. The computational approach is based on a finite volume method with HLLC flux reconstruction with WENO3-Z spatial 3rd order accuracy coupled with Hancock 2nd order predictor-corrector scheme. These have been implemented into the parallel solver SCB. This solver relies on OpenMP and MPI libraries and provides an HPC framework for the efficient treatment of high-order stencil schemes as well as point-wise operations. We have validated our approach on 2D test cases in order to establish an appropriate computational strategy. By doing so, we have demonstrated the importance of very fine computational grid and, as a result, very quickly prohibitive computational cost. Eventually, we have investigated a highly-resolved shock-induced bubble collapse in vicinity to a wall, which has highlighted the high potential of the proposed approach and solver. Our new results are in close agreement with previously computed tests in [23] and provide a solid base for the future studies of pressure loads on solid materials with lower computational cost and higher spacial accuarcy.

An interesting perspective of the present work would be based on introduction of dynamic mesh adaptation techniques based on the location of physical phenomena. Indeed, the simple mesh stretching might have a significant computational cost increase if updated at every iteration. Solvers permitting, more advanced dynamic mesh refinement techniques might be beneficial in this case. These include, but are not limited to, adaptive mesh refinement method [14, 15] and Chimera [45]. Moreover, higher order methods relying on larger computational stencil, such as WENO5-IS mentioned before or even PPM reconstruction, could be examined. Additionally, the pressure loads study similar to those performed in [23] could be performed to examine the pressure loads differences with improved non-uniform mesh third order accurate solver.
Acknowledgments

This research was supported by the French National Research Agency ANR
(project 18-CE46-009).

Computations have been performed on the HPC resources of GENCI under
allocations A0112A1098.

Author contributions

All of the authors were involved in the preparation of the manuscript and
have read and approved the final manuscript version.

References

[1] J.-F. Haas, B. Sturtevant, Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities, J. Fluid Mech. 181 (1) (1987) 44–76.

[2] E. Johnsen, T. Colonius, Numerical simulations of non-spherical bub-
ble collapse, J. Fluid Mech. 629 (2009) 231–262. doi:10.1017/
S0022112009006351.

[3] S. Gong, E. Klaseboer, Interaction between collapsing bubble and viscoelastic solid: Numerical modelling and simulation, Applied Mathematical Modelling 40 (2016) 4746–4764.

[4] E. Goncalves, D. Zeidan, Simulation of compressible two-phase flows using a void ratio transport equation, Commun. Comput. Phys. 24 (1) (2018) 167–203. doi:10.4208/cicp.OA-2017-0024.

[5] M. Baer, J. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials, Int. J. Multi-phase Flow 12 (1986) 861–889. doi:10.1016/0301-9322(86)90033-9.

[6] A. Kapila, R. Menikoff, J. Bdzil, S. Son, D. Stewart, Two-phase mod-
elling of deflagration-to-detonation transition in granular materials: re-
duced equations, Phys. Fluids 13 (10) (2001) 3002–3024. doi:10.1063/1.
1398042.
[7] R. Saurel, F. Petitpas, R. Abgrall, Modelling phase transition in metastable liquids: application to cavitating and flashing flows, Journal of Fluid Mechanics 607 (2008) 313–350.

[8] E. Goncalves, P. Parnaudeau, Comparison of multiphase models for computing shock-induced bubble collapse, Int. J. Numer. Methods Heat Fluid Flow 30 (8) (2020) 3845–3877. doi:10.1108/HFF-05-2019-0399.

[9] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach, J. Comput. Phys. 125 (1) (1996) 150–160. doi:10.1006/jcph.1996.0085.

[10] E. Johnsen, T. Colonius, Implementation of WENO schemes in compressible multicomponent flow problems, J. Comput. Phys. 219 (2) (2006) 715–732. doi:10.1016/j.jcp.2006.04.018.

[11] V. Coralic, T. Colonius, Finite-volume WENO scheme for viscous compressible multicomponent flow, J. Comput. Phys. 274 (2014) 95–121. doi:10.1016/j.jcp.2014.06.003.

[12] B. Wang, G. Xiang, X. Hu, An incremental-stencil WENO reconstruction for simulation of compressible two-phase flows, Int. J. Multiphase Flow 104 (2018) 20–31. doi:10.1016/j.ijmultiphaseflow.2018.03.013.

[13] Y. Paquette, M. Fivel, G. Ghigliotti, E. Johnsen, J.-P. Franc, Fluid-structure interaction in cavitation erosion, in: 10th International Symposium on Cavitation CAV2018, Baltimore, USA, 2018.

[14] M. Berger, O. Joseph, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of Computational Physics 53 (3) (1984) 484–512.

[15] M. Berger, P. Colella, Local adaptive mesh refinement for shock hydrodynamics, Journal of Computational Physics 82 (1) (1989) 64–84.
[16] J. Thompson, Z. Warsi, C. Mastin, Numerical grid generation: foundations and applications, Elsevier Science Ltd, 1985.

[17] M. Vinokur, On one-dimensional stretching functions for finite-difference calculations, Journal of computational physics 50 (2) (1983) 215-234.

[18] B. Pierson, P. Kutler, Optimal nodal point distribution for improved accuracy in computational fluid dynamic, AIAA J. 18 (1) (1980) 49-54.

[19] D. Gough, E. Spiegel, J. Toomre, Highly stretched meshes as functionals of solutions, Proceedings of the fourth international conference on numerical methods in fluid dynamics (1975) 191-196.

[20] R. Wang, H. Feng, R. Spiteri, Observations on the fifth-order WENO method with non-uniform meshes, Appl. Math. Comput. 196 (2008) 433-447.

[21] M. Vinokur, On one-dimensional stretching functions for finite-difference calculations, NASA CR 3313.

[22] E. Goncalves, B. Charriére, Modelling for isothermal cavitation with a four-equation model, Int. J. Multiphase Flow 59 (2014) 54-72. doi:10.1016/j.ijmultiphaseflow.2013.10.015.

[23] E. Goncalves, P. Parnaureau, Numerical study of pressure loads generated by a shock-induced bubble collapse, Physics of fluids doi:10.1063/5.0069332.

[24] O. LeMétayer, J. Massoni, R. Saurel, Elaborating equations of state of a liquid and its vapor for two-phase flow models, Int. J. Therm. Sci. 43 (3) (2004) 265-276. doi:10.1016/j.ijthermalsci.2003.09.002.

[25] B. Charriére, Modelisation et simulation d'ecoulements turbulents cavitants avec un modèle de transport de taux de vide, PhD of Institut Polytechnique de Grenoble.
[26] A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, D. S. Stewart, Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations, Physics of Fluids 13 (10) (2001) 3002–3024.
doi:10.1063/1.1398042.

[27] C. Gouin, C. Junqueira-Junior, E. Goncalves, J.-C. Robinet, Numerical investigation of three-dimensional partial cavitation in a venturi geometry, Physics of Fluids 33(6) (2021) 063312.

[28] R. Abgrall, P. Bacigaluppi, S. Tokareva, A high-order nonconservative approach for hyperbolic equations in uid dynamics, Computers Fluids 169 (2018) 10–22.

[29] T. Gallouet, J. Herard, N. Seguin, A hybrid scheme to compute contact discontinuities in Euler systems, Math. Modell. Numer. Anal. 36 (6) (2002) 1133–1159.
doi:10.1051/m2an:2003009.

[30] K. Kitamura, M.-S. Liou, C.-H. Chang, Extension and comparative study of AUSM-family schemes for compressible multiphase flow simulations, Commun. Comput. Phys. 16 (3) (2014) 632–674.
doi:10.4208/cicp.020813.190214a.

[31] A. Jameson, W. Schmidt, E. Turkel, Numerical solution of the euler equations by finite volume methods using Runge-Kutta time stepping schemes, AIAA J. 81. doi:10.2514/6.1981-1259.

[32] E. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL Riemann solver, Shock Waves 4 (1) (1994) 25–34.

[33] F. Daude, P. Galon, Z. Gao, E. Blaud, Numerical experiments using a HLLC-type scheme with ALE formulation for compressible two-phase flows five-equation models with phase transition, Comput. Fluids 94 (2014) 112–138.
doi:10.1016/j.compfluid.2014.02.008.
[34] B. van Leer, Upwind and high-resolution methods for compressible flow: from donor cell to residual-distribution schemes, 16th AIAA Computational Fluid Dynamics Conference.

[35] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to godunov's method, J. Comput. Phys. 32 (1) (1979) 101-136. doi:10.1016/0021-9991(79)90145-1.

[36] J. Bell, P. Colella, J. Trangenstein, Higher order Godunov methods for general systems of hyperbolic conservation laws, J. Comput. Phys. 82 (1989) 362-397. doi:10.1016/0021-9991(89)90054-5.

[37] R. Borges, M. Carmona, B. Costa, W. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (2008) 3191-3211. doi:10.1016/j.jcp.2007.11.038.

[38] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1) (1996) 202-228. doi:10.1006/jcph.1996.0130.

[39] K. Kozhanova, High-order numerical methods for shock-bubble interaction computations. other, PhD in ISAE-ENSMA Ecole Nationale Supérieure de Mécanique et d'Aérotechique - Poitiers.

[40] G. L. Spina, M. Vitturi, High-resolution finite volume central schemes for a compressible two-phase model, SIAM Journal on Scientific Computing 34 (6) (2012) 861-80. doi:10.1137/12087089X.

[41] N. Bourne, J. Field, Shock-induced collapse of single cavities in liquids, Journal of Fluid Mechanics 244 (1992) 225-240. doi:10.1017/S0022112092003045.

[42] R. Nourgaliev, T. Dinh, T. Theofanous, Adaptive characteristics-based matching for compressible multifluid dynamics, J. Comput. Phys. 213 (2006) 500-529. doi:10.1016/j.jcp.2005.08.028.
[43] A. Savitzky, M. Golay, Smoothing and differentiation of data by simplified least squares procedures, Analytical Chemistry 36 (1964) 1627–1639.

[44] R. Dubois, E. Goncalves, P. Parnaudeau, High performance computing of stiff bubble collapse on cpu-gpu heterogeneous platform, Computers and Mathematics with Applications 99 (2021) 246–256.

[45] F. C. D. J. Benek, J. L. Steger, A flexible grid embedding technique with application to the euler equations, AIAA Paper (1983) 83–1944.

[46] J. Thompson, F. Thames, C. Mastin, Tomcat - a code for numerical generation of boundary-fitted curvilinear coordinate systems on fields containing any number of arbitrary two-dimensional bodies, J. Comp. Phys 24 (3) (1977) 274–302.

## Appendix A. Coefficients for WENO computed on non-uniform mesh

### Spacial stencil

$$ \hat{b}_{22} = \frac{1}{h_1 + h_2 + h_3} + \frac{1}{h_2 + h_3} + \frac{1}{h_3}; \quad \hat{b}_{21} = \hat{b}_{22} - \frac{(h_1 + h_2 + h_3)(h_2 + h_3)}{(h_1 + h_2)h_3h_2} \quad (A.1) $$

$$ \hat{b}_{20} = \hat{b}_{21} + \frac{(h_1 + h_2 + h_3)h_3}{(h_2 + h_3)h_2h_1}; \quad \hat{b}_{12} = \frac{(h_2 + h_3)h_3}{(h_2 + h_3 + h_4)(h_3 + h_4)h_4} \quad (A.2) $$

$$ \hat{b}_{11} = \hat{b}_{12} - \frac{1}{h_2 + h_3} + \frac{1}{h_3} - \frac{1}{h_4}; \quad \hat{b}_{10} = \hat{b}_{11} - \frac{(h_2 + h_3)h_4}{(h_3 + h_4)h_2h_3} \quad (A.3) $$

$$ \hat{b}_{02} = - \frac{h_4 h_3}{(h_3 + h_4 + h_5)(h_4 + h_5)h_5}; \quad \hat{b}_{01} = \hat{b}_{02} + \frac{h_3(h_4 + h_5)}{(h_3 + h_4)h_4h_5} \quad (A.4) $$

$$ \hat{b}_{00} = \hat{b}_{01} - \frac{1}{h_3} - \frac{1}{h_4} - \frac{1}{h_4 + h_5} \quad (A.5) $$

For the case when r = 0 and j = 0, 1, 2.

$$ \hat{b}_{02} = \frac{h_3(h_4 + h_3)}{(h_3 + h_4 + h_5)(h_4 + h_5)h_5}; \quad \hat{b}_{01} = \hat{b}_{02} - \frac{h_3(h_3 + h_4 + h_5)}{(h_3 + h_4)h_4h_5} \quad (A.6) $$

$$ \hat{b}_{00} = \hat{b}_{01} + \frac{(h_3 + h_4)(h_3 + h_4 + h_5)}{h_3h_4(h_4 + h_5)} \quad (A.7) $$
Smoothness measure.

$$B_{r2}'' = \frac{6}{(h_{3-r} + h_{4-r} + h_{5-r})(h_{4-r} + h_{5-r})h_{5-r}} \quad (A.8)$$

$$B_{r1}'' = B_{r2}'' - \frac{6}{(h_{3-r} + h_{4-r})h_{4-r}h_{5-r}} \quad (A.9)$$

$$B_{r0}'' = B_{r1}'' + \frac{6}{(h_{4-r} + h_{5-r})h_{3-r}h_{4-r}} \quad (A.10)$$

$$B'_{22}(x_{i-\frac{1}{2}}) = \frac{2(h_1 + 2h_2)}{(h_1 + h_2 + h_3)(h_2 + h_3)h_3}; \qquad B'_{21}(x_{i-\frac{1}{2}}) = B'_{22}(x_{i-\frac{1}{2}}) - \frac{2(h_1 + 2h_2 - h_3)}{(h_1 + h_2)h_2h_3} \quad (A.11)$$

$$B'_{20}(x_{i-\frac{1}{2}}) = B'_{21}(x_{i-\frac{1}{2}}) + \frac{2(h_1 + 2h_2 - h_3)}{(h_2 + h_3)h_1h_2}; \qquad B'_{12}(x_{i-\frac{1}{2}}) = \frac{2(h_2 - h_3)}{(h_2 + h_3 + h_4)(h_3 + h_4)h_4} \quad (A.12)$$

$$B'_{11}(x_{i-\frac{1}{2}}) = B'_{12}(x_{i-\frac{1}{2}}) - \frac{2(h_2 - h_3 - h_4)}{(h_2 + h_3)h_3h_4}; \qquad B'_{10}(x_{i-\frac{1}{2}}) = B'_{11}(x_{i-\frac{1}{2}}) + \frac{2(h_2 - 2h_3 - h_4)}{(h_3 + h_4)h_2h_3} \quad (A.13)$$

$$B'_{02}(x_{i-\frac{1}{2}}) = -\frac{4(h_3 + h_4)}{(h_3 + h_4 + h_5)(h_4 + h_5)h_5}; \qquad B'_{01}(x_{i-\frac{1}{2}}) = B'_{02}(x_{i-\frac{1}{2}}) + \frac{2(2h_3 + h_4 + h_5)}{(h_3 + h_4)h_4h_5} \quad (A.14)$$

$$B'_{00}(x_{i-\frac{1}{2}}) = B'_{01}(x_{i-\frac{1}{2}}) - \frac{2(2h_3 - 2h_4 - h_5)}{(h_4 + h_5)h_3h_4} \quad (A.15)$$

$$B'_{r2}(x_i) = B'_{r2}(x_{i-\frac{1}{2}}) + \frac{1}{2}h_3 B''_{r2}; \qquad B'_{r2}(x_{i+\frac{1}{2}}) = B'_{r2}(x_{i-\frac{1}{2}}) + h_3 B''_{r2} \quad (A.16)$$

$$B'_{r1}(x_i) = B'_{r1}(x_{i-\frac{1}{2}}) + \frac{1}{2}h_3 B''_{r1}; \qquad B'_{r1}(x_i) = B'_{r1}(x_{i-\frac{1}{2}}) + h_3 B''_{r1} \quad (A.17)$$

$$B'_{r0}(x_i) = B'_{r0}(x_{i-\frac{1}{2}}) + \frac{1}{2}h_3 B''_{r0}; \qquad B'_{r0}(x_i) = B'_{r0}(x_{i-\frac{1}{2}}) + h_3 B''_{r0} \quad (A.18)$$
Weights coefficients.

$$d_2 = \frac{(h_3 + h_4)(h_3 + h_4 + h_5)}{(h_1 + h_2 + h_3 + h_4)(h_1 + h_2 + h_3 + h_4 + h_5)} \qquad (A.19)$$

$$d_1 = \frac{(h_1 + h_2)(h_3 + h_4 + h_5)(h_1 + 2h_2 + 2h_3 + 2h_4 + h_5)}{(h_1 + h_2 + h_3 + h_4)(h_2 + h_3 + h_4 + h_5)(h_1 + h_2 + h_3 + h_4 + h_5)} \qquad (A.20)$$

$$d_0 = \frac{(h_1 + h_2)h_2}{(h_2 + h_3 + h_4 + h_5)(h_1 + h_2 + h_3 + h_4 + h_5)} \qquad (A.21)$$

$$\hat{d}_2 = \frac{(h_4 + h_5)h_4}{(h_1 + h_2 + h_3 + h_4)(h_1 + h_2 + h_3 + h_4 + h_5)} \qquad (A.22)$$

$$\hat{d}_1 = \frac{(h_1 + h_2 + h_3)(h_4 + h_5)(h_1 + 2h_2 + 2h_3 + 2h_4 + h_5)}{(h_1 + h_2 + h_3 + h_4)(h_2 + h_3 + h_4 + h_5)(h_1 + h_2 + h_3 + h_4 + h_5)} \qquad (A.23)$$

$$\hat{d}_0 = \frac{(h_2 + h_3)(h_1 + h_2 + h_3)}{(h_2 + h_3 + h_4 + h_5)(h_1 + h_2 + h_3 + h_4 + h_5)} \qquad (A.24)$$

Appendix B. Non-uniform mesh generation

The generation of the non-uniform mesh follows the idea of [17] which is based on the stretched mesh by using one-sided function. The detailed analysis of appropriate stretching function can be found in [46]. Present study uses a hyperbolic one-sided tangent function which is briefly sketched below.

We define the left and right end of the function A and B, respectively. The general form of the function is $\epsilon(t, s_0, s_1)$, where the following is used:

* $t$ is a normalised parameter $t = \frac{\bar{t}-\bar{t}_A}{\bar{t}_B-\bar{t}_A}$ and

* $s_0$ and $s_1$ are the dimensionless slopes at A and B, respectively,

$$s_0 = \frac{d\epsilon}{dt}(0) \quad \text{and} \quad s_1 = \frac{d\epsilon}{dt}(1). \qquad (B.1)$$

The requirements for the stretching function are:

1. the order of all inverse scales of the $t$-variation with respect to $\epsilon$ should not be higher than one

2. the order of the slope $dt/d\epsilon$ must be equal to the order of minimum length scale of function of unknowns $\phi$ or vector function $r$-variations with respect to $t$ in the region where a rapid variation is located
These conditions are to be satisfied for all $s_0$ and $s_1$. Additionally, the required stretching function should be monotonic and preferably simple, easily invertible and continuous for all $s_0$, $s_1$. Lastly, this function to be constructed as such that $\bar{\epsilon}(\bar{t})$ is independent of the choice of the ends of the curve, A and B, when not normalised.

The following is a one-sided stretching function with zero slope at $t = 1$ and given initial slope at $t = 0$. There are three cases which have to be taken into account.

• if $s_0 > 1$,

$$s_0 = \frac{\sinh(2\Delta y)}{2\Delta y} \quad \text{and} \quad t = 1 + \frac{\tanh(\Delta y(\varepsilon - 1))}{\tanh(\Delta y)} \quad (B.2)$$

with the inverted equation as,

$$\varepsilon = 1 + \frac{\tanh^{-1}((t-1)\tanh(\Delta y))}{\Delta y} \quad (B.3)$$

• if $s_0 < 1$,

$$s_0 = \frac{\sin(2\Delta x)}{2\Delta x} \quad \text{and} \quad t = 1 + \frac{\tan(\Delta x(\varepsilon - 1))}{\tan(\Delta x)} \quad (B.4)$$

with the inverted equation as,

$$\varepsilon = 1 + \frac{\tan^{-1}((t-1)\tan(\Delta x))}{\Delta x} \quad (B.5)$$

• if $s_0 \asymp 0$. The following is an approximate formulations for $t$ and its inversion,

$$t = \varepsilon(1 - 0.5(s_0 - 1)(1 - \varepsilon)(2 - \varepsilon)) \quad (B.6)$$

and

$$\varepsilon = t(1 + 0.5(s_0 - 1)(1 - t)(2 - t)) \quad (B.7)$$
In order to incorporate the above developments, the solution of the non-linear equation has to be computed $y = sinh(x)/x$ for the case when $s_0 > 1$. This solution can be obtained by either solving this equation with one of the available iterative methods or by using approximate analytical formulation. The approximate analytical expression is presented in [17]. This analytical solution leads to a negligible numerical error and, thus, the resulting mesh discontinuity does not have an effect on the final solution. This solution is given in Appendix C.

## Appendix C. Approximate analytical solution to equation $y = sinh(x)/x$

We define $\bar{y} = 1 - y$, $v = \log(y)$ and $w = 1/y - 0.028527431$.

For $y < 2.7829681$:

$$x = \sqrt{6\bar{y}(1 - 0.15\bar{y} + 0.057321429\bar{y}^2 - 0.024907295\bar{y}^3 + 0.0077424461\bar{y}^4 - 0.0010794123\bar{y}^5)} \quad (C.1)$$

For $y > 2.7829681$:

$$x = v + (1 + 1/v)\log(2v) - 0.02041793 + 0.24902722w + 1.9496443 * w^2 - 2.6294547w^3 + 8.56795911w^4 \quad (C.2)$$

The one-sided stretching function requires the inversion of

$$y = \sin(x)/x \text{ for } s_0 < 1.$$

We denote $\bar{y} = 1 - y$ and an approximate analytical representation is following.

For $y < 0.26938972$,

$$x = \pi(1 - y + y^2 - (1 + \pi/6)y^3 + 6.794732y^4 - 13.205501y^5 + 11.72609y^6) \quad (C.3)$$
For $y < 0.26938972$,

$$x = \pi(1 - y + y^2 - (1 + \pi/6)y^3 + 6.794732y^4 - 13.205501y^5 + 11.72609y^6) \quad (C.4)$$

For $0.26938972 < y < 1$,

$$x = \sqrt{6\bar{y}}(1 + 0.15\bar{y} + 0.057321429\bar{y}^2 + 0.048774238\bar{y}^3 - 0.053337753\bar{y}^4 + 0.07584513\bar{y}^5) \quad (C.5)$$