Iterative Local Bézier Reconstruction Algorithm of Smooth Droplet Surface for the Immersed Boundary Method

Empirical formulae of the terminal velocity and shape of a water droplet in microphysics parametrizations are derived from experiments or theoretical works and they are only verified under room temperature and standard atmospheric pressure. A two-phase direct numerical simulation model can be a strong tool to improve those empirical formulae under general conditions. Recently, the global B-spline fitting has been applied and its smooth one dimensional (1D) surface reconstruction of water droplets has enabled stable simulations of falling two dimensional (2D) droplets by the immersed boundary method (IBM). However, an extension of the global reconstruction from 1D to 2D is highly complex and is almost impossible to use in a model. To overcome this limitation, an iterative algorithm is proposed for a local smooth surface reconstruction in this work. One significant advantage is its straightforward extension to 2D surfaces. To test the new method, simulations of an axisymmetric free-oscillating water droplet are compared between the global and local surface reconstructions. A further simulation of a rising air bubble is performed to examine the robustness of the new algorithm for the highly distorted interface. This new method opens a pathway to three dimensional (3D) water droplet simulations by the IBM. (Citation: C. R., Ong, and H. Miura, 2018: Iterative local Bézier reconstruction algorithm of smooth droplet surface for the immersed boundary method. SOLA, 14, 170−173, doi:10.2151/ sola.2018-030.)


Introduction
The cloud microphysics is composed of processes that control the formation and growth of hydrometeors in a cloud.These processes are usually parametrized by empirical relationships that are characterized by parameters such as the terminal velocity, shape, and oscillatory amplitude of droplets.In many cases, these parameters were derived from theoretical and/or experimental works that assumed room temperature and standard atmospheric pressure (Gunn and Kinzer 1949;Beard and Pruppacher 1969;Beard and Chuang 1987;Feng and Beard 1991), and thus the correctness of those parameters is questionable if they are used under different environmental conditions.Therefore, as a method to verify or improve those parameters, two-phase flow direct numerical simulations (DNS) can be a breakthrough.
The immersed boundary method (IBM) is one of those twophase flow DNS models, which has been applied for bubble or droplet simulations (Peskin 1972;Unverdi and Tryggvason 1992;Peskin 2002;Mittal and Iaccarino 2005).In the conventional IBM, the surface of a droplet is represented by a series of connected Lagrangian markers and its curvature is estimated from local cubic polynomial fitting of the locations of the markers.The estimated curvature is then used to compute surface tension force at each marker location to satisfy the interfacial boundary condition.Since the local cubic polynomial fitting is sensitive to small noises, numerical errors introduced in updating the locations of markers are often magnified and lead to a numerical instability that is characterized by a fluctuating droplet surface.The length scale of the fluctuation is two to three marker spacings.The authors have shown that one can settle this instability problem by using the global B-spline fitting instead of the local cubic polynomial fitting (see Ong and Miura 2018, hereafter OM, or Supplement 1).The markers are redistributed along the reconstructed B-spline surface at every time step so that marker intervals are equidistant.
Although a stable 2D water droplet simulation has been achieved by using the global B-spline fitting (OM), its extension to 2D surface reconstruction in 3D water droplet simulations is not straightforward.The difficulty comes from the global fitting.A global spline reconstruction of a 2D surface either requires higher degree of polynomials (Alfeld 1984;Farin 1986) or high computational cost (Pfeifle and Seidel 1994) to maintain C 2 continuity (continuity of the second derivative on the surface), which is needed to evaluate the surface stress.Moreover, a global smoothing of a 2D surface requires the solution of a large system of linear equations.Adjustments of the locations of markers in equal intervals on a 2D surface add an additional complexity.Thus, although the IBM has been improved largely, its applications are still restricted to 2D or axisymmetric smaller water droplets.This strongly limits its future application to the cloud microphysics study.
To avoid these difficulties, we propose an iterative local Bézier reconstruction algorithm to replace the global B-spline fitting.It is expected that this local smooth reconstruction will be more readily extended for 3D large water droplet simulations.The overview of this paper is as follows.The new local smooth surface reconstruction algorithm is explained in Section 2 and is subjected to several test cases in Section 3. Finally, Section 4 concludes this work.

An iterative local Bézier reconstruction algorithm
We explain each step of the iterative local Bézier reconstruction algorithm below and the schematic view of the algorithm is shown in Fig. 1.
If a set of control points (X 0 , Y 0 ), (X 1 , Y 1 ), …, (X d , Y d ) (black open squares in Fig. 1) is given, a Bézier curve in the 2D Cartesian coordinates is defined by where B x d (u) and B y d (u) are the x and y coordinates of the Bézier curve, d i ( ) is the binomial coefficient d!/[i!(di)!], X i and Y i denote the x and y coordinates of control points, u is the curve parameter, and d is the degree of polynomial.We set the number of the control points to five, this is because the five control points correspond to the five marker points of the cubic polynomial fitting.The degree of polynomial of the Bézier curve is therefore four (d = 4).A Béizer curve always passes through the first and last control points, while the intermediate control points do not necessarily lie on the curve.The parameter u is defined over the range of [0, 1].This range is mapped into the piecewise linear interface [-l, l ] (blue lines in Fig. 1a) and this mapping is one-toone.The center of the range [0,1] corresponds to the parametric Ong and Miura, Iterative Local Bézier Reconstruction Algorithm ment of each marker is only performed after the least squares fitting has been finished for all markers.These Bézier curve fitting and location adjustments are repeated several times until the surface becomes sufficiently smooth, in other words, surface fluctuation of length scale smaller than two to three marker spacings is smoothed out.In our tests, five iterations or less are enough to obtain a smooth 1D surface when l = 5∆ s.However, it is possible that a larger number of iterations are required when this algorithm is extended to 3D simulations.We will examine this point in a future work.The number of iterations also depends on the range of least squares fitting.Smaller l tends to require more iterations.It is noted that this iteration step works similarly to the low-pass filter of Gaussian smoothing based on the discrete Laplace-Beltrami operator (Taubin 1995) to filter out high frequency noises.Finally, the curvature κ is computed on the smooth Bézier curve by the following equation: where the prime denotes the derivative with respect to u.An obvious advantage of this newly proposed local reconstruction algorithm over the global B-spline fitting is that it can be easily extended to an arbitrary 2D interface.One just needs to construct a tensor-product Bézier surface at each marker location and the rest of the algorithm remains exactly the same.Unlike the B-spline curves or surfaces, reconstruction by Bézier patches has a disadvantage that it is difficult to impose higher than C 1 continuity between neighboring patches; where C 1 means continuity of the first derivative.However, in the present algorithm, a smooth connection between neighboring patches is not necessary.A pair of patches are allowed to overlap onto or be apart from one another and the iteration of least squares fitting provides smooth transition from one patch to another patch.Thus, this difficulty is not crucial.

Results
Simulations of a free-oscillating droplet and a rising bubble were performed to test the accuracy and the robustness of the new local smooth surface reconstruction algorithm.The fundamental equations of the two-phase incompressible flow problem and the numerical schemes employed in the flow solver are explained in Supplement 2.
In free-oscillating droplet simulations, an axisymmetric water droplet of elliptical shape is suspended freely in an ambient fluid and undergoes oscillatory motion without gravitational pull.We set the equivalent radius of the droplet r = 0.500 cm.Simulations of three different grid resolutions were performed and the diameter of the droplet is non-dimensionalized to one (see Supplement 3 for the numerical settings).The Cartesian grid spacings ∆ x are 0.05, 0.025, and 0.0125, respectively.The marker spacing is equal to the grid spacing, i.e. ∆ s = ∆ x.The distance between two adjacent markers is kept between 1.5∆ s and 0.5∆ s (Singh and Shyy 2007; OM) and the initial number of markers on the droplet surface for these three grid spacings are 33, 64, and 127 respectively.The number of iterations for the new algorithm is fixed to five.The time step is 4.0 ´ 10 -3 for all grid resolutions so that the same number of smoothing iterations is applied to the droplet surface in a fixed period of time.The simulations were run until time t = 800.Note that space and time are also dimensionless in the actual simulations.The unit non-dimensional time corresponds to 1.095 ´ 10 -3 s in its dimensional form (see Supplement 3 for the conversion).
Based on the capillary wave theory (Lamb 1932;Rayleigh 1879), an analytic solution exists for the periodic motion of the droplet when the perturbation is small and its period T is given by location of the Lagrangian marker under consideration (red solid circle in Fig. 1), i.e. the zero point in [-l, l ].There is no theoretical restriction on the determination of l, but larger l causes interfacial motions of larger spatial scale to be damped selectively and the boundary condition of the continuity of velocity across a fluid-fluid interface can break down.Thus, l should be determined accordingly so that only oscillations of spatial scale below two or three marker spacings is removed.Since the fluctuation of length scale below about 0.5 l is expected to be smoothed by the least squares fitting of the Bézier curve, we choose l = 5∆ s, where ∆ s is the marker spacing, in order to eliminate the surface fluctuation of length scale below about 2.5∆ s.To demonstrate that the erroneous surface fluctuation happens when l is not sufficiently large, tests with l = 3∆ s were also performed and compared with l = 5∆ s.The range of the parameter u is fixed to [0, 1] regardless of the values of l because a Bézier curve is invariant under an affine transformation.
We perform the least squares fitting by minimizing δ x and δ y that are defined by the following equations: where the summation about i is taken over all markers lying within the range [-l, l ], u i is the marker curve parameter, and x i and y i are the marker coordinates.This minimization problem is solved by taking derivative of the equations above with respect to the control points X j and Y j and forcing them to be zero, , , .for k = … 0 4 (3) These are linear equations, and can be solved by any linear matrix solver to obtain the control points X j and Y j .In this work, the LU decomposition is used.With these X j and Y j ( j = 0, …, 4), a Bézier curve is determined by eq. ( 1).
The location of the marker under consideration is adjusted to the center of a Bézier curve computed from eq. ( 1).This adjust- where ρ d is the droplet density, ρ a is the ambient fluid density, and n is the mode of oscillation.By substituting n = 2, ρ d = 998 kg m -3 , ρ a = 1.20 kg m -3 , and σ = 0.0728 N m -1 into the above equation, the analytic period is 0.0920 s.Errors in the period of the first oscillation of the numerical solutions are evaluated to compare the convergence rate and accuracy between the new local surface reconstruction algorithm and the global B-spline fitting.
The periods of the first oscillation cycle are summarized in Table 1.The iterative local Bézier fitting achieves as similar accuracy as the global B-spline fitting for grid spacings smaller than 0.025.The numerical solution converges to the analytic solution as the grid size decreases.However, the accuracy of the new algorithm becomes worse when ∆ x = 0.05 and l = 5∆ s.This indicates that the oscillatory motion is quickly damped by the iterative least squares fitting.Furthermore, due to this excessive smoothing, the interfacial boundary condition of continuous velocity field is broken down because the interfacial velocity deviates largely from the interpolated fluid velocity at each marker location.This artificial damping does not seriously affect the oscillatory motion when the grid spacing is finer than 0.025.
Figure 2 shows the time evolution of the minor axis diameter, which is defined by the length of the line segment cutting through the center of the droplet in y-direction.The shape of the droplet when its minor axis diameter approximately attains its maximum and minimum in the first oscillation cycle are also shown.The amplitude decreases with time due to the viscous dissipation and the shape of the droplet slowly restores back to a sphere.Snapshots of curvature at the location of each marker using the global B-spline fitting and the iterative local Bézier reconstruction algorithm are compared in Fig. 3 for the grid size of 0.0125.The curvature remains smooth even though the distance between neighboring markers are not equidistant in the new algorithm when l = 5∆ s.However, the simulation of l = 3∆ s suffers from a small erroneous surface fluctuation.In this case, l = 3∆ s is too small to operate a sufficient smoothing.This explains why the convergence rate is slightly worse when l = 3∆ s in Table 1.
In the above test, the deformation of the droplet was limited to a small amplitude.To test the robustness of the new algorithm, we have performed a free-oscillation simulation with a large amplitude (see Supplement 4).It was found that the new algorithm successfully maintained a smooth interface for a much larger deformation.
Hereafter, we present simulation results of a rising air bubble in water with weak surface tension coefficient.This experiment is designed to demonstrate that a complex surface can be treated with modest accuracy by the new algorithm.The shape of the bubble in the rising bubble problem can be categorized into several regimes such as the spherical, ellipsoidal, and skirted shape regimes when the terminal velocity is reached (Clift et al. 1978, hereafter CGW).Furthermore, the bubble shape and its terminal velocity are characterized by four dimensionless numbers, namely the density ratio, viscosity ratio, Eötvos number, and Morton number (see figure 2.5 in CGW).The Eötvos number E and the Morton number M are defined as Table 1.Comparison of the period of first oscillation cycle with different fitting methods.The analytic period is 0.0920 s which is calculated from eq. ( 5) with the oscillatory mode n equal to 2. The conversion of the dimen sionless period obtained from the numerical solution to its dimensional form is given in Supplement 3. The number of the significant digits of the period is three.

E M (7)
This is the simulation of an air bubble of diameter D = 4.77 ´ 10 -1 mm rising in water at 20°C with surface tension coefficient σ = 2.14 ´ 10 -5 N m -1 .These numbers are chosen because, based on the experimental data (figure 2.5 in CGW), the bubble shape belongs to the complicated skirted shape regime, and thus matches our purpose to test the new algorithm in simulations of complex water-air interface.The simulation was run from t = 0 to t = 3.0 when steady state is approximately attained and the unit nondimensional time corresponds to the dimensional time of 2.19 ´ 10 -2 s (see Supplement 5 for the details of the computational domain and numerical settings).
The time evolution of the shape until the steady state is reached is shown in Fig. 4. The shape changes gradually from a spherical shape to an inverted "U" shape and further to a skirted shape.The final shape of the bubble agrees qualitatively well with the experimental data in CGW.Moreover, the rising speed and Reynolds number at steady state are 47.1 mm s -1 and 22.4, which are also in agreement with the experimental data.There is no erroneous fluctuation either in the flattened area on the bottom or near the strongly curved area along the sides of the bubble.This test demonstrates that the IBM with the local Bézier smoothing algorithm is robust for both large and small deformations of the fluid-fluid interface with a large density ratio.

Conclusion
A new iterative local Bézier reconstruction algorithm has been proposed to enable a smooth surface reconstruction of a fluid-fluid interface in the IBM.By using this new method, we can avoid difficulties in extending the global B-spline fitting to a 2D interface.Numerical tests have shown that there is little difference in the convergence rate and accuracy between the global B-spline fitting and the present local surface reconstruction algorithm.The simulation of a rising bubble was performed to demonstrate that the new algorithm is robust even when the surface is highly distorted.Based on these results, we conclude that the new algorithm can replace the global B-spline fitting in stable water-droplet simulations.However, the collision and coalescence process, which is an important process in cloud formation, is still not enabled by the new algorithm.The criteria for the onset of coalescence of two separate surfaces and the accurate redistribution of markers when they start to merge together are among the remaining tasks.These are left for future study to enhance applicability of the IBM.

Fig. 1 .
Fig. 1.Schematic view of the new reconstruction and adjustment algorithm.(a): A fixed range (blue line [-l, l ]) from a particular marker (red solid circle) along the polygonal interface is first determined.(b): Bézier curve is formed by least squares fitting of the markers lying within the range, and then the red marker is adjusted to the new location in the middle of the curve as shown by the blue arrow.Black open squares denote the control points of the Bézier curve.This is an example of d = 4.
Fig. 3. Plot of curvature along the interface from the first to the last marker with respect to their angles at time t = 4.0 (blue), t = 20.0 (green), and t = 40.0(red).(a): Global B-spline fitting by least squares.(b): Iterative local Bézier surface reconstruction algorithm (l = 5∆ s).(c): Iterative local Bézier surface reconstruction algorithm (l = 3∆ s).The first marker is located at the uppermost point of the droplet at zero degree and the last marker is at 180 degrees.The number of markers is 64 and the grid spacing is 0.0125.