Development of a numerical method for multiphase flows using an electrostatic model in a wire-mesh sensor

ABSTRACT Wire-mesh sensors (WMSs) have been used to measure various characteristics of bubbles, such as volume fraction, size and velocity. The WMS used in this study consists of a grid of nine by nine wires arranged in an orthogonal shape. Transmitter wires are activated sequentially and receiver wires pick up the electric signals and convert them to a current. Multiphase flow, the electric conductivity of each phase, allows characteristic electrical signals to flow and from this information the volume fraction can be calculated. However, deformation of bubbles and the intrusive effects of wires have yet to be considered. To address these issues, we conducted numerical simulations using commercial software. The actual behavior of a single bubble was considered to calculate the electrostatic and two-phase flow simultaneously. First, a bubble was considered stationary in a spherical shape at the center of the WMS. A user-defined function was used to calculate the governing equation of the electric field and then activate electric signals on the transmitter wires. The numerical results were compared with reference values, and the error of a volume fraction was quantified. After that, we considered a bubble rising without a change of shape in the pseudo-dynamic process. However, for the dynamic process, a volume of fluid model was added and a bubble interface that reflects its actual movement was considered. A comparison between the pseudo-dynamic and dynamic processes was conducted using the same positional relations between the wire and bubble surface. When the bubble contacted with and then separated from the wire, differences in the received signal and an overshooting phenomenon were observed owing to deformation of the bubble. Also, the range of locations where the bubble would begin to affect the electric field because of the flattened shape of the bubble was observed.


Introduction
Multiphase flow has relevance to various industrial fields (Akbarian et al., 2018;Huang et al., 2015;Ma et al., 2020;Mosavi et al., 2019;Ramezanizadeh et al., 2019;Yayla et al., 2017). Multiphase flows can be classified according to the spatial distributions of the continuous and dispersed phases, e.g. bubbly or slug flow (L-G), spray flow (G-L), aerosol flow, and different liquid flows such as oil and water (L-L). Among these, bubbly flow, which contains very small sized bubbles that form a dispersed phase in liquid medium, has received particular attention for several decades because of its importance to the nuclear (Kong & Kim, 2017;Paranjape et al., 2010), petrochemical (Gao et al., 2016;Hanafizadeh et al., 2015) and biochemical (Luo & Al-Dahhan, 2008 processing industries. When quantifying the efficiency of industrial reactors that use gas-liquid flow, the characteristics of bubbles, e.g. bubble size, volume fraction and velocity of both phases, play a very important role. Therefore, CONTACT Han Seo Ko hanseoko@skku.edu several techniques have been studied to identify these parameters in flow. If an optical probe is used for measurement, the set-up must be cost and time efficient; however, the results tend to be unreliable, especially in the large scale of fluid. The radiographic method measures the volume fraction of gas inside a liquid using X-rays or gamma rays and is able to produce very accurate images for visualization owing to its advantage of being a representative non-intrusive technique. However, the extremely high energies from nuclear reactions would be harmful to anyone nearby who tried to use this technique (Boyer et al., 2002).
The concept of a wire-mesh sensor (WMS) was devised around 30 years ago (Johnson, 1987). It consists of perpendicularly placed wires at a specific distance that form a grid or mesh. Electric potential is sequentially supplied to certain wires, called the transmitter wires (TWs), and this makes the form of a square or sinusoidal pulse of a certain width. After that, regularly positioned wires, called receiver wires (RWs), read the electric signals from the TWs. Then, the electric signal is converted to form a current or electric potential and it is stored for a time in a steady state for about 0.06-0.1 ms after activation. After activation of all the TWs, a matrix that consists of the measured current values can be recorded to produce a discrete image, which is based on a parameter called the received signal (RS) that comes from the difference in electric conductivity between the phases. Therefore, the volume of non-conductive or less conductive fluid is detected in inverse relation to the RS. As a result, the volume fraction (Bowden et al., 2018;Frank et al., 2008), velocity (Kanai et al., 2012;Krepper et al., 2005) and size (Nuryadin et al., 2015;Prasser et al., 2001) of the bubble in the measurement system can be verified.
Recent research using numerical simulations that looks at the accuracy of WMS has been carried out. Prasser and Häfeli (2018) studied the electric potential field using a computational method to compare various conversion methods. Using a reference system, they made comparisons between the Maxwell method and the linear method, looking at the accuracy of the results. Validation was also carried out using the exclusion of negative values of the volume fraction (CUT and NO-CUT) to identify the effect of non-physical values on the numerical simulation. However, it was found that there was some overestimation of bubble size and volume fraction using the conventional methods. Zhang et al. (2019) analyzed the electric field and quantified the systematic error in the WMS using commercial software, COM-SOL Multiphysics. Uniform sensitive volume and cuboid zero conductivity were adopted to verify the effects of non-symmetric formation of the electric field inside the system. The authors tried to demonstrate the cause of overshoots using the three-dimensional geometry of the WMS in their research. In recent years, both phenomena, electrostatics and motion of the bubble, have been considered at the same time in the numerical simulations performed with the open-source software Open-FOAM (Clifford et al., 2019). In these approaches, the electric field is solved as a steady-state problem, after which the movement of the bubble is calculated separately. Various sizes of bubble, lateral crossing points and spheroidal shapes have also been studied. It has also been identified that the observed overshoots are related to non-uniformly distributed electric fields in the bubbly flows.
The purpose of the present study is a comparison and analysis of how electric characteristics are affected by multiphase flow, which is represented by different motions of the bubble: 'pseudo-dynamic' and 'dynamic' processes. From these processes, numerical results from the WMS were compared according to the shape of actual and ideal bubbles. Before that, static bubbles with various sizes were located at the center of the computational domain to verify the compatibility of the code adopted in this analysis. Through these steps, the accuracy of the numerical simulation was validated against previous research (Zhang et al., 2019). In the pseudodynamic process, the electrostatic equation was solved and the center of the spherical non-conductive volume was relocated according to the computation time. Previous work studied rising non-conductive volumes with a consistent sphere shape and a constant velocity (Clifford et al., 2019;Prasser & Häfeli, 2018). However, the results in those studies did not reflect the intrusiveness of effects such as deformation of the interface induced by the interaction between the wires and bubbles. Since the actual single bubble shows a change in shape, the measurement from the pseudo-dynamic process will have different results and cause errors in the accuracy of a sensor. Therefore, a dynamic process which started from an initially spherical bubble in quiescent liquid was computed. Simultaneous computation of both electrostatic and hydrodynamic forces in the two-phase flow was carried out. As a result, more realistic bubbles were depicted and analyzed. This dynamic process is necessary to evaluate the measurement accuracy in the WMS. In addition, the pseudo-dynamic and dynamic processes were compared using the measured current from the numerical results.

Governing equation
A combined electric and two-phase system was solved to visualize the movement of the bubbles through the WMS. The numerical simulation was carried out using commercial software, FLUENT 14.5. To import the governing equation of the electric field induced by the distribution of the electric potential from the TW, a user-defined function (UDF) was compiled.

Electrostatics
The software supports the transport equation with a userdefined scalar (UDS) and particular UDF for unsteady, convective and diffusion terms, which could be specified for each of the kth scalar equations. In the multiphase flow including volume fraction, the transport equation with UDS can be expressed as where φ k , k , ρ, α and U are the kth UDS variable, diffusivity of UDS variable, density, volume fraction and vector form of the velocity field of a certain phase, respectively (ANSYS FLUENT, 2012b). In Equation (1), several conditions were designated so that the equation for the electrostatic model could be derived. At first, because there was no volumetric force, named the Coulomb force, which is induced by the distribution of charge density exerting on the fluid around the WMS, the convective term in the left-hand side could be neglected. Also, the source term in the right-hand side was defined as zero to implement the form of the Laplace equation. Finally, the phase for the computation of Equation (1) was only water, so that ρ equaled ρ l and variables related to the volume fraction could be defined as 1. Therefore, Equation (1) could be rewritten as follows: Because the diffusivity was constant in all computational domains, Equation (1) could be expressed as where φ k is the electric potential of the kth UDS variable in the electrostatic model. In other words, Equation (3) could be given by Other research related to the numerical simulation about the WMS (Clifford et al., 2019;Prasser & Häfeli, 2018;Zhang et al., 2019) solved Equation (4) without the unsteady term, which means that the steady state of the Laplace equation for each location of the non-conductive volume was defined in the simulation. However, since the actual working process of the WMS was imitated in the present study considering the time-dependent condition, the unsteady Laplace equation should be implemented, and Equation (4) was used as a governing equation for the electrostatic model. As mentioned in Section 1, bubbles can be numerically visualized by the RS, which represents the ratio between the calculated electric currents in accordance with the existence of the bubble. The RS can be expressed as follows: where I and I 0 are the current from the surface of the RW with (two-phase) or without (single-phase) bubbles, and i and j are the order of the RW and active TW, respectively. In Equation (5), the electric current can be obtained using the surface integration of the electric field strength on the surface of the RW, which is given by where E is the electric field strength [V m −1 ], and can be expressed as

Multiphase flow
To compute the accurate dynamics of the bubble interface, a volume of fluid (VOF) model was adopted with the finite volume method (FVM) (ANSYS FLUENT, 2012a). The volume fraction of a certain phase in the computational domain can be expressed as α q . The subscript q represents the primary or secondary phase in the system. In this study, since only liquid and gas phases are included in the computation, q can be replaced with either g (gas) or l (liquid). The constraints on the unit cell for the volume fraction equation are given by The interface between gas and liquid can be calculated by solving the continuity equation for the volume fraction of the certain phase, as follows: where ρ q is the density of the phase for q. For an incompressible and Newtonian fluid, the equation for momentum conservation should include a parameter about the volume fraction, which is given by where P, g and F s are the pressure, gravitational acceleration and the source term for body force as a vector expression, respectively. In Equation (10), ρ and μ are the density and dynamic viscosity of fluids obtained by multiplication of the volume fraction. These are given by where μ g and μ l are the dynamic viscosity for the gas and liquid phases, respectively. Designated values for the characteristics of each phase used in the numerical simulation are clarified in Table 1. The body force term related to the surface tension in Equation (10) is calculated using the continuous surface force (CSF) model (Brackbill et al., 1992). The CSF model is given by where γ is the surface tension (0.072 N m −1 in this study) and κ is the average curvature of the interface, which is defined as wheren is the unit normal vector, which is determined aŝ n = n/|n|.

Geometry analysis
To perform the numerical analysis for the single bubble flow around the WMS, a geometry was generated for computation. Figure 1(a) shows a schematic of the computational domain for the numerical simulation of the WMS. The total width of the domain (W) was designated from the result of the wall effect in the rising bubble. The total height of the domain (H) was controlled by varying the distances between TWs and walls, treated as free surfaces (h a ); and this was done to verify the wall effect on the electrostatic model, which will be discussed later in this section. Details of the mesh and scales around the wires are shown in Figure 1(b). The diameter of all wires (d w ) was about 0.2 mm. The TWs are located above the RWs and the wires run perpendicularly to each other. The lateral distance between the TWs and RWs (Δz) was fixed as 2.8 mm and the axial distance between the wires (Δy) was defined as 1.3 mm. The side walls were distanced from both RWs and TWs as half of the lateral distance between the wires (Δz/2). Except for the periphery of the wires, all meshes in the computational domain were generated as cube-shaped cells with a length determined by the mesh convergence. To connect the surface of wire having a circular shape and other mesh cells, square areas around the wires were composed of a trapezoidal shape of mesh. The line connecting the square face and the circular face of the wire was divided into regular intervals and it was controlled to analyze the mesh convergence. The meshes in the longitudinal direction of the wire in the square area have a consistent dimension, which is the same as the cube-shaped mesh. Clifford et al. (2019) investigated the wall effect using various ratios between the diameter of the wire and the distance between the upper wall and surface of the TW. The wall effect in the electrostatic model describes the phenomenon where the electrically insulated top and bottom surfaces of a domain affect the numerical results. When there is a narrower space between the wires and the wall, the formation of the electric field can become more distorted, which could induce inaccurate electric signals.
In previous work, there was no considerable transition in the electric signal when h a /d w was larger than 100 (Clifford et al., 2019). As shown in Figure 2(a), the sufficient distance for the wall effect of the electrostatic model can be obtained with h a = 9.15 mm and the ratio between the scales being similar to the value from the referred research (Clifford et al., 2019). As a result, the total height of the computational domain can be determined as H = 50 mm. In addition, other wall effects related to two-phase flow should be considered. Figure  2(b) shows the terminal velocity of 15 mm of the bubble with respect to the width of the computational domain. It could be verified that there was convergence of the terminal velocity with larger width. As shown in Figure  2(b), lower terminal velocity was found when the width equaled 15 mm owing to the shear induced by the no-slip conditioned walls inside. However, a discrepancy of less than 5% was found between the cases where the width was same as 27 mm. Therefore, the width of the computational domain could be determined for negligible wall effects exerted on the rising bubble.
Using a larger number of mesh cells results in more accurate numerical results; however, more calculation time is required. To optimize the mesh size for the VOF   model, various sizes of mesh were tested, from 0.2 to 0.4 mm, and the number of generated mesh cells was varied from 280,000 to 2,300,000 according to the mesh size. Figure 3(a) shows the interface of the bubble obtained from the numerical simulation at 0.03 s of flow time. No clear difference was observed between 0.25 and 0.2 mm. Also, as shown in Figure 3(b), there was no large difference less than 0.5% in the surface area of the bubble interface (S b ) at the sampled flow time between the referred cases. Therefore, it could be concluded that a mesh size of 0.25 mm was appropriate for the numerical simulation of the two-phase flow.
Since the size of mesh for the electrostatic model should also be considered, a smaller size of mesh was used. A convergence analysis was carried out in a steady state with only one active TW in the center of the system. Figure 4(a) shows the converging trend of the electric strength along the y and z = 0 axis at the central location between the TW and RW according to the increasing size of the mesh. From Figure 4(b), it can be seen that the optimum mesh size is around 4,000,000, which shows negligible discrepancy from cases with a smaller sized mesh.

Computational set-up
For other cases where there is movement of the bubbles, the distance between the RWs and the bottom of the domain (h b ) is extended to 23.75 mm to be able to describe the natural characteristics of the rising bubbles. The top surface of the domain was defined as the 'outlet', which has zero gauge pressure, and a 'no-slip' condition was applied to other surfaces in the computational domain. Half of the system was designed and symmetry was used about the plane at x = 0 to reduce unnecessary computational costs. To depict the motion of the rising bubbles relatively accurately, the maximum time step size was set to 0.1 ms; this is only defined in the computational process for hydrodynamic phenomena. The SIMPLE scheme was used to couple the pressure with velocity. For spatial discretization for the convergence of the solution during the calculation, the pressure and momentum were computed with PRESTO! and a second order upwind scheme, respectively. The first order implicit was used for the transient formulation. The volume fraction in each cell induced by the advection of fluid was calculated by adopting the PLIC. To describe the actual situation of sequential activation of the TWs in the WMS system, a time-dependent boundary condition was coded into the computation. Two different sorts of boundary condition were considered for the numerical simulation. Tested signals of pulses were named A and B. These are illustrated in Figure 5, where the horizontal axis and the vertical axis represent the time and applied electric potential, respectively, to obtain the RS according to the TWs. At first, for pulse A it was expected that the formation of a negative electric potential signal would describe a more realistic electric signal from the WMS, as depicted in the experimental research of Prasser et al. (Figure 5a). The electric field inside the conductive phase in the actual WMS immediately disappeared when the wire was disabled. However, this phenomenon could not be depicted to imitate the result of operation of the WMS without some numerical modification. Therefore, pulse B, shown in Figure 5(b), can be characterized as the initialization signal (green line). In this step, the electric potential values in the computational domain were reinitialized after a steady state of the current induced by the previous TW had been reached.
Because the accuracy of the numerical results with a bubble should be considered not only with the current but also with the electric field around the wires, the distribution of the electric potential in plain water was verified. Figure 6(a) illustrates the location of the  line drawn at x = 0 mm and z = −1.5 and 1.5 mm to identify the asymmetric distribution of electric potential. Through the measurement of the electric potential, the ratio between the electric potential for lines α and β was found. Asymmetry was clearly observed for the boundary conditions of pulse A, as shown in Figure 6(b). Because denser electric potential was distributed in the space between the TWs and RWs, the weaker effect of the residual of the electric potential could be verified. Therefore, there was almost no difference in a certain local range. However, the negative boundary condition of the previous TW induced an increase in the ratio for the y direction, which is not negligible. On the other hand, in the case of pulse B we see a perfectly symmetric distribution of the electric potential, which means that there is no dependency between the activation of other wires. Therefore, it was determined that pulse B was more appropriate for accurate computation and so it was adopted for the numerical simulations discussed in Section 4.

Computational procedure
The numerical simulations of a moving bubble were carried out for two processes: pseudo-dynamic and dynamic. Figure 7 shows the flowchart of the computational procedure for each process. In the pseudo-dynamic process, the initial location of center of the bubble (y o ), expressed as a sphere with non-conductive cells, was defined at the initial time (t o ) and it varied according to computation time (t i ), as shown in Figure 7(a). The rising velocity with defined function (v p ) was modeled to have no change of the interface, but the location of the center of the bubble for every time step (y i ) was considered. The function of velocity according to time could be found from the experimental results. This method has been used in previous work to investigate the electric characteristics of the WMS. The process for the pseudodynamic case can be summarized as follows: (1) A sphere volume of the gas phase (α = 1) was defined inside the computational domain. (2) The first TW was activated (φ = 1 V) during the time of a pulse width. (3) In the liquid domain (α = 0), the unsteady Laplace equation was solved. However, in cells with a gas phase, there was no computation. (4) After the pulse width had finished, the electric potential value for all cells was set to zero. (5) If the residuals for electric potential did not satisfy the designated criteria (1.25 × 10 −4 ), steps (2)-(4) were repeated. (6) If the flow time reached the designated end time, the calculation was terminated.
However, this method has limitations in describing actual two-phase phenomena in the system. Therefore, the dynamic process should also be considered. To do this, the calculation was separated into two parts for efficient computation: one part for only hydrodynamics with two-phase flow (blue box in Figure 7b) and the other with both phenomena (red box in Figure 7b). First, the numerical simulation for only two-phase flow was performed with the designated size of time step (Δt 1 ) and residual (10 −3 ) from the initial environment for t a . When the bubble was located around the RWs, the numerical data were exported and the UDF was compiled to implement the electrostatic model and boundary conditions. After that, with the exported data including information about two-phase flow, the numerical simulation was performed to check the change in the RS. If there was a maximum difference larger than 1% in the RSs, the electrostatic model started to be solved in earnest and the flow time at this moment was defined as t b . However, if the difference was smaller than 1%, the computation of two-phase flow was continued until the difference in signals larger than the criterion was observed. In addition, the size of time step (Δt 2 ) should be changed to be 10 times smaller for operation of the WMS. With these conditions, the electrostatic model for two-phase flow could be solved at the same time. The process for the dynamic case can be summarized as follows: (1) A sphere volume of the gas phase (α = 1) was defined inside the computational domain.
(2) The models for two-phase flow were computed until the flow time (t i = t b ) when the volume of the bubble affected the RS (blue box in Figure 7). (3) If the difference between the maximal RS for two phases and that for plain water was larger than 1% after compiling the UDF, the simulation including the electrostatic model was initiated. (4) The models for the two-phase flow and electrostatics were computed at the same time, according to the phases (red box in Figure 7). (5) If the flow time reached the designated end time, the calculation was terminated. Figures 8 and 9 show images for the rising motion of the bubble in both processes according to time. It can be observed that there are distinct characteristics at the interface. The bottom part of the computational domain was extended by 30 mm to calculate the natural motion of the rising bubble. Comparing the numerical results between the two cases, the initial center of the bubble was defined at 33.25 mm from the center of the domain, where the bubble was located as far as possible from the RW.
The location of the bubble while initiating the electrostatic model was investigated to verify that the bubble is independent of the electric phenomena. It was found that there was a distance of around 8.5 mm between the origin of the computational domain and the center of mass for both processes; this shows that the maximum difference is less than 1% with the RS from plain water.

Validation
The results of the numerical simulation were validated against a previous study (Zhang et al., 2019) to verify the numerical accuracy of our set-up. An RS was obtained using the same values for bubble diameters of 9, 15 and 21 mm, which was modeled as a sphere made up of nonconductive cells. Images of the defined non-conductive volume are shown in Figure 10.
It is important to measure the sectional area populated by the bubble using the results from the WMS. Two parameters, i.e. α avg_true and α avg , are introduced to verify the characteristics of the bubble. The parameter α avg_true is given by where S total and S non_cond are the total sectional area and non-conductive sectional area, respectively. The parameter α avg is the averaged volume fraction found by the RS for the total sectional area using the linear conversion method, which can be obtained by where i, j and n are the sequence index of each TW and RW and the number of wires, respectively. The numerical simulation for validation with a timedependent boundary condition was carried out. Figure  11 shows the distribution of the electric potential from the TW for the computation time when considered as a steady-state system. It can be seen in Figure 11 that the normal vector of the distributed surface for the electric potential was orthogonal to that of the bubble interface, which means that the conditions for insulation on the bubble were satisfied during the computation.
From the computed results, the current was measured for each RW. This can be used to obtain the RS with the calibrated current from plain water. Figure 12 shows the reduction in current when there is a bubble 15 mm in diameter. As shown in Figure 12(a)-(e), the current value at 0.46 ms, when the center of the TW was activated, decreased owing to the bubble having a non-conductive volume, which did not allow the distribution of electric potential. Figure 13 shows the error between α avg and α avg_true . The numerical conditions in this study not only model the system of the WMS well, but also help us to obtain more accurate results by reducing the error relative to the reference. Figure 14 shows the numerically visualized bubble from the WMS. In this figure, the first and second lines are the visualized images of the volume fraction from the numerical results of Zhang et al. (2019) and the present study, respectively. It is of note that there was a smaller amount of overshooting in cells outside the bubble and the conversed volume fraction inside the bubble was much smaller, by about an order of three. It can be inferred that these two dominant differences induced a more accurate numerical result. As a result, it was confirmed that the numerical conditions in this study made it possible to describe the exact characteristics of the dispersed phase; therefore, the cases that will be discussed in Section 4.3 could be modeled.

Bubble motion
Both pseudo-dynamic and dynamic processes were calculated using the electrostatic model with the computational conditions tested in the previous sections. Figure   Figure 10. Computational domain including a bubble with diameter of (a) 9 mm; (b) 15 mm; (c) 21 mm. 15 shows the contour of the electric potential from the activated TW with 9 mm bubble diameter in the pseudodynamic ( Figure 15a) and dynamic (Figure 15b) processes at a certain moment. Different formations of the electric potential distribution around the bubble can be observed. In other words, there would be a different current. This was found by integration of the electric field on the surface of the RW, according to the calculation method regarding the motion of the bubble.
There was a need to reflect the intrusive effect induced by the WMS, which has been investigated for over a decade (Ito et al., 2011;Mukin, 2016), to depict more accurate bubbles inside the system. In particular, some researchers have investigated the break-up phenomenon caused by the WMS using the visualization method (Prasser et al., 2001). Although there was a different regime during the deformation of the bubble owing to physical characteristics such as bubble size, wettability of wires and hydraulic pressure induced by depth of location, it is obvious that the interaction between the wires and the bubble caused a strain on the interface. Therefore, a VOF model needs to be included in the numerical simulation to calculate more accurate results for two-phase flow in the WMS system.
For the dynamic process, the solved flow field around the rising bubble could be expressed as the vector field, as shown in Figure 16. Because of the rising bubble, buoyancy driven flow was calculated around the interface. It  could also be observed that there were vortices behind the bubble due to the pressure difference.
The velocity of the rising bubble under a dynamic process is compared with the experimental results in Figure  17. The velocity can be obtained from the local center of mass of the rising bubble. Only the vertical direction was considered and the mass of gas in each cell was used for the calculation. As shown in Figure 17, the bubble began to rise and there was a linear increase in rising velocity until 10 ms. However, a fluctuation in the velocity was found after 15 ms owing to the change in bubble shape. The terminal velocity from the numerical simulation was measured to be around 0.2121 m/s for a discrepancy of less than 1% from the experimental results. This is slightly lower than the velocity of a rising bubble in the absence of a WMS (Bozzano & Dente, 2001;Liu et al., 2016).  As a result, the reliability of the numerical simulation with single bubble flow could be verified by the excellent agreement with the experiments. Figure 18 is a discrete image of the detected bubble with a diameter of 9 mm in the dynamic process at y CM ≈ 0. As the interface of the bubble in the dynamic process forms a flattened rim against the WMS, a smaller volume fraction can be observed around the center of the plane. To verify the accuracy of the measured results, the average of the volume fraction obtained using Equation (16) was compared with the results of postprocessing in Figure 16. Velocity field around a 9 mm bubble at 0.06 s.  commercial software, and the discrepancy between the two values was found to be around 7%. In other words, even though deformation of the bubble had to be considered in the simulation, the numerical procedure in this study produced excellent results and precisely measured the spatial characteristics in multiphase flow.

Numerical comparison
To compare and analyze the differences in the calculations of the bubble interface between the two processes, four specific cases were considered: contact with receiver wire (CR), contact with transmitter wire (CT), separate from receiver wire (SR) and separate from transmitter wire (ST). Figure 18 illustrates a moment in these situations, with the pseudo-dynamic process on the left side and the dynamic process on the right side. Even though these cases were set up as the same situation, it can be seen that there is a different center of mass. The details of the center of mass and the aspect ratio found in the deformed bubble are clarified in Table 2, where the origin of the coordinate for the absolute distance is defined in the medium of both wires. The aspect ratio (AR) can be obtained from the value of the horizontal length divided by the vertical length of the bubble in the WMS. The discrete images were obtained at a time that characterized the situation. Figure 20 shows the numerically obtained discrete images of the volume fraction of the bubble and the measured RS according to the dimensionless z direction in the pseudo-dynamic (dashed blue line) and dynamic processes (dashed red line). As shown in Figure 20, there was no detected volume fraction which was larger than 0 in Figure 20(a) and (d). It can be demonstrated that only negative values of the volume fraction were induced by overshooting of the RS, as plotted at the bottom of the figure. The same phenomenon was also observed in the cases with the dynamic process. In addition, a smaller distance between the interface and surface of the RW caused larger overshooting values. Therefore, RS was the largest at rec 1 , which was also verified by Figure 21.
Because the elongated bubble between the wires was detected at the moment of situation CT, a larger area was affected, as observed in rec 1 . In Figure 20(b), there is overshooting in the pseudo-dynamic process in rec 2 . It can be seen in Figure 19(b) that the RW was not totally wrapped around and just touched the interface. On the other hand, there was normal detection of the volume fraction in the case of the dynamic process owing to the non-conductive gas cells existing in the space between the wires. Also, owing to the deformation of the bubble, the asymmetry of the volume fraction of the discrete images along the z * direction could be clearly observed in the dynamic process. There was a little overshooting from rec 3 in the pseudo-dynamic process. From the results of the numerical simulation, it should be noted that there was little effect of the bubble around the RW, even if there was no contact.
There was a larger difference between the two processes in the SR case (Figure 20c). Measurement of the non-conductive volume can be carried out in this case with the pseudo-dynamic approach. There was no overshooting, owing to the totally surrounded TW in that process. However, a little overshooting can be observed in the case of the dynamic process, which was induced by partially covering parts of the wire. Also, even though there was non-conductive volume between the wires, overshooting was not detected owing to the narrower formation of the bubble induced by viscous effects from the interface. The details of the quantified values are plotted in Figure 21(b) and (c).
The continuously varying RS was investigated according to the center of mass of the bubble. Figure 22 shows the RS for different processes when the centered TW was activated in a steady state. In Figure 22(a), which shows results from the pseudo-dynamic process, overshooting can be seen in the CR and CT cases from rec 1 and rec 2 , respectively. After that, continuously small values were measured in the range between CT and SR. It is known that the RS from rec 2 started increasing earlier than that from rec 1 owing to the constant curvature of the bubble. Finally, overshooting was also shown at ST; however, the value in that situation was smaller than that of the previous one.
More complex phenomena were observed in the dynamic process (Figure 22b). The overshooting from rec 1 and rec 2 occurred at CR, which was the same as in the pseudo-dynamic process. On the other hand, the RS from the other RW had already decreased at CT. After that, the small value for RS continued until the center of mass of the bubble was located at 0.2 mm. Overshooting occurred at SR only from rec 1 , while rec 2 showed slight signs of overshooting at ST.
In summary, obvious differences were observed from the calculation of the interface. As a result of the horizontal deformation of the bubble, which was induced by hydrodynamic parameters, the electric signals obtained were changed in various ways. In addition, the deformed shape could result in a shorter period for the rising bubble, with a larger area for sensing. Therefore, if the numerical simulation is performed on the WMS system, for accurate results, it is necessary to describe the intrusive effect using the VOF model considering the characteristics of the system.

Final remarks
In this study, a numerical method for the simultaneous computation of multiphysics was developed. First, the designated conditions mimicking the actual operation of a WMS were identified and visualized images of a static bubble were obtained. Second, to ensure the numerical reliability of the methods used in the simulation, our results were validated against previous work and a reduced error of around 1% was obtained. Third, simulations of a bubble under motion were carried out to verify the critical differences induced by considering how the interface during sensing affected the results. Our pseudo-dynamic process only considered updating the center of mass at every time step, while the spherical shape was assumed to be constant. A dynamic process was also simulated with hydrodynamic phenomena using  a VOF mode. Through the numerical results from characteristic situations for both processes, it was verified that there are several differences in terms of deformation of the bubble, as follows: • Similar electric signals were observed when the bubble came into contact with the RW or detached from the TW.
• Depending on how the bubble surrounded the wire, there was a difference in the RS. Also, the asymmetric discrete images were found to reflect the actual bubble. • As the shape of the bubble is flattened by the WMS, the range of center of mass affecting the RS was much smaller in the case of the dynamic process, and the RW located far from the center also showed variance in the signal. As a result, when modeling a WMS numerically using commercial software, time-dependent changes in the bubble interface should be considered in the computational procedure, along with electrostatic phenomena. Our approach solved these at the same time to produce more effective results.
However, there were some limitations, in that only numerical results were analyzed without verification about the deformation of the bubble. In general, the intrusive effect induced by the interaction between the bubble and sensors for measurement should be considered to improve the accuracy of the WMS. This requires experimental results showing several dynamics of the bubble, such as break-up and coalescence, which can significantly affect the electric field around the activated wires. In addition, it was realized that there is a complex relationship with the path instability of the bubble. In future research, the flow field around the WMS could be analyzed using visualization methods such as a particle image velocimetry. Furthermore, the effect of the dynamics of the bubbles from the actual WMS with the electric circuit could be analyzed in detail to enhance the accuracy of the measurements in two-phase flow.