# Numerical Simulation of Micro-Bubbles Dispersion by Surface Waves
## 3 Numerical Method  
In order to avoid coping with a strong geometric nonlinearity caused by the wavy upper boundary during the integration of the governing equations discussed above, a mapping is introduced transforming the domain with a wavy upper boundary into a domain with a flat upper boundary, and relating the Cartesian coordinates \((x,z)\) to curvilinear coordinates \((\xi,\eta)\) as:  
\[ \begin{array}{l}\xi=x+a\exp(kz)\sin k(x-ct)\\ \eta=z+a\exp(kz)\cos k(x-ct)\end{array}  \tag{9}\]  
Mapping (9) transforms the wavy (upper) boundary at \(z=z_{\rm s}(x,t)\) into a flat-plane boundary at \(\eta=0\). The water surface elevation, \(z_{\rm s}(x,t)\), is defined implicitly by Equation (9) and up to the second order in \(ka\) coincides with the Stokes-wave solution [25]:  
\[z_{\rm s}(x,t)\approx-a\cos k(x-ct)+\frac{a^{2}k}{2}\cos 2k(x-ct). \tag{10}\]  
An additional mapping is also employed for the vertical coordinate in the form:  
\[\widetilde{\eta}=\tanh\eta, \tag{11}\]  
where \(-1\leq\widetilde{\eta}\leq 0\), so that the grid nodes are clustered in the vicinity of the upper boundary, at \(\widetilde{\eta}=\eta=0\), and stretched with increasing depth. Since the mapping, Equation (9), is conformal, the following relations between the derivatives hold:  
\[\frac{\partial\xi}{\partial x}=J\frac{\partial x}{\partial\xi}=J\frac{\partial z }{\partial\eta}=\frac{\partial\eta}{\partial z};\ \ \frac{\partial\xi}{\partial z}=-J\frac{\partial x}{ \partial\eta}=J\frac{\partial z}{\partial\xi}=-\frac{\partial\eta}{\partial x}, \tag{12}\]  
where the Jacobian of the transformation is:  
\[J=\left(\frac{\partial\xi}{\partial x}\right)^{2}+\left(\frac{\partial\xi}{ \partial z}\right)^{2}. \tag{13}\]  
Due to the properties in Equations (12) and (13), the derivatives over the Cartesian coordinates, \(x\) and \(z\), can be related to the derivatives over curvilinear coordinates, \(\xi,\eta\), as,  
\[ \begin{array}{l}\frac{\partial}{\partial x}=\frac{\partial\xi}{\partial x} \frac{\partial}{\partial\xi}+\frac{\partial\eta}{\partial x}\frac{\partial}{ \partial\eta}=J\left(\frac{\partial x}{\partial\xi}\frac{\partial}{\partial \xi}+\frac{\partial x}{\partial\eta}\frac{\partial}{\partial\eta}\right),\\ \frac{\partial}{\partial z}=\frac{\partial\xi}{\partial x}\frac{\partial}{ \partial\xi}+\frac{\partial\eta}{\partial x}\frac{\partial}{\partial\eta}=J \left(\frac{\partial z}{\partial\xi}\frac{\partial}{\partial\xi}+\frac{ \partial z}{\partial\eta}\frac{\partial}{\partial\eta}\right).\end{array}  \tag{14}\]  
The Laplacian operator is also rewritten as:  
\[\frac{\partial^{2}}{\partial x_{j}\partial x_{j}}=I\left(\frac{\partial^{2}}{ \partial\xi^{2}}+\frac{\partial^{2}}{\partial\eta^{2}}\right)+\frac{\partial^ {2}}{\partial y^{2}} \tag{15}\]  
Equations (1) and (3) for the water velocity are discretized on a staggered grid consisting of \(360\times 240\times 180\) nodes, in the \(\xi\), \(y\) and \(\widetilde{\eta}\) coordinate directions using a second-order-accuracy, finite-difference method. At each time moment, \(t_{k}\), the mesh is redefined and adapted to the shape of the water surface, \(z_{\rm s}(x,t_{k})\), according to Equation (9), and all fields computed during the preceding time step are re-defined on the new mesh by a linear interpolation. The integration of Equation (1) is advanced in time by the second-order-accuracy Adams-Bashforth method in two stages to calculate the water velocity at each new time step, \(U_{i}(t_{k+1})\).