A method for high-precision wetting boundary based on color lattice Boltzmann method

Abstract

The present invention discloses a method for high-precision wetting boundary based on color lattice Boltzmann method, comprising the following steps: step 1, obtaining a solid boundary unit normal vector; step 2, obtaining a target phase interface unit normal vector according to the solid boundary unit normal vector and a set contact angle; step 3, introducing the target phase interface unit normal vector into a formula for calculating interface force, introducing a zero interface force condition, realizing the operation of interface force at the boundary, and completing the introduction of boundary wetting conditions. The present invention is based on the existing strict color lattice Boltzmann model, so that the model is applicable to solid surfaces of arbitrary shapes; improves the existing contact angle introduction mode, so that the model is applicable to contact angles of arbitrary sizes; and introduces a zero interface force condition, reduces non-physical mass migration effects, and improves the stability and accuracy of the algorithm.

Description

Technical Field

[0001] The invention relates to the research field of the internal mechanism of displacement oil production, and in particular to a high-precision wetting boundary method based on a color lattice Boltzmann method. Background Art

[0002] In the past two decades, LBM (Lattice Boltzmann Method) has developed into a heavy solver for simulating complex flow problems. Different from the traditional method of solving macroscopic physical quantities, LBM is a pseudo-molecular method based on microscopic operations and macroscopic kinetic equations. It achieves the solution of macroscopic physical quantities by capturing the equilibrium distribution function of particles. Innately, LBM has the advantages of easy parallelization, easy handling of complex boundaries and easy handling of phase interfaces, making LBM very effective for simulating porous media problems at the pore scale. Many multiphase and multicomponent LBM models have been seen in technical literature. Generally speaking, they mainly include the following categories: color model, pseudopotential model, free energy model, and mean field model. The color model is more suitable for simulating immiscible two-phase flow in porous media because it has the advantages of high numerical accuracy, strict guarantee of mass conservation of each fluid, and good numerical stability within a wide range of viscosity ratios.

[0003] With a suitable multiphase flow model framework, an important issue is the modeling of the wall contact angle. Among the existing two-phase flow wetting boundary processing methods, especially in the field of color models, the most commonly used contact angle implementation method is the "virtual density method". This method imposes a virtual fluid density condition on the solid wall node close to the fluid, so that the color gradient on the outermost fluid node can be calculated, and then the interfacial force that affects the interface movement can be obtained. In the virtual density method, there is a mapping relationship between the contact angle and the virtual density. By adjusting the virtual density, different contact angles can be obtained. However, the contact angle cannot be explicitly and accurately adjusted; in particular, when the contact angle is less than 90 degrees, serious non-physical mass migration will occur, making the results inaccurate and unstable. Taking the research work of Huang et al. as an example, this work can only be applied to the contact angle implementation method of straight boundaries, and there is a large non-physical mass migration problem, which leads to interface slip (Huang H, Huang JJ, Lu X Y. Study of immiscible displacements in porous media using a color-gradient-based multiphase lattice Boltzmann method [J]. Computers & Fluids, 2014, 93: 164-172).

[0004] In order to explicitly adjust the contact angle, Ding Hang et al. proposed a "geometric" contact angle implementation method, in which the contact angle is applied to the boundary as a boundary condition, so that it can be explicitly adjusted (Ding H, Spelt PD M. Wetting condition in diffuse interface simulations of contact line motion [J]. Physical Review E, 2007, 75 (4): 046708); however, this algorithm can still only be applied to straight boundaries, and the range of contact angles that can be studied is small.

[0005] Leclaire et al. further improved this method by using the modulus of the color gradient as a reference to introduce the contact angle, thereby increasing the accuracy of the contact angle simulation. The implementation method of Leclaire et al. still has some non-physical mass migration, and symmetric problems will calculate asymmetric results. There is no algorithm content for calculating the binary zigzag solid boundary, which is reflected in his own paper (Leclaire S, Abahri K, Belarbi R, et al. Modeling of static contact angles with curved boundaries using a multiphase lattice Boltzmann method with variable density and viscosity ratios [J]. International Journal for Numerical Methods in Fluids, 2016, 82 (8): 451-470).

[0006] In summary, the existing "geometric" contact angle implementation method has the problems of being unable to handle curved surface boundary strips, excessive non-physical mass migration effects, and non-physical slippage of the contact line. Especially in the simulation of displacement or fingering phenomena in porous media, these problems will make the simulation results less convincing. Summary of the invention

[0007] Since the existing technology can only calculate the contact angle of a flat interface or a surface with a known curvature, and there is a large degree of non-physical mass migration effect, the calculation results are inaccurate, and asymmetric conclusions appear for symmetric problems. Therefore, it is very necessary to develop a method with ultra-high precision, suitable for a wide range of surface tensions, and able to calculate the contact angle of solid surfaces of arbitrary shapes.

[0008] In order to solve the above technical problems, the present invention discloses a method for high-precision wetting boundary based on color lattice Boltzmann method, comprising the following steps:

[0009] Step 1, obtain the unit normal vector of the solid boundary;

[0010] Step 2: obtaining a target phase interface unit normal vector according to the solid boundary unit normal vector and a set contact angle;

[0011] Step 3: Introduce the unit normal vector of the target phase interface into the formula for calculating the interface force, introduce the zero interface force condition, realize the calculation of the interface force at the boundary, and complete the introduction of the boundary wetting condition.

[0012] In one embodiment of the present invention, in step 1, the solid boundary unit normal vector is obtained according to formula (1):

[0013]

[0014] Among them, n s (x) is the unit normal vector of the solid boundary at the x position,

[0015] c k is the direction of the eighth-order isotropic difference,

[0016] w(|c k | 2 ) is the corresponding weight coefficient,

[0017] δ t is the time step,

[0018] s(x) is the attribute indicator function, and its specific value is shown in formula (2):

[0019]

[0020] In one embodiment of the present invention, step 2 includes the following steps:

[0021] First, according to a given contact angle θ, the solid boundary unit normal vector is rotated to obtain a phase interface unit normal vector, and the angle between the phase interface unit normal vector and the solid boundary unit normal vector is θ; there are two phase interface unit normal vectors, namely n1 and n2,

[0022] n1=(n s，x cosθ-n s，y sinθ,n s，y cosθ+n s，x sinθ)

[0023] n2=(n s，x cosθ+n s，y sinθ,ns，y cosθ-n s，x sinθ)

[0024] Among them, n s，x is the projection size of the unit normal vector of the solid boundary in the x direction of the coordinate axis,

[0025] n s，y is the projection size of the unit normal vector of the solid boundary in the y direction of the coordinate axis,

[0026] θ is the contact angle;

[0027] Then, the target interface unit normal vector is determined according to the Euler distance between the two interface unit normal vectors and the predicted unit normal vector, as shown in equations (3) and (4).

[0028] D1=|n1-n * |, D2 = |n2-n * | (3)

[0029]

[0030] Among them, n * is the predicted unit normal vector,

[0031] n f is the unit normal vector of the target phase interface,

[0032] n1 and n2 are two unit normal vectors of the phase interface,

[0033] D1 and D2 are the Euler distances between the unit normal vectors of the two phase interfaces and the predicted unit normal vector.

[0034] In one embodiment of the present invention, the predicted unit normal vector n * According to formula (5),

[0035]

[0036] in, is the predicted color gradient.

[0037] In one embodiment of the present invention, the predicted color gradient is obtained according to the following method:

[0038] Firstly, the virtual fluid density distribution is introduced into the outermost solid nodes by weighted averaging;

[0039] Then, the predicted color gradient is calculated from the density distribution of the solid node and its neighboring points.

[0040] In one embodiment of the present invention, in step three, the zero interface force condition is introduced by introducing a weighted average virtual point method.

[0041] In one embodiment of the present invention, the zero interface force condition is the gradient of the unit normal vector of the target phase interface of the fluid node at the boundary.

[0042] Based on the existing strict color lattice Boltzmann model, this invention proposes a method for calculating the normal vector of a solid curved boundary or a binary zigzag boundary, making the model applicable to solid surfaces of arbitrary shapes; improving the existing contact angle introduction mode, making the model applicable to contact angles of arbitrary sizes; and introducing the zero interface force condition, reducing the non-physical mass migration effect, and improving the stability and accuracy of the algorithm. The implementation of this series of algorithms constructs a set of highly accurate two-phase flow wetting boundary condition methods, providing a powerful research tool for the field of surface tension-dominated flows.

[0043] The above-mentioned technical features can be combined in various technically feasible ways to produce new embodiments, as long as the purpose of the present invention can be achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] The present invention will be described in more detail below based on embodiments and with reference to the accompanying drawings, wherein:

[0045] Figure 1 is a flow chart of a method for high-precision wetting boundary based on color lattice Boltzmann method in the present invention;

[0046] Figure 2 Schematic diagram of the rotation angle relationship between the unit normal vector of the solid boundary and the unit normal vectors of the two phase interfaces in the first embodiment of the present invention;

[0047] Figure 3 This is a comparison diagram of the non-physical mass migration effect applied to the color grid Boltzmann model in the first embodiment of the present invention and the non-physical mass migration effect between the original geometric model and the virtual density model;

[0048] Figure 4 is a geometric initialization model diagram for simulating the capillary penetration problem in the second embodiment of the present invention;

[0049] Figure 5 is a comparison diagram of the theoretical solution of the capillary penetration problem under different wettability and viscosity ratios in Example 2 of the present invention, the solution of this model and the solution of the Lecalire model;

[0050] Figure 6a and Figure 6b It is a comparison between simulation and experimental data of the drainage phenomenon in the third embodiment of the present invention.

[0051] In the drawings, like reference numerals are used for like parts. The drawings are not necessarily to scale. DETAILED DESCRIPTION

[0052] The present invention will be further described below in conjunction with the accompanying drawings.

[0053] like Figure 1 , which is a flow chart of a method for high-precision wetting boundaries based on the color lattice Boltzmann method in the present invention. The method for high-precision wetting boundaries based on the color lattice Boltzmann method in the present invention comprises the following steps:

[0054] S1. Obtain the unit normal vector of the solid boundary;

[0055] S2. Obtaining a target phase interface unit normal vector according to the solid boundary unit normal vector and a set contact angle;

[0056] S3. Introducing the target phase interface unit normal vector into the formula for calculating the interface force, introducing the zero interface force condition, realizing the calculation of the interface force at the boundary, and completing the introduction of the boundary wetting condition.

[0057] Embodiment 1

[0058] In this embodiment, the problem of a stationary droplet adhering to a horizontal solid wall surface is taken as an example to specifically introduce the implementation steps of the present invention:

[0059] S1. Obtain the unit normal vector of the solid boundary;

[0060] Specifically, for the binary sawtooth wall, its solid boundary unit normal vector can be given by the formula in step (1):

[0061]

[0062] Among them, n s (x) is the unit normal vector of the solid boundary at the x position,

[0063] c k is the direction of the eighth-order isotropic difference,

[0064] w(|c k | 2 ) is the corresponding weight coefficient,

[0065] δ t is the time step,

[0066] s(x) is the attribute indicator function, and its specific value is shown in formula (2):

[0067]

[0068] For horizontal walls, the solid boundary unit normal vector of the solid wall is the same everywhere, which is n s(x) = (0, -1).

[0069] S2. Obtaining a target phase interface unit normal vector according to the solid boundary unit normal vector and a set contact angle;

[0070] First, according to a given contact angle θ, the solid boundary unit normal vector is rotated to obtain a phase interface unit normal vector, and the angle between the phase interface unit normal vector and the solid boundary unit normal vector is θ; there are two phase interface unit normal vectors, namely n1 and n2,

[0071] n1=(n s，x cosθ-n s，y sinθ,n s，y cosθ+n s，x sinθ)

[0072] n2=(n s，x cosθ+n s，y sinθ,n s，y cosθ-n s，x sinθ)

[0073] Among them, n s，x is the projection size of the unit normal vector of the solid boundary in the x direction of the coordinate axis,

[0074] n s，y is the projection size of the unit normal vector of the solid boundary in the y direction of the coordinate axis,

[0075] θ is the contact angle.

[0076] Then, according to the Euler distance between the two phase interface unit normal vectors and the predicted normal unit vector, the target phase interface unit normal vector is determined, as shown in equations (3) and (4).

[0077] D1=|n1-n * |, D2 = |n2-n * | (3)

[0078]

[0079] Among them, n * is the predicted unit normal vector,

[0080] n f is the unit normal vector of the target phase interface,

[0081] n1 and n2 are two unit normal vectors of the phase interface,

[0082] D1 and D2 are the Euler distances between the unit normal vectors of the two phase interfaces and the predicted unit normal vector.

[0083] Preferably, the predicted unit normal vector n * According to formula (5),

[0084]

[0085] in, To predict the color gradient.

[0086] Preferably, the predicted color gradient is obtained according to the following method: first, a virtual fluid density distribution is introduced into the outermost solid node by weighted averaging; then, the predicted color gradient is calculated from the density distribution of the solid node and its adjacent points.

[0087] In this implementation, we first need to calculate two phase interface unit normal vectors n1 and n2 to introduce the contact angle θ. Figure 2 Taking θ=60° as an example, we calculate It can be seen that for the two contact line positions on the left and right, their target interface unit normal vector is one of the two interface unit normal vectors n1 and n2.

[0088] In order to identify which of the two phase interface unit methods is the target phase interface unit normal vector, we need to introduce a predicted unit normal vector. Taking the simulation of the working conditions with a calculation area size of 200*100, an initial droplet radius of 40, a surface tension coefficient of 0.02, and two fluid viscosities of 0.35mPa*s as an example, we select the time step t=500 and the point (59, 2) in the three-phase contact line area as the observation object. First, the virtual fluid density distribution is introduced into the solid nodes close to the fluid by weighted averaging, and then the color gradient can be calculated from the density distribution of this point and its nearby points: Therefore n * =(0.9711, -0.2386), and then calculate the two Euler distances to get D1 = 1.8556 and D2 = 0.2817. Therefore, the unit normal vector of the target interface is

[0089] S3. Introducing the target phase interface unit normal vector into the formula for calculating the interface force, introducing the zero interface force condition, realizing the calculation of the interface force at the boundary, and completing the introduction of the boundary wetting condition.

[0090] Through rigorous mathematical derivation, it is proved that when a given contact angle condition is introduced, the fluid should satisfy a zero interface force condition in the three-phase contact region. This condition can be converted into that the fluid node in the three-phase region and in contact with the solid should satisfy the condition that the gradient of its target unit normal vector It is further proved that the method of introducing weighted average virtual points to calculate the gradient actually implements a zero interface force condition. Therefore, when calculating the gradient of the unit normal vector of the target phase interface of the outermost fluid node, When the contact angle is introduced, the method of introducing weighted average virtual points is directly adopted. The zero interface force condition is implemented when calculating the gradient of the unit normal vector of the target phase interface adjacent to the solid wall after the contact angle is introduced. It is a necessary step in calculating the interface force.

[0091] The method of high-precision wetting boundary based on the color lattice Boltzmann method in the present invention is applied to the problem of simulating a stationary droplet sitting on a horizontal solid wall surface based on the existing color lattice Boltzmann model framework. The present invention is based on the existing color lattice Boltzmann model framework, and the boundary conditions are improved. The color lattice Boltzmann model is already known. In the embodiments, there is no need to elaborate on the calculation of other parts of the color lattice Boltzmann model. It is believed that those skilled in the art can apply the method of the present invention thereto.

[0092] In this embodiment, five contact angles are set, namely 30°, 60°, 90°, 120° and 150°, and the viscosity ratio is calculated to be 1 and 100. The simulated contact angles are shown in Table 1.

[0093] Table 1 Simulation results of different contact angles of a droplet sitting still on a horizontal wall

[0094]

[0095] It can be found from Table 1 that the simulation results are in good agreement with the theoretical results. Figure 3 The results of the patent model), the original geometric formula, and the virtual density formula were compared for non-physical mass migration, and five contact angles were simulated.

[0096] from Figure 3 It can be seen that the present invention (solid line in the figure) effectively reduces the non-physical mass migration effect.

[0097] Embodiment 2

[0098] In this embodiment, the method of the present invention is used to simulate the capillary permeation problem, and the initialization modeling of the problem is as follows: Figure 4 As shown, where L is the capillary length and W is the capillary width.

[0099] In this embodiment, two working conditions with viscosity ratios of 1 and 100 were calculated, and four contact angles were discussed, namely 10°, 30°, 45° and 60°. This is because when the contact angle is 90°, the surface tension reaches equilibrium and capillary force cannot be caused on the interface; when the contact angle is greater than 90°, the capillary force is similar to the case when the contact angle is less than 90°.

[0100] Figure 5 is a comparison diagram of the theoretical solution of the capillary penetration problem under different wettability and viscosity ratios in Example 2 of the present invention, the solution of this model and the solution of the Lecalire model. Figure 5 It can be seen that the method of high-precision wetting boundary based on the color lattice Boltzmann method in the present invention, the solution obtained based on the existing color lattice Boltzmann model framework (hereinafter referred to as the present model solution) can accurately match the theoretical solution, and has higher accuracy than the Leclaire model solution obtained by the algorithm of Leclaire et al., especially when the viscosity ratio is high and the contact angle is close to 90°. This shows that the present invention has higher accuracy.

[0101] In summary, the method of the present invention can accurately calculate the contact angle in a very wide range, for example, the contact angle in the range of 10° to 170° can be calculated), and the more extreme contact angle will be limited by the number of grids. At the same time, this method can calculate the problems of curved surface boundaries and zigzag boundaries with high accuracy, and effectively reduce non-physical mass migration, especially for high viscosity specific dynamic problems, the method of the present invention also shows high accuracy.

[0102] Embodiment 3

[0103] The present invention is used to simulate the suction process and the displacement process in the pore network micromodel. The geometry used in the simulation is consistent with the microchannel chip used in the experiment by Dr. Wu Rui et al. (see Wu R, Kharaghani A, Tsotsas E. Two-phase flow with capillary valve effect in porous media. Chemical Engineering Science, 2016, 139: 241-248). The numerical simulation is performed on a calculation area of ​​900×1000, and the resolution of the grid is 10 microns. Compared with the microchannel used in the experiment, the area of ​​the numerical simulation is extended by 100 grid sizes at the outlet to reduce the influence of the outlet boundary conditions on the numerical simulation results. The displacement fluid is injected from the inlet channel with a width of 500 microns at the bottom of the microchannel chip, and the viscosity ratio of the displacement fluid and the driven fluid in the simulation is consistent with the viscosity ratio in the actual simulation. The viscosity of the displacement fluid and the driven fluid is 1 / 6 and 0.003 grid units respectively. According to the above parameters, the capillary number Ca in this simulation is calculated to be 10 -4 , and the capillary number 10 in Dr. Wu Rui's experiment -8 However, since the viscous pressure drop at this capillary number is also negligible, further reducing the capillary number has almost no effect on the simulation results. Figure 6a and Figure 6b As shown, Figure 6a This is a comparison chart of the simulation results and experimental results during the inhalation process. Figure 6b The figure is a comparison of the simulation results and experimental results in the displacement process, where the displacement fluid, displaced fluid and solid are represented by black, white and gray respectively. By comparing the experimental and simulation results, we found that the two-phase fluid distribution and interface morphology obtained by simulation and experiment are almost the same, which verifies the effectiveness and high accuracy of the algorithm of the present invention in simulating the two-phase displacement process at the pore scale.

[0104] Although the present invention has been described with reference to preferred embodiments, various modifications may be made thereto and parts thereof may be replaced by equivalents without departing from the scope of the present invention. In particular, the various technical features mentioned in the various embodiments may be combined in any manner as long as there are no structural conflicts. The present invention is not limited to the specific embodiments disclosed herein, but includes all technical solutions falling within the scope of the claims.

Claims

1. A method for high-precision wetting boundaries based on the color lattice Boltzmann method, comprising the following steps: Step 1, obtain the unit normal vector of the solid boundary; Step 2: obtaining a target phase interface unit normal vector according to the solid boundary unit normal vector and a set contact angle; Step 3: Introduce the unit normal vector of the target phase interface into the formula for calculating the interface force, introduce the zero interface force condition, realize the calculation of the interface force at the boundary, and complete the introduction of the boundary wetting condition; In step 1, the unit normal vector of the solid boundary is obtained according to formula (1): Among them, n s (x) is the unit normal vector of the solid boundary at the x position, c k is the direction of the eighth-order isotropic difference, w(|c k |2) is the corresponding weight coefficient, δ t is the time step, s(x) is the attribute indicator function, and its specific value is shown in formula (2):

2. The method according to claim 1, characterized in that Step 2 includes the following steps: First, according to a given contact angle θ, the solid boundary unit normal vector is rotated to obtain a phase interface unit normal vector, and the angle between the phase interface unit normal vector and the solid boundary unit normal vector is θ; there are two phase interface unit normal vectors, namely n1 and n2, n1=(n s,x cosθ-n s,y sinθ,n s,y cosθ+n s,x sinθ) n2=(n s,x cosθ+n s,y sinθ,n s,y cosθ-n s,x sinθ) Among them, n s,x is the projection size of the unit normal vector of the solid boundary in the x direction of the coordinate axis, n s,y is the projection size of the unit normal vector of the solid boundary in the y direction of the coordinate axis, θ is the contact angle; Then, the target interface unit normal vector is determined according to the Euler distance between the two interface unit normal vectors and the predicted unit normal vector, as shown in equations (3) and (4). <h2 style=";text-align:left;direction:ltr">D1 = |n1-n<h2 style=";text-align:left;direction:ltr"> * <h2 style=";text-align:left;direction:ltr"> |,D2 = |n2-n<h2 style=";text-align:left;direction:ltr"> * <h2 style=";text-align:left;direction:ltr"> | (3) Among them, n * is the predicted unit normal vector, n f is the unit normal vector of the target phase interface, n1 and n2 are two unit normal vectors of the phase interface, D1 and D2 are the Euler distances between the unit normal vectors of the two phase interfaces and the predicted unit normal vector.

3. The method according to claim 2, characterized in that The predicted unit normal vector n * According to formula (5), in, To predict the color gradient.

4. The method according to claim 3, characterized in that The predicted color gradient is obtained according to the following method: Firstly, the virtual fluid density distribution is introduced into the outermost solid nodes by weighted averaging; Then, the predicted color gradient is calculated from the density distribution of the solid node and its neighboring points.

5. The method according to any one of claims 1 to 4, characterized in that: In step three, the zero interface force condition is introduced by introducing a weighted average virtual point method.

6. The method according to claim 5, characterized in that The zero interface force condition is the gradient of the unit normal vector of the target phase interface of the fluid node at the boundary

Images