Dynamic simulation and experimental verification of foam transport in porous media based on level set method

The properties of the foam fluid in porous media depend on its dynamic structure. However, few studies have focused on the dynamic structural change in foam. In this research, a method based on the level set method is proposed for the simulation of aqueous foam transport in porous media. As a first step, we focus on the Jamin effect, and a concise relation is established which reveals that the resultant pressure is proportional to the capillary pressure. Second, the snap‐off process of foam in porous media is studied, and both the two stages, including partial disintegration, are revealed based on the results. Third, the conclusion that foam automatically selects a larger channel in which to continue flowing is obtained by studying the foam selectivity. Finally, experiments are carried out to verify the simulation results using different models. The simulation results are in good agreement with the experimental results. These results are expected to be helpful for further understanding the foam characteristics and their effective applications.

Previous studies 4,5 have shown that the ability of foam to trap gases in porous media and the stability of foam are very important for its successful application. This is because the stronger the ability of foam trap gas, the better the control of gas mobility. The stronger the foam stability is, the better the foam plugging will be; the higher permeability layer will be blocked, the production profile will be adjusted, and the volume will be increased, so that oil recovery can be enhanced. Therefore, the dynamic migration mechanism of foam transport in porous media is worth exploring. Many studies, mainly including experimental research and model simulation, have been carried out on the mechanisms of foam transport in porous media.
Displacement experiments using stainless-steel tubes filled with silica sands or glass beads consistently have been used to explore the law of the overall migration of foam transport in porous media. The results from such studies have shown that [6][7][8][9] because of an added capillary resistance, foam can reduce the gas mobility by trapping gas in porous media. Therefore, foam fluid has great potential for profile control and resolving interactions between interfaces. In recent years, CT (computed tomography) has been used to study the dynamic change mechanism of foam during core flooding. [10][11][12][13][14][15] The texture and distribution of foam in porous media have been thoroughly studied, and some results have been obtained. For example, the liquid injected after a foam fluid does not follow the path of the mobile gas in the foam. Pore-level investigations [16][17][18][19][20][21][22][23] consistently have been carried out to explore the law of the local migration of foam in porous media, including the generation, coalescence, and rupture of foam. Lamella leave-behind, gas bubble snap-off, and lamella division are the mechanisms of foam generation. The breakage of pseudoemulsion films 24 can be used to explain the cause of foam decay.
Two types of models are used for simulating the overall migration of foam transport in porous media, namely local equilibrium models and population balance models. 25 Local equilibrium models 10 are usually established based on the resistance factor, which is used to explain the reduction of gas mobility. Population balance models 26,27 are established based on the balance of foam generation and coalescence; furthermore, foam texture has been introduced into these models. The two types of simulation models 10,28-33 were built by introducing some indexes, and reducing the gas relative permeability and increasing the gas viscosity were two main methods used to alter the mobility.
Although foam texture is considered in the population balance models, foam density, the key index of foam texture, was also proposed based on some empirical parameters, such as the flow rate of the flowing gas and the liquid. The models for the micromechanism of foam are mainly based on mechanical analysis. [34][35][36][37][38][39][40][41] Visual simulation results of the foam micromigration mechanism in porous media, such as the micromorphology mechanism, the mechanism of interaction between bubbles and the Jamin effect, are rare. Saye and Sethian 42 proposed an effective simulation model based on foam texture, which is a multiscale model built with describing the evolution of foam membrane rearrangement, drainage, and rupture. Unfortunately, this model was studied in an open space, and the results were not suitable for a narrow space. Nevertheless, we can learn from the level set method mentioned in this model to achieve foam interface tracking.
The simulation of foam transport in porous media is difficult to execute because of the complex boundary conditions and changeable foam interface required. In recent years, a method based on fractal theory was proposed for the evaluation of foam in porous media but failed to give dynamic simulation results. 43,44 There are many ways to trace the foam interface, such as by using the VOF (volume of fluid) method 45,46 and the level set method. [47][48][49] The level set method is advantageous to use when dealing with the boundary conditions, especially with a change in topological structure, such as a bubble splitting into two. Therefore, using level set method to simulate foam flow can directly understand the change of foam morphology and structure, and guide the field application. However, there are few reports on the simulation of complex structure channels using this method, especially on the verification of experiments.
In this paper, we attempt to establish a model of foam fluid based on the level set method. Several simplified porous media structures are designed, and some major mechanisms, such as the Jamin effect, snap-off, and selective transport, are analyzed. As a verification, experiments were carried out by using a microscopic glass model with the same pore structure as the simulation. Our results are of potential interest regarding the dynamic study of foam transport in real formations based on digital core.

ON THE LEVEL SET METHOD
In the standard level set method, the level set function (⃗ x,t) is defined as a signed distance function using Equation (1) where i is the interface of the foam, (⃗ x,t) > 0 represents one side of the interface, (⃗ x,t) < 0 represents the other, and (⃗ x,t) = 0 represents the interface.
To ensure the conservation of phase mass in the calculation process, it is usually necessary to improve the signed distance function. 50,51 In general, the Heaviside function is used to replace the initial signed distance function as the new level set function where corresponds to half the thickness of the interface. Let the new (⃗ x,t) = � (⃗ x,t).
Normal ⃗ n and curvature ⃗ k can be obtained from our phase field function as The level set function is affected by the fluid velocity field: where ⃗ u is the velocity field, and t is the time. In addition, after a certain number of iteration steps, (⃗ x,t) will no longer be a distance function. Generally, to maintain (⃗ x,t) as a distance function, a reinitialization equation proposed by Sussman 52 is adopted. This is accomplished by solving the following equation for a steady state: where 0 is a level set function at any computational time, and is the virtual time in a reinitialization step and denotes the smoothed sgn( 0 ) function with appropriate numerical smearing to avoid any numerical difficulty.
The density and viscosity vary smoothly over the interface by letting where 1 and 2 and 1 and 2 are the dimensionless densities and viscosities of the two fluids, respectively.
The Navier-Stokes equations including surface tension are where the continuity condition is ∇⃗ u = 0, I is the unit matrix, and ⃗ P is the pressure field. ⃗ F sv is given by transforming surface tension to a volume force spread over a few layers of cells using the model (CSF) proposed by Sussman 44 : where is the surface tension, and is the Dirac function.
The simulation of foam was based on the abovementioned principles of the level set method.

DISCUSSION
The mechanisms of foam transport in porous media are difficult to study because of the complex boundary conditions. In this paper, we tried to simplify this problem, and the complex boundary problem of porous media is separated into a combination of some representative models with simple boundaries. Therefore, we can achieve the goal from part to whole.

| Jamin effect
The Jamin effect is the gas resistance effect caused by the deformation and stretching of gas bubbles through a narrow pore throat. This is the principal reason for many special properties of the foam, such as blocking and diversion ability. Therefore, a simple pore throat model is established to simulate the process. The geometric model (2-D) is as follows: As shown in Figure 1, the geometric model is composed of two narrow channels: The first is a circular tube with a width of 1 mm, and the second is a circular tube with a width of 0.5 mm. We studied foam bubbles with a radius of 0.45 mm near a pore throat. The gas phase is nitrogen, and the liquid phase is water; furthermore, the surface tension of the foam is 5 mN/m. The flow velocity at the inlet is 1 mm/s, and the pressure of the outlet is 0 Pa. The capillary number is 2 × 10 −4 . In all the simulation models in the work, the compressibility of the gas phase was neglected.
As shown is Figure 2, we used triangular mesh generation method in this work. The number of elements is 681, and the minimum element quality is 0.5494. Average element quality is 0.9322. Element area ratio is 0.05693. Maximum growth rate is 2.325. Average growth rate is 1.305. Edge elements are 141. Vertex elements are 12.
The simulation results are shown in Figure 3, and the foam deforms as it passes through the pore throat, stretching at 0.1 seconds and passing through the pore throat within 0.8 seconds. Through the change in lamella (Figure 4), we can clearly observe this process.
During the driving period, the lamella gradually changes from a left depression to a right depression. There are a stack jump and asymmetric jump driven by the minimum surface area. In the subsequent migration, the direction of the curvature of the lamella also changes.
For comparison, we added two cases with different contact angle. One is 60 degree; another is 120 degree. Some simulation results are shown in Figure 5 as follows: As shown in Figure 5, the change of contact angle will affect the direction of the lamella or the dynamic shape of the foam, but it will not affect the advancing speed of the foam. In both cases, foam completely passes through the throat within 0.80 seconds, which is the same as the case above.
As shown in Figure 6, we selected point 1 (green) as an analysis point to study the pressure variation as the foam moves through the pore throat, and the pressure change curve and corresponding foam structure are shown as follows. Figure 6 shows the variation in pressure at the pore throat with time, and the shape of the foam is added on the graph. The pressure at the pore throat changed dynamically as the foam passed through the throat because foam will change its shape according to the difference between the radius of the pore and the throat. Figure 6 shows that the maximum pressure corresponds to the maximum deformation (half the volume of the bubble entering narrower channel); once the foam passed through the throat (0.11-0.26 seconds), the pressure dropped rapidly, even when a negative pressure occurred because of the change in the lamella direction. The pressure at the pore throat is affected by the characteristics of both the foam and pore throat; therefore, the influence of the surface tension and throat size on the pressure was thoroughly analyzed next. Figure 7 shows the variation in pressure at the pore throat as a function of time for different surface tensions. As shown in Figure 7, the greater the surface tension, the greater the pressure. In fact, the greater the surface tension, the harder the foam stretches, and the greater the driving force of the foam passing through the same throat. Besides, after 0.6 seconds, the pressure at the observed point is sometimes negative and sometimes positive with the variation of the surface tension. According to Bernoulli's principle, the total pressure in fluid flow is constant, and the total pressure is equal to the sum of dynamic pressure and static pressure. The greater the surface tension, the greater the dynamic pressure of the foam due to the change of the shape of the lamella under the same flow rate. In order to ensure the balance of the total pressure, the static pressure will be smaller, and negative pressure is more likely to occur. Figure 8 shows the curve of the relation between the surface tension and the maximum pressure corresponding to the maximum deformation.
The curve in Figure 8 shows that the surface tension is proportional to the maximum pressure P of the pore throat. We can express this relation with the following formula: ZHANG et Al. Figure 9 shows the variation in pressure at the pore throat as a function of time for different ratio of the pore throat's radii. As shown in Figure 9, the smaller the throat size is, the greater the pressure is. In fact, the smaller the throat size, the greater the degree of deformation and the capillary force, the greater the energy consumed when foam passing through the throat, and the greater the resistant force of the foam passing through the same throat produced by Jamin effect. Figure 10 shows the curve of the relation between the reciprocal of the throat radius and the maximum pressure corresponding to the maximum deformation. There is a linear relationship between the maximum pressure at the throat and the reciprocal of the throat radius. We can express this relationship with the following formula: From Equations (9) and (10), we can obtain that: In fact, as we all know, the capillary force formula is: So, we can deduce that the maximum pressure is proportional to the capillary force.

| Snap-off
When foams are transported in porous media, snap-off occurs easily when they pass through narrow channels. Here, this physical phenomenon was modeled and simulated, and the mechanism of foam snap-off was analyzed. As shown in Figure 10, a three-dimensional geometric model was established by rotating a two-dimensional surface.
As shown in Figure 11, the maximum radius of the model is 3 mm, the radius of the narrow pore throat is 0.1 mm, and the initial radius of the foam is 2 mm. The gas phase is nitrogen, and the liquid phase is water; furthermore, the surface tension of the foam is 10 mN/m. The flow velocity at the inlet is 20 mm/s, and the pressure of the outlet is 0 Pa. The capillary number is 2 × 10 −3 . The simulation results are as follows: As Figure 12 shows, at the beginning of the simulation, the foam gradually compressed with increasing pressure; when the pressure exceeded the resistance from the Jamin effect at 0.300 seconds, the foam deformed. The phenomenon of snap-off occurred at 0.660 seconds because of the lamella behind the narrow pore throat. Then, a small part of the foam flowed out as a liquid (yellow circle marked), while most of the foam passed through the pore throat and regrouped into a new foam body by coalescence.
Based on the simulation above, the snap-off process of foam in porous media was divided into two stages. In the first stage (before 0.618 seconds), the foam is partially disintegrated or ruptures at the pore throat, forming small bubbles after passing through the throat. At this stage, the pressure increases gradually along the pore throat, resulting in the deformation of the foam. In the second stage (after 0.618 seconds), the small bubbles that have passed onto the other side of the pore throat form large bubbles or flow with the liquid. At this stage, the series of small bubbles reach a critical distance and do not coalesce.
Besides, we discussed at what flow velocity, the bubble would stop at the narrow throat, which is of great significance for field application. From Figure 13, we can see that when the flow velocity is 17 mm/s, the bubble will be stuck and stop at the narrow throat.

| Selectivity
In the application of foam, selective blocking is one of its most important properties. Here, the selective migration mechanism caused by this property was simulated and analyzed, and a basic analysis of the law of foam seepage was performed. The geometric model (2-D) is shown in Figure  13.
As shown in Figure 14, the geometric model is composed of three parts to simulate coalescence, deformation, and selective transport. The width of the part near the entrance is 3 mm; the narrowest width of the middle part is 0.4 mm; and the part near the outlet is composed of two channels, in which the width of the upper part is 1.0 mm, and the width of the lower part is 1.6 mm. In addition, the radius of each of the two foam bodies is 1.2 mm. The gas phase is nitrogen, and the liquid phase is water; furthermore, the surface tension of the foam is 10 mN/m. The flow velocity at the inlet is 20 mm/s, and the pressure of the outlet is 0 Pa. The capillary number is 2 × 10 −3 . The simulation results are shown in Figure 15.
As shown in Figure 15, after coalescence and deformation, the foam automatically selected the larger channel to continue its flow. Accordingly, foam selectively entering wide passages is the fundamental reason for the temporary blocking property of foam; furthermore, additional resistance would develop because of the Jamin effect.

| EXPERIMENTAL VERIFICATION
The previous simulation results were all based on the theory of the level set method and gas-liquid two-phase flow. However, the reality of the simulation results needs to be verified. Therefore, we verified the validity of the model through microscopic experiments.

| Apparatus
Foam image acquisition is the basis of the analysis of the foam fluid, and its realization requires a highly accurate imaging device. A microscopic visualization device was used for aqueous foam image acquisition. A schematic of the microscopic experimental apparatus is shown in Figure  16.
The glass-etched micromodel used in this study was made by etching a two-dimensional network of pores and pore throats by using a photochemical method. Then, a high-definition camera was employed to accurately record the changes in foam during the processes of foam transport in porous media. The size of the micromodel was 10 cm × 10 cm, the size of the flow area was 5 cm × 5 cm, and three pore throat models were designed, as shown in Figure 17.
In addition, a syringe pump was used to inject liquid (surfactant solution), and a mass flow controller was used to control the rate of gas flow. Then, a series of bubbles were generated before they entered the etching model. A BPR (back pressure regulator) was used to control the back pressure.

| Materials
The surfactant used was SDS (sodium dodecyl sulfate). The gas used was nitrogen with a purity of 99.99%.

| Experimental procedure
Micromodel experiments were performed as follows. First, several pore volumes of SDS surfactant solution (0.3 wt%) were injected at a suitable rate into the pore model using a syringe pump. Next, nitrogen was injected into the model at 10 bars, from a cylinder to the mass flow controller, and a row of bubbles was generated. Finally, the SDS surfactant solution was injected again to force the bubble forward.

| Experimental and simulation results
The experiments and simulations were carried out using the model above. The geometric models (2-D) were established to mirror the glass model. We kept the experiment and simulation processes as consistent as possible by controlling the flow fluid rate.

| Changing diameter model
The experimental results and the simulation results of foam flowing in the tube with a changing diameter are shown in  Figure 17 at different times. The tube with a changing diameter was built as a very small Laval pipe with a narrowest diameter of 0.1 mm and a widest diameter of 1 mm. During the simulation, the initial foam radius was 0.45 mm, with a nitrogen gas phase and a water liquid phase. Furthermore, the surface tension of foam is 35 mN/m. The flow velocity at the inlet is 8 mm/s, and the pressure at the outlet is 0 Pa. The capillary number is 2.3 × 10 −4 . In addition, the wetting angle of the wall is 15 degrees.
As shown in Figure 18A, the shape of the foam transformed with the changing tube diameter. Through the narrowest part of the tube, the foam accelerated and rapidly passed through the channel. Figure 18B shows the simulation results of the foam flowing in the tube with a changing F I G U R E 1 6 Schematic of the microscopic experimental apparatus F I G U R E 1 7 Three pore throat models diameter. When the foam passed through the variable diameter model, the speed was slower in the beginning stage but faster in the middle stage as it passed through the narrowest part of the tube. The foam took 0.50 seconds to pass through the first 1/3 of the total distance but less than 0.20 seconds to pass through the middle 1/3 distance. A contrastive analysis between the experimental result and the simulation result at 0.60 seconds, as shown in Figure 18C, shows that the simulation result is in good agreement with the experimental results.

| Double-cone model
The experimental results and the simulation results of foam flowing in the double-cone tube at different times are shown in Figure 19. The double-cone tube was built to replicate the model used in the simulation, with a narrowest diameter of 0.1 mm. During the simulation, the initial foam radius is 0.50 mm, with a nitrogen gas phase and a water liquid phase. Meanwhile, the surface tension of the foam is 35 mN/m. The flow velocity at the inlet is 10 mm/s, and the pressure at the outlet is 0 Pa. In addition, the wetting angle of the wall is 15 degrees. The capillary number is 2.9 × 10 −4 .
As shown in Figure 19A, the shape of foam transformed as the pore throat changed. Figure 19B shows the simulation results of the double-cone model at different times. In fact, foam is not a simple single-phase or multiphase fluid. It consists of continuous bubbles; because of existence of bubbles, Jamin effect will happen in the narrower channel that will increase the additional resistance to flow. Therefore, with the constant change of foam shape, the Jamin effect is constantly occurring, so the flow has additional resistance all along, resulting in the reduction of flow velocity. A contrastive analysis between the experimental result and the simulation result at 0.30 seconds, as shown in Figure 19C, shows that the simulation results are in good agreement with the experimental results.

| Complex structure model
The experimental results and the simulation results of foam flowing in the tube with a complex structure are shown in Figure 20 at different times. The complex structure tube was built to the same specifications as the model used in the simulation, with a narrowest diameter of 0.2 mm and a widest diameter of 1 mm. During the simulation, the initial foam radius is 0.50 mm, with a nitrogen gas phase and a water liquid phase. Meanwhile, the surface tension of the foam is 35 mN/m. The flow velocity at the inlet is 5 mm/s, and the pressure at the outlet is 0 Pa. The capillary number is 1.4 × 10 -4 . In addition, the wetting angle of the wall is 15 degrees.
As shown in Figure 20A, the shape of the foam transformed as the pore throat changed and was influenced by the Jamin effect. Figure 20B shows the simulation results of the complex structure model at different times, indicating that the shape of the foam changed when it flowed into a smaller pore throat influenced by the Jamin effect; this behavior was in good agreement with the experimental results.
F I G U R E 2 0 Pictures and simulation results of foam flowing in the complex structure tube 5 | CONCLUSIONS In this work, the level set method was used to study the dynamic simulations of foam transport in porous media. When foam passes through a pore throat, the greater the degree of foam deformation, the greater the pressure at the pore throat. The pressure is proportional to the capillary pressure. The snap-off process of the foam in the porous media is divided into two stages: (a) Foam appears partially disintegrated or ruptures at the throat, forming small bubbles after passing through the throat, or (b) the small bubbles that pass through to the other side of the throat form large bubbles or flow with the liquid. Foam automatically selects the largest channel available to continue its flow. Accordingly, foam selectively entering wide passages is the fundamental reason for the temporary blocking property of foam. Finally, experiments were carried out to verify the simulation results using different models. The simulation results are in good agreement with the experiments.