Simulation of Fluid Particle Cutting: Validation and Case Study

In this paper we present the comparison of experiments and numerical simulations for bubble cutting by a wire. The air bubble is surrounded by water. In the experimental setup an air bubble is injected on the bottom of a water column. When the bubble rises and contacts the wire, it is separated into two daughter bubbles. The flow is modeled by the incompressible Navier–Stokes equations. A meshfree method is used to simulate the bubble cutting. We have observed that the experimental and numerical results are in very good agreement. Moreover, we have further presented simulation results for liquid with higher viscosity. In this case the numerical results are close to previously published results.


Introduction
Fluid particle cutting plays an important role in gas-liquid and liquid-gas contactors. In gas-liquid contactors, the bubble size distribution, determining the mass transfer area, is B S. Tiwari influenced by the local hydrodynamics, but also by measuring probes such as needle probes [1][2][3] and mesh based conductivity probes [4]. The shape of the probes are mainly cylindrical, while the probe may be in flow direction but also in a rectangular angle to it. The rising bubbles approach the immersed object and starts to change its shape. Depending on the position of the bubble to the wire, the bubble will pass the object or be cutted in two fragments (daughter bubbles). Beside the unwanted cutting at probes, a wire mesh can be used to generate smaller bubbles and homogenise the flow structure. Furthermore, in liquid-gas contactors, phase separation is often a problem. Demisters are then frequently used to prevent a phase slip (entrainment) of fine dispersed phase droplets in the continuous product phase. A loss of the total solvent inventory within one year is reported causing costs and environmental hazards. Entrainment can cause a significant reduction in separation efficiency. Demisters are based on wire meshes, where the small droplets should accumulate. Using an optimal design, the small droplet separation efficiencies can be up to 99.9%. Nevertheless, bigger droplets tend to break up in the rows of wires.
Hence, particle cutting is a frequently observed phenomena in various separation processes ranging from low viscosity to high viscosity of the continuous fluid. Nevertheless, it can be hardly investigated under operation conditions due to the complex insertion of optical probes into the apparatus or the complex mesh structure e.g. of the demister, but also the operation conditions as high pressure, high dispersed phase hold ups make an experimental investigation challenging.
In this study, we focus on the simulation of particle cutting at a single wire strengthened by experimental investigations to generate the basis for further numerical studies at complex geometries and fluid flow conditions such as demister simulations. For the simulation of bubble cutting, a meshfree approach is applied. It overcomes several drawbacks of classical computational fluid dynamics (CFD) methods such as finite element method (FEM) , finite volume method (FVM). The main drawback of the classical methods (FEM, FVM) is the relatively expensive geometrical mesh/grid required to carry out the numerical computations. The computational costs to generate and maintain the grid becomes particularly high for complex geometries and when the grid moves in time, as in the case of fluid particles with a dynamic interface or in case where the interface between fluids changes in time.
For such problems meshfree methods are appropriate. Here, we use a meshfree method, based on the generalized finite difference method, called finite pointset method (FPM). The two phase flow is modeled by using the continuous surface force (CSF) model [8]. Each phase is indicated by the color of the respective particles. When particles move, they carry all the information about the flow with them such as their color, density, velocity, etc. The colors, densities and viscosity values of all particles remain constant during the time evolution. The fluid-fluid interface is easily determined with the help of the color function [9]. In [12] an implementation of the CSF model within the FPM was presented to simulate surface-tension driven flows. We have further extended the method to simulate wetting phenomena [11].

Experimental Setup
Bubble cutting is investigated in a Plexiglas column filled with reversed osmosis water up to a level of 10 cm. The column has a width and depth of 46 mm. A syringe pump (PSD/3, Hamilton) is used to inject air of known volume at the bottom of the column. The injection diameter is 8 mm. The schematic setup is given in Fig. 1. In a distance of 60 mm from the bottom, a wire is mounted in the middle of the column. The wire has a diameter of 3 mm. Two cameras (Imaging Solutions NX8-S2 and Os 8-S2) are mounted in an angle of 90 • to track the bubble motion over an image sequence, respectively, over time. Both cameras are triggered and allow a synchronous detection at 4000 fps and a resolution of 1600 × 1200 px 2 . By tracking the bubble motion from two sides, it is possible to analyse the side movement and to detect the exact position of bubble contact with the wire. Also, the bubble deformation can be analyzed in two direction and therefore leads to more precise results compared to single camera setups.

Bubble Motion Analyses
For bubble motion analyses, the tool box ImageJ (https://imagej.nih.gov/ij/) is used. The raw images (Fig. 2a, b) are therefore binarized, followed by a watershed segmentation. The tracks of the single bubble and cutted particles are analysed using the Plugin Mtrack2 (http://imagej. net/MTrack2). Two particle tracks are tracked, one by each camera and reconstructed using Matlab software toolbox (Fig. 2c) These are the basis for three dimensional reconstruction of the bubble motion. The conversion of pixels to metric length is done by a afore performed calibration.
Matlab is used to reconstruct the bubble in a three dimensional domain. The images are converted to greyscale and further to binary images. Possible holes (white spots in a surrounded black bubble structure) are filled to get a better identification of the bubbles. To detect the bubble position, a distance transform is performed, followed by a watershed segmentation to separate the bubble from the pipe structure. Finally, the bubble size and shape is determined from each image. In a next step, the basic grey scale images are again converted to binary images, followed by a watershed algorithm [5]. Finally, the resulting structures are transformed into 3D space. Therefore, the detected structures are extruded into the third dimension resulting in overlapping structures. The overlapping structures represent the bubble and the wire and are visualized in Fig. 2d. Applying the assumption of an ellipsoidic structure for the bubble, results finally in Fig. 2e.

Mathematical Model
We consider two immiscible fluids, gas and liquid. The equation of motion of these immiscible fluids are modeled by a one-fluid model with the incompressible Navier-Stokes (a) 1. camera raw image.
(b) 2. camera raw image. (e) Visualization based on the assumption of ellipsoidic structure. equations. The incompressible Navier-Stokes equations in the Lagrangian form is given by where x ∈ Ω ⊂ R 2 is the position, v ∈ R 2 is the fluid velocity, ρ is the fluid density, p is the pressure, τ is the stress tensor given by τ = 1 2 (∇v + (∇v) T ), g is the external force and F S is the surface tension force which is the force density and acts in the surrounding of the interface between liquid and gas. The surface tension force is computed using the CSF model of Brackbill et al ( [8]), given by Here, σ is the surface tension coefficient, n is the unit normal at the interface and its surrounding, κ is the curvature and δ S is the surface delta function, which is strong in the interface and its vicinity. The detail computation of these quantities, we will describe in the following section. The Eq. (3) are solved together with initial and boundary conditions.

Numerical Methods
We solve the Eq. (3) by a meshfree Lagrangian particle method. In this method, we first approximate a computational domain by discrete grid points. The grids points are divided into two parts as interior and boundary particles. The boundary particles approximate boundaries and we prescribe boundary conditions on them. The interior particles move with the fluid velocities. Particles may come very close to each other or can go far away from each other leading to very fine or very coarse approximations. This problem has to be tackled carefully due to stability reasons. To obtain a uniform distribution of particles in each time step one has to add or remove particles, if necessary. We refer to [10] for details of such a particle management.

Computation of the Quantities in Surface Tension Force
For meshfree particle methods the interfaces between fluids are easily tracked by using flags on the particles. Initially, we assign different flags or color function c of particles representing the corresponding fluids. We define the color function c = 1 for fluid type 1 and c = 2 for fluid type 2. On the interface and its vicinity, the Shepard interpolation is applied for smoothing of the color functions usingc where m is the number of neighbor of arbitrary particle having position x, c i are the color values at neighboring particle i and w i is the weight as a function of distance from x to x i given by where α is a positive constant After smoothing the color function, we compute the unit normal vector n = ∇c |∇c| .
Finally, we compute the curvature by The quantity δ s is approximated as δ s ≈ |∇c|.
Here, δ s is non-zero in the vicinity of the interface and vanishes far from it.

Numerical Scheme
We consider Chorin's projection method [7] in the framework of a meshfree particle method. Let dt be a time step and set t n = n dt, n = 0, 1, 2, . . .. We denote, for example, x n as the position of a particle at time level n. Chorin's projection scheme consists of two steps, where in the first step we compute the intermediate velocity v * implicitly from the momentum equation without pressure term Due to the Lagrangian formulation we do not have to deal with the nonlinear convective term.
In the second step, called the projection step, we compute the velocity at time level (n + 1) by solving the equation with the constraint that v n+1 satisfies the continuity equation In order to compute v n+1 we need the knowledge of p n+1 . This is obtained by taking the divergence of equation (11) and making use of the constraint (12). Then we get the Poisson equation for the pressure The boundary condition for p is obtained by projecting equation (11) on the outward unit normal vector n at the boundary Γ . Thus, we obtain the Neumann boundary condition where v Γ is the value of v on Γ . Assuming v · n = 0 on Γ , we obtain In addition, we compute the new particle positions at the (n + 1)th level by We note that we have to approximate the spatial derivatives at each particle position as well as solve the second order elliptic problems for the velocities and the pressure. The spatial derivatives at each particle position are approximated from its neighboring clouds of particles based on the weighted least squares method. The weight is a function of a distance of a particle position to its neighbors. We observe that in Eq. 10 there is a discontinues coefficient μ inside the divergence operator since the viscosities of two liquid may have the ratio of up to 1 to 100. Similarly, the density ratio also has 1 to 1000, which can be seen also in Eq. 13. This discontinuous coefficients have to be smoothed for stable computation. This is done using a similar procedure as for smoothing the color function. We denote the smoothed viscosity and density byμ andρ, respectively. We note that we smooth the density and viscosity while solving Eqs. 10 and 13, but keep them constant on each phase of particles during the entire computational time. If the density and viscosity has larger ratios, we may have to iterate the smoothing 2 or 3 times. Finally, Eqs. 10 and 13 can be re-expressed as Note that, for constant density, the first term of Eq. 19 vanishes and we get the pressure Poisson equation. Far from the interface we haveμ = μ andρ = ρ. The momentum and pressure equations have the following general form where A, B and C are known quantities. This equation is solved with Dirichlet or Neumann boundary conditions Remark: For the x component of the momentum equations we have A = 1, B = − dt ρ ∇μ, C = − dt ρμ and f is equal to the right hand side of Eq. 17. Similarly, for the pressure equation Eq. 19 we have A = 0, B = ∇ρ ρ , C = 1 and f =ρ ∇·v * dt . In the following section we describe the method of solving equations Eqs. 20-21 by a meshfree particle method, called the finite pointset method (FPM).

A meshfree particle method for general elliptic boundary value problems
In this subsection we describe a meshfree method for solving second order elliptic boundary value problems of type Eqs. 20 and 21. The method will be described in a two-dimensional space. The extension of the method to three-dimensional space is straightforward. We consider the computational domain Ω ∈ R 2 . Approximate Ω by particles of positions x i , i = 1, . . . , N , whose distribution can be irregular. These particles serve as numerical grid points. Let ψ(x) be a scalar function and ψ i = ψ(x i ) its values for i = 1, . . . , N . We consider the problem to solve the Eqs. 20 and 21 at an arbitrary point x ∈ {x i , i = 1, . . . , N }, in terms of the values of a set of its neighboring points. In order to restrict the number of neighboring points we define a weight function w = w(x i −x, h) with small compact support of size h. The value of h has to be chosen such that we have at least a minimum number of particles, for example, in 2D, we need at least six neighboring particles. In practice we define h as 2.5 to 3 times the initial spacing of particles, keeping in mind that this is a user defined factor. We consider the Gaussian weight function defined in (6), where α is equal to 6.25. In general, the value of α has to be chosen according to the choice of h such that the approximation of spatial derivatives is accurate. In this paper, we have chosen h equal to 3 times the initial spacing of particles, so this choice of α gives an accurate approximation of spatial derivatives. Let P(x, h) = {x j : j = 1, 2, . . . , m} be the set of m neighboring points of x in a circle of radius h. The point x is one of x j .
Consider m Taylor expansions of ψ( for j = 1, . . . , m, where e j is the residual error. We denote the coefficients We add the constraint that at particle position (x, y) the partial differential equation (20) should be satisfied. If the point (x, y) lies on the boundary, also the boundary condition (21) needs to be satisfied. Therefore, we add Eqs. 20 and 21 for the Neumann boundary condition to the m equations (22). Equations 20 and 21 are re-expressed as where B = (B 1 , B 2 ) and n x , n y are the x, y components of the unit normal vector n on the boundary Γ . The coefficients a i , i = 1, . . . , 6 are the unknowns. We have six unknowns and m + 1 equations for the interior points and m + 2 unknowns for the Neumann boundary points. This means, we always need a minimum of six neighbors. In general, we have more than six neighbors, so the system is overdetermined and can be written in matrix form as with the vectors given by a = (a 1 , a 2 , . . . a 6 ) T , b = (ψ 1 , . . . , ψ m , f , ψ N ) T and e = (e 1 , . . . , e m , e m+1 , e m+2 ) T and dx j = x j − x, dy j = y j − y. For the numerical implementation, we set n x = n y = 0 and ψ Γ N = 0 for the interior particles. For the Dirichlet boundary particles, we directly prescribe the boundary conditions, and for the Neumann boundary particles the matrix coefficients are given by Eq. 26. The unknowns a i are computed by minimizing a weighted error over the neighboring points. Thus, we have to minimize the following quadratic form The minimization of J with respect to a formally yields ( if M T W M is nonsingular) In Eq. 28 the vector (M T W )b is explicitly given by Equating the first components on both sides of Eq. 28, we get where Q 1 , Q 2 , . . . , Q 6 are the components of the first row of the matrix (M T W M) −1 .
Rearranging the terms, we have We obtain the following In matrix form we have where R is the right-hand side vector, is the unknown vector and L is the sparse matrix having non-zero entries only for neighboring particles. The sparse system (33) is solved by an iterative method. In this paper we apply the method of Gauss-Seidel. In the time iteration the initial values of ψ for time step n +1 are taken as the values from time step n. We note that, for solving the equations for intermediate velocities and the pressure Poisson equation will require more iterations in the first few time steps. After a certain number of time steps, the values of pressure and velocities at the old time step are close to those of new time step, so the number of iterations required gets reduced.
The iteration process is stopped if the relative error satisfies

Validation: Low Viscosity
In a first step, we validate the simulations with the experimental results of the single bubble cutting in reversed osmosis purified water. We apply the same bubble diameter from the experiments (6.5 mm) and the wire diameter of 3 mm in the simulation. The viscosity of the fluid (water) is μ l = 0.001Pa s and the interfacial tension between water and air is σ = 0.072 N/m. The density of water is ρ l = 998.2kg/m 3 and the density of air is approximated by ρ g = 1kg/m 3 and the dynamic viscosity of air is μ g = 2e −5 . For the numerical simulations we consider a two-dimensional geometry of size 36mm × 63mm. The initial bubble position has the center at x = 18 mm and y = 10 mm and the wire has the center at x = 19.5 mm and y = 45 mm as shown in Fig. 3. The bubble is approximated by red particles, the liquid is approximated by blue particles. The white circular spot is the position of the wire. We have considered the total number of boundary particles, (including 4 walls and the wire) equal to 527 and the initial number of interior particles equal to 18, 310. The constant time step t = 5e −6 is considered. Here the horizontal distance between the initial center of bubble and the wire is d x = 1.5 mm. In all four walls and the wire we have considered no-slip boundary conditions. Initially, the velocity and pressure are equal to zero. The gravitational force is g = (0, −9.81)m/s 2 .
The comparison between simulation and experiment is depicted in Figs. 4 and 5. We extracted a time sequence from the experiments and the corresponding simulations, starting at 0.19 s simulated time to 0.27 s. The temporal distance between each image is 0.02 s. The rising bubble approaches the wire and starts to deform. There is no direct contact during this phase between the bubble and the wire. Due to the non central approach to the wire, the bubble is cut in a smaller daughter bubble (right) and a larger bubble (left). The larger bubble has three times the diameter of the smaller bubble. The comparison of the cutting process gives a qualitatively good agreement between the experiment and the simulation. Also the shape and size of the mother and daughter bubbles have qualitatively very good agreement. A detailed comparison of the bubble path from experiment and simulation is shown in Fig. 6. The bubble position during first contact is important, which agrees well between experiment and simulation. Nevertheless, the path of the larger bubble in the simulation shows after the cutting a slightly different behaviour than in the experiment. In the experiment, the larger bubble moves inwards again, while the bubble in the simulation moves horizontally away from the wire, which may arise from a slight horizontally movement of the bubble in the experiment. The cutting of the bubble also depends on the bubble velocity, plotted in Fig. 7. The bubble accelerates in the simulation and finally reaches the same end velocity that is observed in the experiment. After the splitting into two daughter bubbles, the larger daughter bubble raises faster than the smaller one. The experimental results are governed by higher fluctuations especially for the smaller bubble, which results from very short temporal distances between the images and fluctuations by detecting the bubble interface. Neverthless, the average velocity for the smaller bubble fits with 0.12 m/s quit well to the simulated result.

Case study: high viscosity
The computational domain is the same as in the previous section. However, the position of the wire is changed. The data has been taken from chapter 6 of [6]. The liquid has density ρ l = 1250 kg/m 3 , dynamical viscosity μ l = 0.219 Pa s. Similarly the gas density ρ g = 1kg/m 3 and the viscosity μ g = 2e −5 Pa s. The surface tension coefficient σ = 0.0658N/m. We consider a bubble of diameter 9.14 mm with its initial center at (18mm, 9mm). We consider 14 The rising velocities of the mother and the daughter bubbles a wire (in 2D a circle) of diameter 3.1 mm with different centers at y = 45 mm and x = 18mm, 18.5mm, 19mm and 19.5mm. This means, we consider the initial distance d x between the center of the bubble and the center of the wire equal to d x = 0mm, 0.5mm, 1mm, 1.5mm. Figure 8 shows the initial geometry with d x = 0mm. The initial number of particles and the time step are the same as in the low viscosity case. The initial and boundary conditions and the rest of other parameters are same as in the previous case.
In Figs. 9, 10, 11 and 12 we have plotted the positions of the bubble and the wire for d x = 0 mm, 0.5 mm, 1 mm and 1.5 mm at time t = 0.288, 0.320, 0.352 and t = 0.384 s, respectively. For d x = 0 we observe the wire located in the middle of the bubble as expected. When we increased the distance d x from 0.5 mm to 1.5 mm, we observed that the left part of the bubble is increasing and the right part becomes smaller. We clearly observe that the daughter bubbles are symmetric for d x = 0 mm in contrast to the other cases. We further observe a small layer between the wire and the bubble. After t = 0.352 s we observe the cutting of the bubble, see Figs. 11 and 12. Two daughter bubbles arise, a larger one on the left and a smaller one on the right side of the wire. The overall numerical results are comparable with the results presented in [6].
In Fig. 13 we have plotted the trajectories of the mother and bubble droplets. We observed that the mother droplet is cutted into two daughter bubbles slightly below the wire, compare with Fig. 11. The trajectories are plotted up to time t = 0.5 s. We see that when the size of the daughter bubble is increasing, it travels longer than the smaller bubbles. The reason is that the rising velocity of the larger bubble is larger than the smaller ones, see Fig. 14 for the velocities of mother and daughter bubbles.

Concluding Remarks
The cutting of bubbles at a single tube (wire) was investigated experimentally and numerically. For the simulations, a meshfree method was applied. The method enables a description of the deforming interface and the hydrodynamics of bubble cutting. For a first validation, we compared the solver to experimental data using the system air and water. A sufficiently good agreement could be found in regard to bubble shape, bubble movement and cutting process itself. To study the effect of higher viscosity and bubble position, a case study was done. One observes that the initial position of the bubble to the wire has a high impact on the final daughter bubble size ratio. A centric approach of the bubble to the wire leads to a cutting of the bubble in two equally sized daughter bubbles. By increasing the initial distance to the wire, the daughter-bubble size ratio increases and the deviation between the velocities of daughter bubbles increases. Also, the movement of the bubbles directly behind the wire changes. While the bubbles split behind the wire at the centric approach, with increasing unsymmetry, the bubbles start to move inwards after the initial separation. In future, further studies with overlapping wires are planned.