Lost circulation control in fractured reservoirs is a hot and difficult problem in the field of drilling engineering. Lost circulation not only consumes a large amount of working fluid and lost circulation materials, and causes huge economic losses directly, but also increase the non-production time, extends the drilling cycle, affects the exploration and development process, and even causes sticking, borehole collapse, blowout and other accidents. The drilling fluid loss in reservoir sections would seriously impede the timely discovery of oil and gas and cause significant reduction of oil and gas production^{[1, 2]}. The high temperature, high pressure and high formation stress of reservoirs in deep- ultra deep wells further add difficulties to lost circulation control and formation damage prevention^[3,4,5].

In order to meet the requirements of controlling working fluid loss and preventing formation damage of naturally fractured reservoir, domestic and foreign researchers have done a lot of studies by means of laboratory experiments, theoretical modeling, and numerical simulation. In line with induced fracture type, extended fracture type and large and medium fracture type fluid loss, methods such as regulating stress around the well, plugging the loss channel and improving the rock mass strength have been put forward^[6,7,8,9]. The method of regulating stress around well prevents formation fracture and thus avoid further opening of natural non-leaking fractures by improving circumferential tangential stress, fracture closure stress and fracture extension pressure. This method is more suitable for the loss of induced fracture type and fracture expanding type. The design of the particle size of loss control materials (LCMs) emphasizes the formation efficiency of the plugging zone. The method of plugging loss channel relies on dense and high-strength plugging zone formed by physical or chemical materials to seal the fracture, and is mainly adaptable for the fluid loss of fracture expansion type and large and medium-sized fracture type. The design of the particle size of LCMs emphasizes more on the strength and stability against pressure of the sealing layer. The method of improving the strength of rock mass uses chemical materials to form high-strength structure and is suitable for fluid loss of large and medium-sized fracture type. The effective implementation of these lost circulation control methods all relies on the fracture plugging zone to establish the balance between wellbore liquid column pressure and ground stress field and formation pressure field. Using the physical materials to form plugging zone that can be removed is the most common way for plugging the naturally fractured reservoirs. To meet the requirements of reservoir protection and fluid loss control at the same time, the plugging materials shall be able to form bridge rapidly, and the plugging zone formed must be tight and able to withstand high pressure^[10,11,12]. Otherwise, excessive loss before plugging zone formation or repetitive loss caused by low strength plugging zone would aggravate the damage degree and enlarge the damage scope, resulting in difficulty to remove the fluid loss damage to reservoir^[13,14]. The plugging efficiency and strength of fractures are the key factors affecting the effects of working fluid loss control and reservoir protection in fractured reservoirs. The formation and evolution mechanism of fracture plugging zone structure determine the efficiency and strength of fracture plugging, but there is no systematic research on the formation and evolution mechanisms of fracture plugging structure.

Coupled CFD-DEM (computational fluid dynamics - finite element method) is adopted to simulate the formation and evolution of fracture plugging zone structure. A photo-elastic system using mesoscale force chain network to characterize plugging zone is developed to simulate the formation and evolution of fracture plugging zone structure, reveal the formation and evolution of the plugging zone geometric and mechanic structure. A new method of material selection and formula design for plugging is proposed. A new type of high efficient and strong retention plugging materials are designed and selected, and the plugging formula is optimized. The research results provide theoretical and technical basis for improving the control effect of formation damage and lost circulation of fractured reservoir.

The CFD-DEM method coupling fluid mechanics and particle discrete element is an effective way to simulate the formation process of fracture plugging zone structure. According to the geometric characteristics of fractures in deep fractured reservoirs of the Tarim Basin, wedge- shaped fracture was selected as the geometric morphology of the fracture model. The inlet and outlet width of the fracture were set at 3 mm and 1 mm respectively, and the fracture length was set at 50 mm, close to the size of fracture modules commonly used in laboratory experiments. The fracture geometric model is shown in Fig. 1.

In the CFD-DEM simulation of this work, the relevant parameters of the LCMs were set according to calcite type LCM most commonly used in fractured reservoirs. The friction coefficient of the plugging material was measured by THE COF-1^[5] measuring device. The fluid-related parameters were set according to the parameters of deep well drilling fluid and the drilling fluid loss rate of fractured formation in Tarim Basin, NW China. The basic parameters of the plugging materials and drilling fluid are shown in Table 1.

Coupling CFD-DEM simulation is divided into resolved method and unresolved method. The particle size of resolved method was obviously larger than that of fluid computing unit (Fig. 2). The dynamic mesh was used to refine the local mesh (Fig. 2), to simulate the flow field and particle force around the particle in more detail and calculate the particle motion behavior^[15]. In studying formation mechanism of the plugging zone structure, special attention should be paid to the microscopic bridging and sealing behavior of the LCMs. The resolved method was adopted to carry out the simulation. The computational domain of resolved CFD-DEM method was composed of fluid domain and particle domain^[16].

For the fluid domain part, the Navier-Stokes equations, boundary conditions, initial conditions and fluid-particle coupling conditions of the resolved CFD-DEM simulation method are as follows^[17].

Fluid continuity equation:

Fluid motion equation:

Boundary conditions:

Initial conditions:

Fluid-particle interface coupling conditions:

For the particle domain, the trajectory of each particle was calculated under the condition of considering the interactions between the particle and surrounding other particles and solid boundaries. The particle motion includes translation and rotation, and the translation acceleration and angular acceleration of the particle were calculated based on the corresponding momentum balance. According to Newton's second law, the translation governing equation of particle i is

The torque of particle j acting on particle i consists of tangential moment and rolling friction moment. According to Euler's second law of motion, the rotation governing equation of particle i with moment of inertia I_i is expressed as:

In the process of coupling calculation, the fluid computing domain occupied by particles requires detailed grid resolution. In order to obtain more accurate simulation results, the ratio of characteristic length of fluid computing unit to particle diameter should not exceed 1/10^[18]. In order to improve the computational efficiency, dynamic grids were used to refine the grids occupied by particles. After the particles moved, the previously occupied meshes were coarsened again. The coupling of CFD and DEM was accomplished by mass, momentum and energy exchange between the fluid and particles. The specific process of CFD-DEM coupling is shown in Fig. 3.

In order to verify the reliability of the CFD-DEM coupling model, a series of particle settlement issues were simulated by using the established simulation method. The size of the settlement area was 25d_p×25d_p×250d_p. The initial position of the sphere was located at the geometric center of the top surface, with a height of about 240d_p. Particles of ten different sizes in total were simulated. The numerical simulation results of particle settlement were compared with the calculation results of particle settlement velocity equation derived by Concha^[19] and the experimental data of Lapple and Shepherd^[20], which verified the reliability of the simulation method presented in this paper (Fig. 4).

1.3.1. Effects of size and volume concentration of plugging material on fracture plugging efficiency

Bridging of plugging materials is the key to the plugging of fractures. After the bridging occurs, subsequent LCMs gradually fill and accumulate to form fracture plugging zone. The volume concentration and particle size distribution of LCMs are important factors determining bridging and formula design. Fig. 5 shows that at a given concentration of plugging material, the ratio of the material diameter to fracture width (R) directly determines the plugging efficiency. At the plugging material concentration of 5% (volume fraction) and friction coefficient of 0.8, at the R value of 0.5, 0.6 and 0.7, simulations showed particles constantly flew out of the crack outlet during the plugging process, and no effective plugging was formed inside the fracture. As R value increased to 0.8, after a period of time, no particles flew out of the fracture outlet, and a plugging zone was formed inside the fracture. The larger the R value, the shorter the plugging time and the higher the plugging efficiency was.

The concentration of LCMs is an important parameter in the design of plugging formula, which was represented by volume fraction in the simulation, that is, the ratio of the volume of plugging material to the volume of plugging slurry. The concentration mainly affects the interaction between particles. With the increase of plugging material concentration, the interaction between particles becomes stronger and the possibility of bridging increases. The bridging behavior of particles is a probabilistic event, which is dependent on the ratio of particle size to fracture width and concentration of the plugging material^[21,22]. The efficiency of fracture plugging is reflected in the time required for particles to form bridging, which is essentially dependent on the probability of bridging, namely bridging probability. Fig. 6 shows that at the R value of 0.7 and friction coefficient of 0.8, when the material concentration is 5%, particles constantly flow out of the fracture outlet, and bridging cannot occur. When the material concentration increases to 10% or more, the interaction between particles enhances and the probability of bridging increases. The larger the plugging material concentration, the shorter the time needed for forming plugging zone and the higher the plugging efficiency are.

1.3.2. Critical and absolute bridging concentration

The bridging behavior of the LCMs in the fracture is random, and the bridging probability is used to characterize the difficulty of bridging in the fracture during the plugging process. Based on CFD-DEM simulation method, bridging probability is defined as:

Based on the relationship between bridging probability and material concentration, the concepts of critical and absolute bridging concentration are proposed (Fig. 7). When the material concentration is low, bridging does not occur, in other words, the bridging probability is zero. When the material concentration increases to a certain critical value, there is a certain probability of bridging, and this material concentration is the critical bridging concentration. As the material concentration increases further, the bridging probability gradually increases. When the bridging probability reaches 100%, the corresponding bridging concentration is the absolute bridging concentration. When the material concentration is higher than the absolute bridging concentration, bridging is bound to occur in the fracture, that is, the bridging probability is 100%.

R value is the key factor affecting the critical and absolute bridging concentrations. Fig. 8 shows that at the friction coefficient of 0.8, both the critical and absolute bridging concentrations tend to decrease with the increase of R value. Based on the critical and absolute bridging concentrations, Fig. 8 can be divided into three regions. When the material concentration and R are located in the red region, P_b=0, bridging does not occur. When the material concentration and R are located in the yellow area, 0<P_b<100%, there is a certain probability of bridging. When the material concentration and R are located in the blue region, P_b =100%, bridging is bound to occur. The critical and absolute bridging concentrations provide a theoretical basis for the design of bridge material concentration. Considering the efficiency of fracture plugging and the cost of LCMs, the lower limit of bridging material concentration should be larger than the critical bridging concentration, and the upper limit should be slightly larger than the absolute bridging concentration (in consideration of additional consumption of bridging materials, such as particle size degradation).

1.3.3. Relationships between shape and particle size and absolute bridging concentration of plugging material

Shape is an important parameter to characterize the shape of material. Shape affects the friction coefficient of the material and the plugging efficiency of the fracture^[23]. Xie et al.^[24] analyzed the bridging behavior of non-spherical particles in narrow passages by changing the friction coefficient of materials. Fig. 9 shows that at different R values, the absolute bridging concentration decreases with the increase of friction coefficient of the plugging material. When R value is greater than 0.8, friction coefficient has little influence on absolute bridging concentration and material bridging behavior. With the decrease of R, the influence of friction coefficient becomes more and more significant. By increasing the irregularity degree of material and sliding/rolling friction coefficients, the absolute bridging concentration can be effectively reduced and the fracture plugging efficiency can be enhanced.

1.3.4. Formation of fracture plugging zone structure

Bridging is the first step in the formation of fracture plugging zone structure. After bridging is formed, materials will further accumulate and fill on the basis of bridging, and the plugging zone will gradually improve in compaction and strength. By analyzing the influences of material concentration, R value, shape/friction coefficient and other factors, this study reveals the mechanism of material bridging in wedge- shaped fracture (Fig. 10). When R≥1, single-particle bridging occurs. When 0.7≤R<1, two particles at the front first bridge within the fracture at an angle, and the subsequent particles are further retained to form a stable structure, that is, double-particle sequential bridging. The occurrence of single particle bridging and double-particle sequential bridging is mainly governed by particle size. As particle size decreases further, when 0.5≤R<0.7, parallel two particle bridging and multi-particle bridging are dominant. When R is less than 0.5, only multi-particle bridging occurs in the fracture. Parallel two particle bridging and multi-particle bridging are mainly governed by the combination of material diameter and friction force. With the decrease of R and the increase of absolute bridging concentration, the governing factor of bridging changes from material diameter to the combination of material diameter and friction force. In the design of plugging formula, the R values that can result in single particle bridging and double particle sequential bridging are selected preferentially. However, when dealing with fractures millimeters to centimeters wide, too large bridging particles often make the plugging slurry difficult to be pumped or downhole tools unable to be used. In this case, the R value that satisfies the pumping requirements of the slurry and the normal use of the downhole tool is determined first. Then, according to the R value, the absolute bridging concentration is determined, and the concentration or the friction coefficient of the plugging material is increased to realize high efficient plugging.

After entering into the fracture, the LCM particles migrate, bridge, pile and finally form the plugging zone. The formation speed of the plugging zone geometric structure is dependent on the bridging efficiency. Whereas the strength and pressure evolution mechanism of the plugging zone after formation is dependent on the mechanical structure, namely the mesoscale force chain network (Fig. 11). In this study, photo-elastic system accounting for mesoscale force chain network characterization was developed to analyze the evolution process of fracture plugging zone.

The fracture plugging zone has a multi-scale structure. During the plugging of fracture, the plugging materials contact each other and form a contact force network, namely a mesoscopic force chain network, which constitutes the mesoscopic scale of the plugging zone^[25]. The stability of the mesoscopic force chain network is affected by the properties of the plugging material at the microscopic scale and determines the strength of the macroscopic plugging zone. The mesoscopic force chain network is the best entry point to analyze the pressure evolution mechanism of the structure of fracture plugging zone.

The photo-elastic experiment reflects the distribution of contact force in particle system by optical interference principle. In the field of particle mechanics, rock-soil mechanics, rock mechanics and so on, it is used to study the stress distribution in the material system. It is an effective method to characterize the mesoscopic force chain network of fracture plugging zone^[26]. Shear instability is the most common structural instability mode of fracture plugging zone^[27].

Based on the principle of photo-elastic experiment, a photo- elastic experiment system was developed to characterize the evolution behavior of mesoscopic force chain network of the fracture plugging zone under vertical and shear loads (Fig. 12). The fracture surfaces of the system are transparent perspex plates. The fracture width is 5 mm, the seam size is 260 mm×260 mm, the horizontal shear band width is 20 mm, and the maximum shear depth is 10 mm. The system can apply vertical load and horizontal shear load, and can simulate the shear instability mode of plugging zone. The experimental material is polycarbonate, which has the advantages of high optical sensitivity and transparency and small creep at room temperatures^[28].

2.2.1. The pressure evolution mechanism of fracture plugging zone structure

Fig. 13 shows the pressure variation in the shear process of the fracture plugging zone. Four points A, B, C and D were selected for analysis. The optical elastic images of the mesoscopic force chain network corresponding to the four pressure points were shown in Fig. 14. The brighter part in Fig. 14 indicates the greater the contact force, and it can be seen that the structure change of the force chain network is the mesoscopic response of the change of exerted force on the macroscopic plugging zone. The contact force between two particles can be calculated based on the average square gray gradient of the pixels of the photo-elastic image of the particles^[29]. The force chain with contact force greater than the average is defined as the strong force chain. The photo-elastic experiment results show that the proportion of strong force chain is positively correlated with the stress exerted on the plugging zone (Fig. 13). The proportion of strong chains is an important index to evaluate the structural stability of mesoscopic force chain network.

2.2.2. Optimization of mesoscopic fracture plugging zone structure

The proportion of strong chain and maximum pressure bearing capacity of the plugging zone in the mesoscopic force chain network corresponding to different material shapes, friction and fluid environments, material types and porosities of the plugging zone during the shearing process are further compared (Table 2). The results show that higher proportion of strong chains under different parameters corresponds to higher strength of the fracture plugging zone. By reducing the roundness and increasing the friction coefficient of the material, and adding elastic material to enhance the tightness of the plugging zone, the structure of mesoscopic force chain and the strength and stability of the fracture plugging zone can be effectively improved. A new method to evaluate LCMs has been established by using the photo-elastic experimental system, and taking the ratio of strong chain as the evaluation index.

By combining the CFD-DEM simulation and photo- elastic experiment method proposed in this paper with the systematic property evaluation method of LCMs^[30], a new method of selecting plugging materials and optimizing plugging formula has been worked out, which can effectively guide the design and selection of plugging materials and optimization of plugging formula. According to the formation and evolution mechanism of fracture plugging zone structure, LCM-K3, a new type of high efficient plugging strong retention LCM, was selected. This material had low spherical degree, high friction coefficient (1.45), and high temperature and pressure resistance capacity (with D₉₀ degradation rate of less than 10% at 200 °C and 30 MPa). The strength of fracture plugging zone formed by the selected LCM was tested with the high-temperature high-pressure (HTHP) fracture plugging apparatus. The fracture module had an inlet 8 mm wide and an outlet 5 mm wide.

The experimental results show that compared with the calcite bridging material LCM-D4, the plugging zone of LCM-K3 has much higher maximum pressure bearing capacity, which effectively improves the fracture plugging strength (Table 3, Fig. 15). It should be noted that in Table 3, the concentration in each formula is the additive amount expressed by the ratio of the material mass to the volume of the drilling fluid in the unit of g/mL. For example, 5% represents 0.05 g/mL.

On the basis of selecting the new type of LCMs, the bridging material concentration was further optimized according to the absolute bridging concentration to improve the fracture plugging efficiency. Simulation was carried out to optimize the formula 1-1# based on the material and fracture parameters, and the C_av was optimized at 3.1%. According to formula (9), the C_av was converted to the corresponding C_am of 0.052 g/mL, which was represented by 5.2% in the formula. Gradient concentration formulas 2-0# - 2-3# were set. The experimental results show that when the amount of bridging material in the formula is too low, plugging zone cannot be formed efficiently, and the slurry would lose continuously; when the material concentration is 0.055 g/mL, that is, higher than C_am, the cumulative loss significantly reduces and the fracture plugging efficiency enhances considerably (Table 3, Fig. 16).

Particle bridging is the key to the formation of fracture plugging zone structure, and bridging probability is the essential factor affecting fracture plugging efficiency. Critical and absolute bridging concentrations can be taken as the key basis for loss control formula design. With the increase of absolute bridging concentration, the governing factor of bridging changes from material diameter to the combination of material diameter and friction force.

Mesoscale force chain network is the intrinsic factor that affects the evolution of fracture plugging under external force and determines the macroscale strength of the fracture plugging zone. Geometric and mechanic parameters of loss control material affect the stability and strength of fracture plugging zone through affecting the ratio of strong force chain.

Results of fracture plugging experiments show that selecting LCMs with higher efficiency and stronger retention and optimizing loss control formula according to results of this study can enhance the fracture plugging efficiency and strength of the plugging zone, and improve the control of working fluid loss in deep fractured reservoirs.

C_am—absolute bridging concentration expressed by the ratio of material mass to base mud volume, g/mL;

C_av—absolute bridging concentration expressed by the ratio of material volume to plugging mud volume, %;

d_p—particle diameter, m;

D₉₀—particle size corresponding to the cumulative frequency of 90% on the cumulative particle size distribution curve, m;

f_pf,i—force vector between particle i and fluid, N;

f_pp,ij—force vector between particle i and particle j, N;

g—acceleration vector of gravity, m/s²;

i—particle number;

I_i—moment of inertia of particle i, kg·m²;

j—number of all particles interacting with particle I;

m_i—mass of particle i, kg;

M_pf,i—moment of fluid acting on particle i, N·m;

M_r,ij—rolling friction moment of particle j acting on particle i, N·m;

M_t,ij—tangential moment of particle j acting on particle i, N·m;

n—total number of particles interacting with particle i;

$\hat{\boldsymbol{n}}$—external normal vector of particle;

N—total number of cases (events);

N_b—number of cases (events) of bridging;

p—pressure, Pa;

P_b—bridging probability, %;

R—ratio of particle size to fracture width;

t—time, s;

t_p—traction vector of fluid acting on particle surface;

u₀—fluid initial velocity vector, m/s;

u_f—fluid velocity vector, m/s;

u_p—particle velocity vector, m/s;

u_b—boundary velocity vector, m/s;

v_i—translation velocity vector, m/s;

x—position coordinates, m;

μ_f—fluid viscosity, Pa·s;

ρ_f—fluid density, kg/m³;

ρ_p—bridging material density, g/mL;

σ—stress tensor in fluid;

ω_i—angular velocity of particle i, rad/s.