Due to the non-volatility of heavy oil, we consider a three-phase condition (oil phase, liquid water phase, and gaseous water phase). When liquid water and gaseous steam are present simultaneously in the three-phase area, the pressure and temperature should obey Dalton’s law of partial pressure. Assumptions include the following: 1) oil, steam, and water are insoluble in each other; 2) flow in the formation follows Darcy’s law; and 3) the effects of capillary forces and gravity are neglected.

The model’s three primary governing equations (mass conservation in the water phase, mass conservation in the oil phase, and energy equation) are as follows (Ping et al., 2021).

Water mass conservation equation: ∂ ( ρ w ϕ S w + ρ g ϕ S g ) ∂ t + ∇ ( ρ w v w + ρ g v g ) = q w , g (1)

Oil mass conservation equation: ∂ ( ρ o ϕ S o ) ∂ t + ∇ ( ρ o v o ) = q o (2)

Energy conservation equation: ∇ · ( λ ∇ T ) + ∇ · ( ∑ i = o , w , g ρ i e i v i ) = ∂ ∂ t [ ( 1 − ϕ ) ρ i e r + ϕ ∑ i = o , w , g S i ρ i e i ] (3) Where, ϕ , λ , S i , v i , e i , i = o , w , g are rock porosity, heat conductivity coefficient, saturation, flow velocity, internal energy, and enthalpy of gaseous water respectively. q w , g , q o Present the mass source terms of oil and water. T is temperature.

When liquid water and steam are both presented, it should meet: T = T s a p ( P ) (4) With the constraint relationship between saturations: S o + S w + S g = 1 (5)

Along with these two relations in addition to the three main governing equations mentioned above, this is a closed mathematical model.

The meshless method can be theoretically justified using the weighted residual method. The moving least-squares approach is utilized to construct the approximation function in this article, followed by the least-squares method to discretize the control equations and establish the related meshless method, namely the weighted least squares meshless method (Lin et al., 2020).

Let there be a field function u(x) and a set of randomly distributed discrete nodes x _j (j = 1, 2, ..., n) in the solution domain Ω j . Denote by Ω j the influence domain of node x _j . In a two-dimensional problem, the influence domain is often taken as a circle or a rectangle. The global approximation using the moving least squares method has an approximation function for any node x in the solution domain as: u h ( x , x ¯ ) = N ( x , x ¯ ) u = N I u I (6) Where N _I is shape function, x ¯ = [ x , y , z ] T refers to the spatial coordinates of the points in the neighborhood of the calculated point x, u _I is the local best approximation.

The weight function is an important component of the moving least-squares approximation, there are no theoretical specific rules with some arbitrariness, generally satisfying the weight function ω(r) non-negative.

In terms of the variational principle and the weighted residual approach, which is used to build generalized functions and whose weight function or test function is the residual itself, the least-squares method can be used to solve partial differential equations of any type.

The weighted least squares meshless technique is generated by applying the above moving least squares-based approximation function to the specific problem and utilizing the above integral format’s discrete point summation form.

The pressure solution equation is first constructed by combining Eqs. 1, 2, disregarding capillary force and differencing the time term, as explained in Ref. (Desheng et al., 2013). After determining the pressure at the subsequent instant, it is replaced into an equation (Guo et al., 2016). The next moment’s saturation of the seepage field is solved and incorporated into the energy conservation equation (Almao, 2012). Assume that the solution domain and boundary include n discrete nodes, such that W s = ϕ d t ( ρ o S o n + 1 C r o + ρ w S w n + 1 C r w ) + ( 1 − ϕ ) d t ρ r e r (7) Q t = ρ o λ o ∇ 2 p h o + ρ w λ w ∇ 2 p h w + ρ g λ g ∇ 2 p h g (8)

The conservation of energy equation can be transformed into: λ ∇ 2 T n + 1 − W s · T n + 1 + Q t − ϕ Δ t S g e g ρ g + W s · T n = 0 (9)

The residuals of the energy conservation equation are: δ R T = λ R ∇ 2 N T − W s · N T (10)

Then the MWLS form of its energy conservation equation is: k T T n + 1 = f T (11) Where, k T = ∑ s = 1 N δ R T ( λ R ∇ 2 N − W s · N ) + ∑ s = 1 N 1 N T N , f T = ∑ s = 1 N δ R T ( − Q t + ϕ Δ t S g h g ρ g − W s · T n )

Solving Eq. 11 gives the temperature distribution for the next moment.

Successive iterations can be used to overcome the three-phase saturation and temperature problems associated with the steam flooding replacement process.

We demonstrate the strategy described in this research by using a one-dimensional steam flooding example. We created a single-injection model with a model size of 400 m*10 m*10 m and a 5 m model distribution interval. At 4,000 mD, the permeability is uniform, the porosity is 0.3, and the first water saturation is 0.15. On the model’s far left side, an injection well with continuous injection of 80 percent wet steam at 7,000 kPa and production well at 300°C are depicted. On the right is a production well capable of continuous output at 1,000 kPa. The reservoir’s initial temperature is 50°C. Calculated over a period of 1,000 days.

From Figures 1– 5, we can see that the temperature and water saturation leading edges are nearly identical. The temperature leading edge is 70°C at 65 m from the injection well site, and the water saturation is 0.175. The steam leading edge is 35 m from the injection well, at a gas saturation of 0.175, a temperature of 285°C, and a pressure of 6,997 kPa, all of which meet the parameters for steam generation and are consistent with objective facts.

It is feasible to divide the reservoir into three regions: steam flooding, hot water flooding, and cold water flooding, which corresponds to prior knowledge of steam flooding calculation and suggests that the method described in this study is reliable and effective.