In the geologic model for fractured horizontal wells within shale oil reservoirs, a binary porosity medium is proposed, comprised of two parts: the matrix substance and the fracture web. The matrix substance functions as the primary storage space, while the fracture network operates as the main channel for shale oil. Both oil and water have the capacity to move within the fractures, however, only oil has the ability to flow within the matrix. The dual-porosity model encompasses both the matrix and fracture systems: the former is the chief repository for shale oil, and the latter is the primary passage for it. Three basic assumptions are taken into account:

(1) The flow of oil-water takes place in the fracture system, while solely single-phase shale oil migrates towards fractures in the matrix, disregarding the flow within the matrix system;

The meshless method, which is based on the weighted residual technique, utilizes the moving least squares (MLS) strategy to construct the approximation function. The governing equation is discretized using the least square method, resulting in the development of MWLS, a meshless technique. For a comprehensive understanding of the MWLS, please refer to the cited literature. (Yu-kun, 2007).

By reorganizing Eqs 11–13, we can derive the subsequent Eq. 18. ∇ λ f g t ∇ p f t + 1 + ∇ λ f w t ∇ p f t + 1 + q m g t − ϕ f C t f t p f t + 1 − p f t ∂ t = 0 ① ∇ λ f w t ∇ p f t + 1 − ϕ f t S w t + 1 − S w t ∂ t = 0 ② σ k m μ g p m t + 1 − p f t + 1 + ϕ m C t m p m t + 1 − p m t ∂ t = 0 ③ (18)

In the process of executing MWLS for inferential purposes, the problem domain, represented by Ω, and its boundary, denoted by Г, undergo discretization through the use of n points. The fracture system’s pressure and water saturation are then modeled by meshless approximate functions, established courtesy of MLS, which can be expressed following the template of Eq. 19. p f ≈ ∑ i = 1 n N i P f i = N d f (19) where N = N 1 x   N 2 x N 3 x ⋯ N n x is a shape function, d f = p f 1   p f 2   p f 3 ⋯ p f n T .

By replacing Eq. 18-①) with Eq. 19 from above, the residual can be derived, as depicted in Eq. 20. R p f = ∇ λ f w t ∇ N d f t + 1 + ∇ λ f g t ∇ N d f t + 1 + σ k m μ g p m t − N d f t + 1 − ϕ f C t f Δ t t N d f t + 1 − N d f t = 0 (20)

Thus, the computational scheme of MWLS for determining the pressure in the fracture system unfolds as follows: Π p f = ∫ Ω R p f T R p f d Ω + ∫ Γ 1 α 1 N d f t + 1 − p f c N d f t + 1 − p f c d Γ + ∫ Γ 2 α 2 ( ∂ p f t + 1 ∂ x ) ∂ p f t + 1 ∂ x d Γ (21)

In Eq. 21, α 1 and α 2 symbolize the penalty functions for the boundary conditions. By identifying the least value of Eq. 21 and transposing it into a discrete format, we arrive at Eq. 22: δ Π p f = ∑ s = 1 n δ R p f R p f | x = x s + ∑ s = 1 n 1 α 1 δ p f t + 1 N T N d f t + 1 − p f c | x = x s + ∑ s = 1 n 2 α 2 δ p f t + 1 ∇ N T ∇ N | x = x s = 0 (22)

By transforming Eq. 22 into matrix system and taking into account the freedom of p f t + 1 , we derive Eq. 23. K 1 p f t + 1 = F 1 (23)

Among them: K 1 = ∑ s = 1 n δ R p f ∇ λ f t ∇ ∇ N − σ k m μ g N − ϕ f C t f Δ t t N + ∑ s = 1 n 1 α 1 N T N + ∑ s = 1 n 1 α 2 ∇ N T ∇ N F 1 = ∑ s = 1 n δ R p f − σ k m μ g p g m t − ϕ f C t f Δ t t d f t + ∑ s = 1 n 1 N T p f c

Eq. 24 provides the fracture’s pressure at the t+1 time step. This value is then utilized in Eq. (18-②) for substitution. Subsequently, the MLS method is employed to construct an approximate equation for the Sw in the fracture. Ultimately, Eq. 24 is derived: K 2 S w t + 1 = F 2 (24)

Eq. 24 provides the fracture’s Sw at the next time step.

The variables and parameters involved are as follows: K 2 = ∑ s = 1 n N T N + ∑ s = 1 n 1 α 1 N T N + ∑ s = 1 n 2 α 2 ∇ N T ∇ N F 2 = ∑ s = 1 n N T e w f t + Δ t ϕ f ∇ λ w t ∇ p f t + 1 + ∑ s = 1 n 1 α 1 N T S w f c t

The matrix system’s pressure can be determined using a similar approach. By inserting the fracture system’s pressure at the next time step into formula (18-③), the final MWLS equation can be derived: K 3 p m t + 1 = F 3 (25)

Eq. 25 provides the matrix’s pressure at the next time step.

The variables and parameters involved are as follows: K 3 = ∑ s = 1 n δ R p m σ k m μ g N + ϕ m C t m N d t + ∑ s = 1 n 1 α 1 N T N + ∑ s = 1 n 2 α 2 ∇ N T ∇ N F 3 = ∑ s = 1 n δ R p m σ k m μ g d f t + 1 + ϕ m C t m d t p m t + ∑ s = 1 n 1 α 1 N T p m Γ 1 t

Where s represents the number of points used for calculations; Let n represent t the total number of points used for calculations within the neighborhood Ω, and let ni (i = 1,2) represent the total number of points used for calculations on boundary Γi.

A field case was constructed to validate the model’s accuracy, considering the fractal features in a shale oil reservoir where the horizontal well is fractured. The reservoir area are 500 m × 400 m × 10m, with the model and layout depicted in Figure 1. The corresponding physical properties of the reservoir can be found in Table 1. Due to its superior continuity and higher computational accuracy, we used the Gaussian formulation as the weighting function in the grid-independent approach. The layout scheme used is 50 × 40, and a node spacing of 3 times (30 m) is used for the node influence domain.

This study incorporated Table 1 data into the oil-water flow fractal shale reservoir model to obtain daily oil production variations. This paper set a result calculated using EDFM as the reference solution. The results were then compared to those of a reference solution, as shown in Figure 2. The proposed method closely matched the reference solution, demonstrating its reliability, accuracy, and high computational efficiency. Furthermore, Figure 3 presents the pressure distribution at different times, showing that the matrix pressure gradually decreases during production. Although there were some discrepancies between the reference solution and the MWLS results, these errors fall within acceptable bounds and can be attributed to inherent differences between the approaches.

After applying the theories and solutions mentioned earlier, a sensitivity analysis of the nodal influence domain is conducted to evaluate its impact on the accuracy of meshless method computations. The solution of the flow of oil and water in horizontally fractured wells. is used to verify the validity of the method. While the MWLS method is inherently precise, the accuracy of its results is also dependent on the process of choosing the nodal influence domain. The optimal parameter values vary depending on the specific problem at hand.

Consequently, the nodal influence domain was chosen to be 2, 3, 4, and 5 times the nodal spacing (10 m), with the initial water saturation is 0.25 while keeping other parameters constant. By doing so, the established mathematical model is further validated and analyzed to ensure its reliability and applicability in solving real-world problems associated with the flow of oil and water in horizontally fractured wells.

Figure 4 demonstrates the variation curves of daily oil and water production for different nodal influence domains. The calculated results demonstrate varying levels of consistency when the nodal influence domain is set at 2, 3, 4, and 5 times the nodal spacing. This suggests that the method is characterized by good stability and convergence in its calculations. This validates its effectiveness and robustness for simulating the flow of oil and water in horizontally fractured wells. Figure 5 shows the crack pressure distribution at 100 days for different nodal influence domains. The calculated crack pressure distributions do not differ significantly under different conditions. When the influence domain size is 5 times the nodal spacing, the pressure drop area is marginally larger than that for the influence domain size of 2, 3, and 4 times the nodal spacing. However, the overall difference is less than 2% because the area of pressure drop increases slightly with the size of the influence domain.

The initial Sw level is a crucial factor affecting shale oil production. To investigate the impact of different initial Sw levels on reservoir recovery and fractured well production while keeping other parameters constant, data from Table 1 were utilized. The initial Sw values were set at 0.35, 0.45, 0.55, and 0.65. Figure 6 illustrates the pressure in the reservoir for different initial levels of Sw after 100 days of production. It is evident that an increase in the initial Sw level leads to a reduction in the pressure difference between the matrix and fracture in the reservoir. This underscores the importance of accurately modeling and understanding initial Sw when simulating and predicting shale oil production. At an initial Sw of 0.35, the matrix pressure in the upper part of the model is about 30 MPa at 100 days, while at this value of 0.65, the matrix pressure in the upper part of the model is only about 26 MPa.

In Figure 7, the Sw distribution in the reservoir at 100 days are shown for different initial Sw levels. When the initial Sw is set at 0.35, the Sw at the fracture increases slightly to 0.3506 after 100 days. However, as the initial Sw level increases, the Sw at the fracture decreases after 100 days. This phenomenon is mainly due to the increase in Sw, which results in a significant rise in water flow capacity. Accurately modeling and predicting the behavior of oil-water flow in shale reservoirs requires an understanding of the relationship between initial Sw and its impact on fracture Sw over time.

In Figure 8, the fluid production under different initial Sw conditions is displayed. As the initial Sw increases, there is less free oil stored in the fractures and reservoirs. Consequently, the initial oil production is significantly reduced when the initial Sw is higher. Moreover, higher Sw reduces the relative permeability of the oil phase, hindering oil flow within the reservoir. Thus, larger initial Sw levels result in lower daily oil production and higher daily water production. For instance, when the initial Sw is set at 0.65, the initial daily well water production rate can reach 1024 m³/day. Understanding these relationships is crucial for optimizing reservoir management and production strategies.

We investigate the effects of different fracture half-lengths on the proposed model by setting the values to 100 m, 150 m, 200 m, 250 m, and 300 m, while keeping other parameters unchanged. This analysis helps to understand how variations in fracture half-length and SRV influence the production performance and recovery efficiency in shale oil reservoirs.

Figure 9 displays the pressure distribution field for different fracture half-lengths at 100 days. Increasing the fracture half-length results in a more rapid and extensive drop in matrix pressure. For instance, with a fracture half-length of 100 m, the matrix pressure in the upper part of the reservoir was maintained at 30 MPa, while at a fracture half-length of 300 m, it decreased to 27 MPa.

Figure 10 depicts the comparison curves of daily and cumulative oil production for various fracture half-lengths. Daily oil production rises with an increase in the fracture half-length, as shown in Figure 10A. Figure 10B indicates a positive correlation between cumulative oil production and fracture half-length, mainly due to the larger stimulated reservoir volume and pressure ripple area as the fracture half-length increases. Despite this, the incremental gain in cumulative oil production diminishes as the fracture half-length increases. For example, when the fracture half-length is 100 m, the cumulative oil production for 100 days is 8075 m³. When the fracture half-length is extended to 150 m, the cumulative oil production increases by 1556 m³ to reach 9631 m³. However, when the fracture half-length is increased from 250 m to 300 m, the increase in cumulative oil production is only 602 m³.

There are multiple factors affecting oil well productivity, among which the number of fractures is one of the key factors. The optimization of fracture numbers is also an important aspect of design of hydraulic fractures for horizontal wells. By fixing fracture half-length at 200 m, setting the initial Sw at 0.65, we change the number of fractures to 2, 3, 4, and 5, keeping other parameters consistent with the previous example.

The pressure distribution at 100 days for different fracture numbers is displayed in Figure 11. The relationship between the number of fractures and the pressure in the reservoir is evident, with an increase in fractures leading to a more rapid decline in pressure and a larger affected area. For example, when there are 2 fractures, the matrix pressure at the upper end of the reservoir can still be maintained at 25 MPa. However, when there are 5 fractures, the matrix pressure at the upper end of the reservoir is only 20 MPa. Figure 12 presents the changes in production capacity under different fracture numbers. Figures 12A, B display the wopr and wwpr curves, respectively, indicating that as the number of fractures increases, the corresponding WOPR and WWPR also increase. Figure 12C shows the FOPT curve, which is also shows a positive correlation with the number of fractures, but the rate of increase is diminishing. Increasing the number of fractures from 2 to 4 results in a 1.28-fold increase in cumulative oil production. However, increasing the number of fractures from 4 to 5 only leads to a 1.05-fold increase in cumulative oil production.