Modeling and Simulation of Coalescence in the Context of Oil Mist Filtration

Today, aerosols are generated in many industrial processes like machining or in pneumatic compressors. The harmfulness of these droplets often lies in the liquid itself, as it can be toxic or cancerous while exhaled natural aerosols can carry pathogens. For health protection, the filtration of these aerosols is an important task. The wetting properties of the fibers have a large influence on the filtration efficiency and coalescence of droplets on the fibers. Oleophilic fibers are better in capturing the oil particles, while drainage increases using oleophobic fibers. For improved efficiency and drainage properties, a mixed-wet fiber system could therefore be advantageous. Describing droplet capture and coalescence in filter-scale simulations is computationally extremely demanding. We present a new approach for simulating oil mist deposition and droplet coalescence, specifically accounting for the wettability of the fibers, especially for mixed-wet filter media. The aim is to gain a more realistic pressure loss of a fibrous filter structure which takes the presence of the filtered phase into account as well as the coalescence of the liquid phase on the fibers depending on the wettability and the fluid parameters. Our approach is to model the shape of deposited droplets via weighted distance maps, according to the contact angles. Starting from a sphere, the shape is fitted on the solid surface such that the contact angle is met. For validation purposes, volume-of-fluid-based simulations were set up using the same droplet-fiber system as a reference to compare the droplet shape and pressure loss in the presence of oil droplets on the fiber surface.


Introduction
Liquid aerosol particles are generated in many industrial processes like machining or in pneumatic compressors. The harmfulness of these industrially produced droplets often lies in the liquid itself as it can be toxic or cancerous, while exhaled natural aerosols can carry viruses or other pathogens.
For health protection, the filtration of these aerosols is an important task. The investigation of solid aerosol particles has been in the focus of research for a long time, and the deposition mechanisms are well known, Liew et al. (1985), Jackson and James (1986), Frising et al. (2004), Mullins et al. (2006Mullins et al. ( , 2007Mullins et al. ( , 2009, Jaganathan et al. (2008). This is different for liquid aerosol particles which have only been investigated in depth since about 20 years, King et al. (2010), Piller et al. (2014), Hall et al. (2016). This is due to the complex nature of droplet separation. Gac and Gradoń (2011) simulated a 2D binary coalescence of droplets with an lattice-Boltzmann method to specify the influence of the droplet size and the interfacial tension on binary coalescence time. Frising et al. (2005) presented a phenomenological model describing the pressure drop and penetration evolution for liquid aerosol deposition. Kampa et al. (2014) introduced a phenomenological model to explain the increase in pressure drop of air filters during steady operation with oil mist, based on semi-quantitative conclusions obtained from measurements of liquid distribution patterns. The wetting properties of the fibers have a large influence on the filtration efficiency and coalescence of droplets on the fibers. Oleophilic fibers are better in capturing the oil aerosol particles, while drainage increases using oleophobic fibers. Investigating the deposition of liquid aerosols, many experimental works were done. Often, the increase of pressure loss over the filter lifetime and the evolution of fractional filtration efficiency and its dependence on parameters like flow velocity and fluid properties were investigated, see, e.g., Payet et al. (1992), Raynor and Leith (2000), Charvet et al. (2008Charvet et al. ( , 2010, Shin et al. (2005), Lehmann et al. (2005), Frising et al. (2005). Mead-Hunter et al. (2013) have shown that CFD simulations can provide a viable alternative to experimental testing of filter materials, by simulating flow through a section of a fibrous porous medium on the micro-scale, especially when high-performance computational resources are available. They modeled both the deposition of droplets on the fibers using a Lagrangian approach as well as the coalescence of the liquid using a volume-of-fluid (VOF) method in the software OpenFOAM. Abishek et al. (2019) evaluated the influence of equilibrium contact angle and contact angle hysteresis, on the dynamics of a microdroplet colliding with an isolated fiber during mist filtration via VOF simulations in OpenFOAM.

3
Simulating wetting behavior of fluids on complex structures bears a problem in many different application areas. One example is the simulation of the Mercury Intrusion Capillary Pressure (MICP) experiment, which is used to determine information about a pore space morphology, e.g., in rock samples. In this experiment, mercurywhich has most un-wetting behavior-is pushed into the pore space, while pressure and saturation of the available volume are measured. Several methods have shown, see Hilpert and Miller (2001), Schulz et al. (2015), that it is possible to simulate this process on a geometrical basis. Simulating on a regular 3-dimensional voxel grid makes it possible to use efficient distance map algorithms, i.e., Saito and Toriwaki (1994). The distance map contains the distance to the nearest solid voxel from each voxel inside the grid. This distance will serve as a basis for many further calculations, like the geometric pore size distribution, the maximum percolation path diameter or even the fluid phase wetting behavior. Using the distance map as a basis is much faster and easier than solving a complex two-phase flow problem. Starting from a sphere, the shape is fitted on the solid surface such that the contact angle is met using the WDM approach.
Although there have been efforts to model liquid aerosol particle separation using CFD, there is still no simulation approach that models aerosol particle separation, coalescence and the resulting pressure loss for realistic fibrous filter structures.
Therefore, the purpose of this work is to 1. Develop a weighted distance map (WDM) approach to simulate the coalescence and wetting of liquid particles in a filter medium 2. To perform a two-step validation by comparing simulation results using the WDM approach to results obtained using the volume-of-fluid (VOF) method 3. To compare fluid distribution and pressure drop as a function of time for a complex fiber structure using the WDM approach and the volume-of-fluid model.
This work is structured as follows: In Sect. 2, the basic equations are discussed. Section 3 describes the WDM method and how it is implemented. Then, in Sect. 4, the validation of the VOF model is presented. The two-step validation of the WDM method with the VOF model is described in Sect. 5, while in Sect. 6, the results are discussed and Sect. 7 provides a conclusion.

Theory and Basic Equations
Although this work deals with the filtration of aerosol particles, it does not consider the aerosol particle transport in the gas phase and the particle collection on the fibers. Instead, our considerations start once the aerosol particles have attached to the fibers. For the aerosol transport and separation mechanisms along with the modeling of aerosol particle deposition, we refer to Hoch et al. (2020). In the following, we describe the numerical and analytical methods used for the simulation and validation of the VOF model as well as the distance map-based method. For model development and since fluid properties are well known, we use water as the liquid phase in the following.

Volume-of-Fluid (VOF) Method
The VOF method is a common free-surface modeling technique and belongs to the Eulerian methods. This numerical technique allows for tracking the fluid-fluid interface of Eulerian phases like the surface of a droplet. Since the VOF method is not a standalone solving algorithm, the conservation equations have to be solved in addition. The phase mass conservation equation of a multi-phase system is given by: where i is the volume fraction and i the density of phasei . Furthermore,v is the mixture velocity, v d,i the diffusion velocity, S i is a user-defined phase mass source term, and D i Dt is the material or Lagrangian derivative of phase density i .
The momentum equation is given by: with p as the pressure, I the unity tensor, T the stress tensor and f b the vector of body forces. The volume fraction is tracked in each cell and sums up to 1 for all phases. All fluids share a single set of momentum equations and one velocity field. Therefore, the pressure loss of the gas phase alone cannot be determined in VOF simulations.
In order to obtain a sharp fluid-fluid interface, a 2nd-order discretization scheme is used without a high-resolution interface capturing scheme. For modeling the surface tension force, a volumetric force using the continuum surface force approach by Brackbill et al. (1992) is used. For stability purposes, the momentum source term due to surface tension is formulated in a semi-implicit way. Thus, a larger time step can be realized in simulations where the surface tension force is the decisive force. In addition, the simulation becomes more stable and the interface sharper.

Surface Evolver
The surface evolver is a text-based simulation tool for modeling liquid surfaces shaped by various forces and constraints. It aims at minimizing the total surface energy and includes the surface tension, gravitational energy, squared mean curvature and user-defined surface integrals. The gravitational potential energy is described as with g as the gravitational acceleration and k the face dimension, while z is the vertical coordinate of the facet. The energy on a facet due to surface tension is governed by the following equation: where is the surface tension and s 0 as well as s 1 are the edges in counterclockwise order around a facet. The surface evolver can be used for steady-state simple geometrical setups including two fluid phases and allows for a fast calculation of fluid-phase surfaces. The well-documented open-source tool by Brakke et al. (2021) is used in this work.

Analytical Solutions
The capillary rise h of a liquid in a cylindrical tube is described as follows: with as the contact angle between liquid and solid and a the radius of the tube.

WDM Approach Describing Droplet Coalescence
The new geometric method for modeling coalescence of droplets is implemented in the software GeoDict® by Math2Market GmbH (2022). Within GeoDict, all calculations are based on a regular 3-dimensional voxel grid. GeoDict has various capabilities of modeling and simulating microstructures of all kinds. Distance map calculation and 3D structure handling are well-tested basic functionality of GeoDict.
Similar to the VOF method, the 3D volume field in GeoDict can represent the liquid volume fraction of each voxel. If a voxel is completely filled with fluid, the value will reach 1.0, while voxel containing no fluid will show value 0.0. Besides this volume fraction field that describes the fluid distribution, GeoDict holds another 3D voxel grid-the voxel geometry-that represents the information about the solid microstructure with which droplets may be interfering.
Coalescence procedure will follow several distinct steps, that will only change the liquid volume fraction field, while the voxel geometry will stay unchanged. Nevertheless, we need the voxel geometry to account for fluid-solid interface interactions, thus respecting the contact angle. Note that it is possible to have several different materials with different contact angles in the voxel geometry.
The first calculation steps concern the merging of droplets that are in contact with each other: As a starting point (see Fig. 1, first image), droplets are assumed to be perfectly spherical and are placed in the volume fraction field. Voxels on the boundary In the next step (see Fig. 1, second image), we identify droplet volume which is in contact with other droplets by labeling connecting components. For this, we assume a threshold of 0.5 as a value to differentiate between liquid and empty voxels (air). Once all connecting components have been labeled, we sum up all volumes of each label and find its center of mass. Next, the center of mass and volume are used to re-distribute fewer, but larger, droplets again in the volume fraction field.
The resulting volume fraction field (see Fig. 1, right image) will now contain coalesced droplets, but still with perfect spherical shape. Thus, a further step is needed to find a correct shape according to contact angle and solid structure surrounding the droplets.
In further steps, the newly created liquid volume fraction field is morphed in such a way that suitable contact angles on solid surfaces are obtained, while the volume of all droplets is conserved. As an example, we take two droplets with different size and different contact angles sitting on a solid surface (see Fig. 2, left image). The volume is subdivided, such that each distinct droplet has its own section (see Fig. 2, right image). With this, we can ensure volume conservation and thus, mass conservation, for each droplet.
Basis for further steps is two distance maps calculated on solid voxels, see Fig. 3. Specifically, these are the wall distance field (Fig. 3 left image) and droplet centers distance field (Fig. 3 right image).
The idea behind is that a threshold of a distance map of the droplet center would result in a perfectly spherical object, while thresholding the wall distance would lead to a thin layer on the whole surface. Correspondingly, thresholding those two distance fields would resemble a perfectly hydrophobic ( = 180 • , spherical) or hydrophilic ( = 0 • , surface-layer) liquid shape. Finally, we are searching for some shape in between these two solutions by combining both distance maps.
In order to find a suitable combination method, we investigated how well the analytical solution of droplets on a flat surface could be resolved by several approaches. Potential equation parameters were found using an automated minimization algorithm. Further work could be done to investigate whether there is a physically motivated approach. Similar problems concerning contact angle and fluid surface shapes have already been solved in the SatuDict module within the GeoDict software (Silin and Patzek 2006;Saxena et al. 2021). A practical approach was found in combining the distance maps by applying an exponent and simply adding them up. A suitable center distance exponent of 1.25 was found, mostly independently of the choice of the wall distance exponent. For the wall distance exponent function E( ) a simple linear fit was quite close to the analytic solution and analysis of the remaining error revealed some missing sinus correlation.
The two distance fields are combined in the next step, using the contact angle as a balancing factor, with where all parameters have been chosen using a least-squares fit. The resulting combined distance (Fig. 4, left image)4 is then thresholded, such that the resulting volumes correspond to the initial droplet volumes (Fig. , right image).
The resulting binary field will have droplets with correct contact angle and approximate correct shape. In a next step, it is also possible to add more droplets in another step and rerun the process.

Validation of the VOF Model
To ensure that the results of the VOF simulation and the distance map method are correct, they must be validated. The validation is split into two main steps. Concerning the volumeof-fluid method, we use the capillary rise in a tube as a test case in the first substep because of the existing analytical solution presented in Sect. 2.4. This first substep is necessary to check the capillary forces in the simulation and the wall interaction as well as the behavior of the VOF phase under gravitational influence.

3
The second substep is needed to check the droplet shape as a function of different contact angles as well as the surface tension between air and liquid. Therefore, a single droplet is initialized on a plate and on a single fiber. The results are compared to the results given by the surface evolver.
Since the surface evolver and the analytical solutions are only valid for very simple cases, we use a two-step validation for more realistic scenarios in order to analyze the performance of the WDM method in Sect. 5. This means that after the validation of the VOF model, this model is used to validate the distance map method on more complex structures in turn. Table 1 summarizes the physical and mesh properties of the validation cases of the VOF method.

Capillary Rise in a Tube (Analytical Solution)
This setup is used for validating the VOF simulations with respect to capillary forces and wall interaction under gravitational influence by looking at the steady state. Specifically, a 3D capillary tube with a diameter of 2 mm is set up in the VOF model of the commercial CFD software StarCCM+. Detailed information on the setup is given in Table 1. A 2D section of the resulting height of the capillary rise is shown in Fig. 5.
For stabilizing the simulation, the water reservoir was simplified to a stagnation inlet with no static pressure. Therefore, this inlet represents the free surface of the water reservoir. In order to imitate the open end of the tube, the boundary condition was set to atmospheric pressure. The rising height of the liquid in the tube is governed by the intermolecular forces between the liquid and the surrounding wall of the tube as well as between gas and solid phases and between gas and liquid. The resulting height of the liquid is governed by the equilibrium between the gravitational forces and the capillary effects. For numerical reasons, viscosity was initially increased to counteract overshooting and damp unphysical oscillations in the simulation.
Looking at the volume fraction of water ( red = 1 , blue = 0 ) in Fig. 6, it can be seen that there is a sharp interface between water and air. This is important for interpolating the resulting shape of the free surface, as well as measuring the rise height. The resulting height of 7.6 mm is slightly larger than the analytical solution given by Eq. (5) in Sect. 2.4, yielding h = 7.4 mm . This is due to the fact that semi-empirical equations are used for the capillary forces acting on the VOF phase which tend to map the contact angle slightly larger. The difference in the resulting height results in an error in the contact angle of less than one degree. Considering the resolution, it becomes clear that the contact angle cannot be reproduced exactly. Accounting for these circumstances, the results are in good agreement.

Droplet Shape Validation (Surface Evolver)
Since the coalescence of droplets deposited on fibers is strongly dependent on the wetting properties, it is the shape of the single droplets that must be checked. Detailed information about the physical and mesh properties of the test cases is given in Table 1. For checking the droplet shape depending on the surface tension and the contact angle, two simple test cases were chosen. The first one is a droplet sitting on a plate (Sect. 4.2.1), while the second case is a straight fiber with a droplet sitting exactly in the middle (Sect. 4.2.2).

Single Droplet on a Plate
We assume that no gravity acts on the fluid, only surface tension and volume conservation will influence the droplets shape. Thus, the droplet will build a spherical cap (Berthier , Brakke 2012).

Geometric method
As stated in Sect. 3, all parameters for the weighting exponents E( ) were found using a least-squares fit. The analytical spherical cap solution was used as a target function in this least-squares fit, making the geometric method produce results as similar as possible to the spherical caps.
Finally, the resulting fluid shape using the geometric method shows only marginal error when comparing to the spherical cap shapes, see Fig. 6. In particular, droplets with larger contact angles form an almost perfect spherical cap.

VOF model
The droplet was initialized as a sphere with a diameter of 1.5 × 10 −4 m , the center exactly on the plate. The resulting shape, depending on the contact angle, is shown in Fig. 7. The 2D slices on the left-hand side of Fig. 7 show the volume fraction of liquid phase ( red = 1 , blue = 0 ).
The volume fraction is characterized by a sharp interface between the phases as before. Considering the limited resolution of the droplets, the contact angles are reproduced correctly although the contact angle at 60 • seems slightly larger using the VOF model in comparison with the solution obtained using the surface evolver.

Single Fiber with One Droplet
A spherical droplet is initialized with a diameter of 1.5 × 10 −4 m , the center exactly in the middle of the fiber. The fiber has a diameter of 4 × 10 −5 m . The resulting shape, depending on the contact angle, is shown in Fig. 8.
The 2D slices on the left-hand side of Fig. 9 show the volume fraction of liquid phase ( red = 1 , blue = 0 ). This volume fraction is characterized by a sharp interface between the phases. Because of the limited resolution of the droplets, the contact angle cannot correspond to the exact value, but is approximated in the simulation. Given the resolution, the contact angle is reproduced well. This can also be seen when the surface of the droplet given by the VOF model is compared to that provided by the surface evolver. The surface evolver's solution is shown as a black line in the 2D slice.
On the right-hand side of Fig. 9, a 3D image of the droplet is shown. The yellow area represents the droplet shape predicted by the surface evolver, while the pink color represents the contour of the droplet predicted by the VOF model. As the 2D view has already shown, the solutions match very well, although the VOF solution appears more choppy due to mesh effects. Since comparing complex structures using the surface evolver is not feasible, the two-step validation with the VOF model is used and described in Sect. 5.

Setup of Numerical Test Cases for the Two-Step Validation
The purpose of the following test cases is to compare the droplet shapes predicted by the VOF model and the WDM approach for more complex fibrous structures. Since the comparison of complex structures is not feasible using the surface evolver, results obtained 1 3 using the WDM approach are compared to results obtained using the VOF model. First, a setup consisting of two crossing fibers will be considered in Sect. 5.1. Next, a small section of a realistic filter structure will be studied in Sect. 5.2. Fig. 7 Droplet on a plate: Contact angle of 120 • /90 • /60 • : 2D slice showing the volume fraction of oil with the surface obtained using the surface evolver with the bold black line (left-hand side) and 3D geometry obtained using the surface evolver (yellow) and the VOF model (pink threshold), right-hand side 1 3

Crossing Fibers
The crossing fiber simulation allows the first comparison of the two methods on a slightly more complex structure. Table 2 summarizes the values used in this simulation.
Both the resulting pressure drop as a function of time and the resulting fluid distribution predicted by the VOF and the distance map model will be compared. In the VOF model, the area of interest is resolved using a finer mesh compared to the rest of the Fig. 8 Droplet on a fiber: Contact angle of 120 • ∕90 • ∕60 • : 2D slice showing the volume fraction of oil with the surface obtained using the surface evolver with the bold black line (left-hand side) and 3D geometry obtained using the surface evolver (yellow) and the VOF model (pink threshold), right-hand side domain. The simulation runs until no shape change occurs any more. Then, an.stl file is saved from the resulting water threshold of 0.5. The result is shown in Fig. 9.
As described in the validation case, a 2D slice showing the volume fraction of water can be seen in Fig. 9 on the left-hand side. The distribution of the volume fraction shows a sharp interface and well-matching contact angles on the fiber surfaces. An isometric 3D view of the droplet shape can be seen on the right side.

Realistic Structure of a Filter Material
In the next step, the same principle is used to investigate droplets in a more complex and realistic filter structure characterized by different contact angles of the fibers at the bottom ( 120 • ) and top layer ( 60 • ) ( Table 3).
In the VOF model, the droplets are not initialized as round droplets, since this becomes geometrically too complex. Instead, velocity inlets for the liquid phase are set up directly on the fibers, which are marked by red circles in Fig. 10. Using a velocity of 0.05 m/s perpendicular to the surface the inlets mimic a constant mass inflow due to deposited liquid particles without necessitating a separate Lagrangian simulation for droplet transport and collection on the fibers.
As the diameter of the droplets increases as a function of time, coalescence starts once droplets touch, see Fig. 11.
The illustration of the coalescing droplets, see Fig. 11, shows the isosurface of the water threshold. The isosurface represents the interpolated area between the water droplet and the air with a value of 0.5. A threshold is also shown in the figure, which shows the volume fraction of water between 10 and 90%. As the threshold coincides exactly with the isosurface, a sharp interface is present.

Results and Discussion
Next, we present and discuss results for the test cases introduced in Sects. 5.1 and 5.2, i.e., the cases of crossing fibers and of a realistic filter structure.  1 3

Crossing Fibers
Simulation results on crossing fibers are shown in Fig. 12. The blue voxels in the left image denote the result using the VOF model; yellow voxels in the right image show the result using GeoDict with the WDM approach. In the middle picture, the magenta voxels mark the overlapping voxels of both results (Fig. 13).

Fig. 11
Coalescing droplets in a complex filter structure obtained using the VOF model Comparing different contact angle configurations (on the 45 • crossing fibers) and checking the overlapping voxels resulted in only 65% overlap on average. This is significantly less than the overlap with the analytical solutions for droplets on a flat surface, but fortunately global flow parameter, such as total permeability, is not as dependent on the volumetric overlap of droplets, see Fig. 14. Considering the mentioned contact angle configurations, there is a maximum of 2% difference in total permeability, when performing a single-phase flow simulation considering both fibers and droplets as solid. Flow simulations have been conducted on both droplet structures using the LIR solver within GeoDict (Linden et al. 2020(Linden et al. , 2015.

Realistic Structure of a Filter Medium
Simulation on a small cutout of a filter medium was performed using the VOF method (left image) and GeoDict using the new geometric method (right image). In both simulations, liquid was injected directly on the fibers at fixed positions. Both simulations were run for 25 iterations, meaning that in each iteration, at each of the 37 inlets, approximately 75 μm 3 of liquid is injected.
The resulting geometry of the interfaces in both cases clearly shows some differences, while the fluid still is distributed similarly. Since the final application of this method is to Fig. 13 Comparison of the flow permeability on the 45 • crossing fibers obtained using the VOF model and the WDM approach with different contact angles (on red/green fiber) Fig. 14 Comparison between the coalescing VOF droplets (left-hand side) and the coalescing droplets using the WDM approach (right-hand side) in the realistic structure simulation simulate coalescence in filter media, the induced pressure drop of the deposited fluid is of major interest. Simulating air flow through all resulting liquid laden structures at a constant volumetric flow rate yields very similar pressure drop values when comparing results using both simulation methods, see Figure 15. This period will describe building up of fluid volume in the filter without drainage effects. The difference in slope at about 25% might be due to the onset of drainage, which cannot be described by the WDM approach. Another reason might be that the interface in StarCCM+ becomes increasingly blurred above 25%. Again, flow simulations have been conducted on both filter structures using the LIR solver within GeoDict (Linden et al. 2020(Linden et al. , 2015. Concerning computational times on the considered small cutout of a filter medium, the WDM method yields results in the order of minutes, while the VOF results are obtained in the order of days.

Conclusions
A new approach for modeling the coalescence of liquid droplets in a fibrous porous medium was introduced which is based on distance maps. Since for complex fiber systems, the new model cannot be validated by means of analytical solutions, a two-step validation via a VOF model was used. After validating the VOF simulation based on analytical solutions as well as results obtained using a surface evolver, results using the WDM approach were compared to results obtained using the VOF model for the case of two crossed fibers. With an average voxel overlap of 65% the results would not look promising at first, but comparison of generated pressure drop shows less than 2% deviation between the two methods. In the next step, results obtained using the WDM approach were compared to results obtained by the VOF model for a small, but complex fiber structure with different contact angles for the fibers at the top and bottom of the domain. Comparing the location of the liquid phase in the filter structure, the distribution obtained using both models looks similar. Although the surface shape is different, the pressure loss as a function of the volume fraction of the liquid phase is in very good agreement with the VOF reference solution. Further validation of the WDM approach is needed by comparing the pressure drop to experimental results. In order to reproduce the long-term behavior of the flow and transport processes in real filters for liquid aerosol particles, surely the movement of droplets in with the flow and the drainage of the liquid phase must be included. GeoDict already has the ability to simulate particle/droplet movement, such that fluid droplets can be injected with the flow. Simulating filtration of solid dust particles is very similar and has been published within GeoDict some years ago and has been improved since (Becker et al. 2016;Azimian et al. 2018). GeoDict is capable of simulating millions of particles within the flow, since a Lagrangian approach is used instead of simulating via a VOF approach. Simulating the droplets in the flow with StarCCM+ is also possible but yields very high computational effort. Since the VOF method will resolve droplet shape and movement even, while they move in free flow, it is very limited when it comes to a realistic case where there are thousands of small droplets injected with the flow. Drainage of fluid along the fibers is currently investigated using the VOF method, and further development of the WDM approach will follow. Finally, the re-entrainment of liquid-phase particles is also an important mechanism that needs consideration, especially at large flow velocities.