Langmuir’s isothermal adsorption equation is (Langmuir, 1917; Chen et al., 2016; Zhao et al., 2016; Xu et al., 2020). V = V L p p L + p (1) Where: V _L is Langmuir volume, sm³/m³; p _L is Langmuir pressure, Pa; p is Reservoirs pressure, Pa; V is Shale gas adsorption volume, sm³/m³.

In order to find the linear governing equation, pseudo-pressure is defined as: ψ = 2 ∫ p 0 p p μ Z d p (2) Where: μ is CBM viscosity, cp.; Z is gas deviation factor, dimensionless; ψ is pseudo-pressure, Pa/s.

For the convenience of mathematical model solving, the following dimensionless variable definition is given in Table 1.

Where: r is the radial distance in x-y plane, m; h is reservoirs thickness, m; T _sc is temperature under standard conditions, K; T is reservoirs temperature, K; p _sc is pressure under standard conditions, Pa; ϕ _m is matrix porosity, decimal; ϕ _f is natural fracture porosity, decimal; C _mt is matrix system compressibility coefficient, Pa⁻¹; C _ft is natural fracture system compressibility coefficient, Pa⁻¹; C _FD is Dimensionless fracture conductivity, dimensionless; μ is the CBM viscosity, cp; q ˜ ¯ is the continuous unite length source strength, m³/s; k _f is the permeability of natural fracture system, m²; k _m is the permeability of matrix system, m²; k _F is the permeability of hydraulic fracture, m²; θ is the angle from linear source to center of reservoirs, degree; L _ref is reference length, m; L _F is fracture length, m; α _m is matrix shape factor, dimensionless; R _e is closed boundary radius, m; r _o is off-center distance, m; x _asmy is distance from the center of hydraulic fracture to wellbore, m;ω is the storativity-ratio, dimensionless; λ is interporosity coefficient, dimensionless; subscript D is the dimensionless.

In cylindrical coordinate system r-θ-z, the location of the line source is assumed to be (r’,θ′,z’). According to the work of Guo et al. (2016) and Xu et al. (2020), he two-dimension governing equation in the radial cylindrical system (Ozkan, 1994; Deng et al., 2017) is 1 r D ∂ ∂ r D ( r D ∂ Δ ψ ¯ f ∂ r D ) + 1 r D 2 ∂ ∂ θ ( ∂ Δ ψ ¯ f ∂ θ ) + 1 r D q ˜ ¯ L ref π k f h T p sc T sc δ ( θ − θ ′ ) δ ( r D − r ′ D ) = f ( s ) Δ ψ ¯ f (3) Where: f ( s ) = { ω s + s β m λ ( 1 + σ ) ( 1 − ω ) β m λ + s ( 1 + σ ) ( 1 − ω ) ,   T r a n s i e n t   s t a t e   ω s + β m λ 5 [ 15 s ( 1 + σ ) ( 1 − ω ) β m λ coth ( 15 s ( 1 + σ ) ( 1 − ω ) β m λ ) − 1 ] , P s e u d o - s t e a d y   s t a t e

Where: r’ is the distance from linear source to the center of reservoir, m; θ’ is the angle from linear source to the center of reservoir, degree; β _m is matrix apparent permeability coefficient, dimensionless; σ is adsorption gas desorption coefficient, dimensionless.

The general solution of Eq. 3 is (Ozkan, 1994). Δ ψ ¯ f = P + E (4) Where: P is a center linear source solution of coal gas reservoirs, this solution can be obtained easily. E is chosen such that P + E satisfies the boundary condition at r = R _e and the contribution of E to the flux vanishes as θ → θ ′ and r → r ′ .

The position relationship between the calculation point and the vertical line source is shown in Figure 3.

According to the previous research work, the center linear source solution P is (Ozkan and Raghavan, 1991) P = μ q ˜ ¯ 2 π k h D L r e f [ K 0 ( R p D f ( s ) ) + D 0 I 0 ( R p D f ( s ) ) ] (5) Where: R p D = r D 2 + r ′ D 2 − 2 r D r ′ D ⁡ cos ( θ − θ ′ )  

According to the addition theorem of the Bessel function (Carslaw and Jaeger, 1959) I k ( a R p D ) = ∑ k = − ∞ + ∞ ( − 1 ) k I k ( a r D ) I k ( a r ′ D ) cos ⁡ k ( θ − θ ′ )     r D < R e D (6) K 0 ( a R p D ) = { ∑ k = − ∞ + ∞ I k ( a r D ) K k ( a r ′ D ) cos ⁡ k ( θ − θ ′ )        r D < r ′ D ∑ k = − ∞ + ∞ I k ( a r ′ D ) K k ( a r D ) cos ⁡ k ( θ − θ ′ )        r D ≥ r ′ D (7) Where: I _k (x) is the k order first type Bessel function; K _k (x) is the k order second type Bessel function.

Therefore, E can be written as E = μ q ˜ ¯ 2 π k f h D L r e f 2 ∑ k = 1 ∞ D k , 0 I k ( a r D ) I k ( a r ′ D ) cos ⁡ k ( θ − θ ′ ) (8)

Finally, when the vertical linear source is not the center of coal gas reservoirs, the general solution can be written as: Δ ψ ¯ f = μ q ˜ ¯ 2 π k h D L r e f [ K 0 ( R p D f ( s ) ) + D 0 I 0 ( R p D f ( s ) ) + 2 ∑ k = 1 ∞ D k I k ( a r D ) I k ( a r ′ D ) cos ⁡ k ( θ − θ ′ ) ] (9)

For circle closed boundary. ∂ Δ ψ ¯ f ∂ r D | r D = R e D = 0 (10)

For circle constant-pressure boundary. Δ ψ ¯ f | r D = R e D = 0 (11)

Substituting Eq. 9 into Eqs 10, 11, coefficient D _k (k = 0, 1, 2, … ,∞) can be obtained for circle closed and constant-pressure boundary. D k = { − K ′ k ( f ( s ) R e D ) I ′ k ( f ( s ) R e D ) ,   Closed − K k ( f ( s ) R e D ) I k ( f ( s ) R e D ) ,   Constant pressure ( k = 0,1,2 , ... , ∞ ) (12)

Substituting Eq. 12 into Eq. 9 and using with the dimensionless variable definition in Table 1. The linear source is integrated along the fracture direction and a uniform flux surface source solution can be obtained. ψ ¯ f D = 1 2 ∫ L D q ˜ ¯ D [ K 0 ( R p D f ( s ) ) + D 0 I 0 ( R p D f ( s ) ) + 2 ∑ k = 1 ∞ D k I k ( a r D ) I k ( a r p D ) cos ⁡ k ( θ − θ p ) ] d s D (13)

However, integration of Eq. 13 is very difficult directly. Therefore, according to the triangle similarity and Pythagorean law, there is the following mathematical relationship. r ′ D = r o D sin ⁡ θ F sin ( θ F − θ ′ ) ,   ( θ F ≠ 0 ,   π ) (14) Where: θ _F is the angle between hydraulic fracture and horizontal line, degree.

In order to improve the calculation speed and accuracy, the grid needs to be divided into unequal distances (Figure 4). Dividing the hydraulic fracture into 2N segments, and the segment numbers of wellbore left fracture and wellbore right fracture is set as N.

The grid end-point coordinates can be expressed as in local coordinate system (r _F-θ) r F D i = { L F R D ( N − i + 1 ) N , ( 1 ≤ i ≤ N + 1 ) L F L D ( i − N − 1 ) N , ( N + 2 ≤ i ≤ 2 N + 1 ) (15) Where: L _FRD is the dimensionless right fracture length; L _FLD is the dimensionless left fracture length; N is the fracture segments number. r _FDi is the i-th segment end-point in the local coordinate system.

The grid mid-point coordinates can be expressed as local coordinate system (r _F-θ) r m F D i = r F D i + r F D i + 1 2    ( 1 ≤ i ≤ 2 N ) (16) Where: r _mFDi is the i-th segment mid-point in the local coordinate system.

The grid end-point coordinates can be expressed as in coordinate system (r-θ) { r D j = r o D 2 + r F D i 2 − 2 r o D r F D i ⁡ cos ( π − θ F ) θ j = tan − 1 r F D i ⁡ sin ⁡ θ F r o D + r F D i D ⁡ cos ⁡ θ F (17)

The grid mid-point coordinates can be expressed as in coordinate system (r-θ) { r m D j = r o D 2 + r m F D i 2 − 2 r o D r m F D i ⁡ cos ( π − θ F ) θ m j = tan − 1 r m F D i ⁡ sin ⁡ θ F r o D + r F D i D ⁡ cos ⁡ θ F (18) Where: r _mDi is the i-th segment mid-point in r-θ coordinate system; r _Di is the i-th segment end-point in r-θ coordinate system. θ _mi is the i-th segment mid-point degree in r-θ coordinate system; θ _i is the i-th segment end-point degree in r-θ coordinate system.

However, the pressure drop of the i-th segment is. ψ ¯ f D ( r m D i , θ m i ) = ∑ j = 1 2 N q ˜ ¯ D j ∫ r D j r D j + 1 C θ ( K 0 ( ε 0 [ r m D i ⁡ cos ( θ m i − θ ′ ) − r ′ d ] 2 + r p D 2 ⁡ sin ( θ m i − θ ′ ) 2 ) + D 0 I 0 ( ε 0 r m D i ) I 0 ( ε 0 r ′ D ) + 2 ∑ k = 1 + ∞ D k I k ( ε 0 r m D i ) I k ( ε 0 r ′ D ) cos ⁡ k ( θ m i − θ ′ ) ) d r ′ D (19) Where: C θ = r oD ⁡ sin ⁡ θ F sin ( θ F − θ ′ ) 2 ; θ ′ = θ F − sin − 1 r oD ⁡ sin ⁡ θ F r ′ D

Wellbore pressure of the off-center vertical fractured well with asymmetric infinite conductivity fracture can be written as. ψ ¯ f D = ∑ i = 1 M ψ ¯ f D ( r m D i , θ m i ) (20) Where: M is the fracture number.

The fluid flow of fracture is only considered as linear flow. The two wing lengths of hydraulic fracture are unequal. Fracture tips have no fluid supplement. According to the work of Wang (Wang and Wang, 2014), the solution of hydraulic fracture can be written as: ψ ¯ f D ( r F D , s ) = ψ ¯ f D , a v g + π C F D ∫ − 1 1 G ( r ′ , r F D ) q ˜ ¯ D ( r ′ , s ) d r ′ − 2 π C F D G ( x a s y m D , r F D ) (21) Where: G ( r ′ , r F D ) = { − 1 4 [ ( r ′ + 1 ) 2 + ( r F D − 1 ) 2 − 4 3 ]   - 1 ≤ r ′ ≤ r F D − 1 4 [ ( r ′ − 1 ) 2 + ( r F D + 1 ) 2 − 4 3 ]   - 1 ≤ r ′ ≤ r F D ; ψ ¯ f D , a v g = 1 2 ∫ − 1 1 ψ ¯ f D ( r ′ , s ) d r ′ Where: ψ ¯ fD,avg is the dimensionless average reservoir pressure; C _FD is the dimensionless conductivity. G(x) is the Green function.

The left term of Eq. 21 is the wellbore pressure of infinite conductivity fracture and it can be replaced by Eq. 20. Therefore, the discrete linear equations of Eq. 21 can be written as ∑ j = 1 2 N q ˜ ¯ D j ∫ r D j r D j + 1 C θ ( K 0 ( ε 0 [ r m D i ⁡ cos ( θ m i − θ p ) − r p D ] 2 + r p D 2 ⁡ sin ( θ m i − θ p ) 2 ) + D 0 I 0 ( ε 0 r m D i ) I 0 ( ε 0 r p D ) + 2 ∑ k = 1 + ∞ D k I k ( ε 0 r m D i ) I k ( ε 0 r p D ) cos ⁡ k ( θ m i − θ p ) ) d r p D = ψ ¯ f D , a v g + π C F D ∑ j = 1 2 N ∫ r F D j r F D j + 1 G ( r ′ , r m F D i ) q ˜ ¯ D j ( r ′ , s ) − 2 π C F D G ( x a s y m D , r m F D i ) 1 ≤ i ≤ N (22)

In addition, according to the mass conservation. 1 2 ∑ j = 1 2 N q ˜ ¯ D j = 1 s (23) Where: s is the Laplace variable.

Combining with Eqs 22, 23, the 2N + 1 linear equation can be obtained and solved by the Gauss elimination method. However, Wellbore pressure still cannot be obtained. Taking surface flux and average pressure into Eq. 21 and letting the r _FD = x _asymD, the wellbore pressure of off-center vertical well with asymmetric finite conductivity fracture can be obtained.

Wellbore pressure skin effect and wellbore storage are calculated by Eq. 24 (Van Everdingen and Hurst, 1949). ψ ¯ w D = s ψ ¯ f D + S s + C D s 2 ( s ψ ¯ f D + S ) (24) Where: C _D is dimensionless well storage; S is the skin, ψ ¯ wD is the dimensionless wellbore pressure.

Figure 7 displays the wellbore pressure curve is affected by dimensionless conductivity. The dimensionless conductivity has an obvious influence on wellbore pressure during the linear and bilinear regimes. Large dimensionless conductivity indicates small seepage resistance in hydraulic fracture. Therefore, the larger dimensionless conductivity is the smaller dimensionless wellbore pressure during bilinear and linear flow regime is. Therefore, in the actual fracturing process, it is helpful for improving the production of the CBM wells to increase the fracture conductivity.

Figure 8 displays the rate distribution curve is influenced by dimensionless conductivity. The solid line represents the rate distribution of short fracture and the dots lines represent the rate distribution of the long fracture. The dimensionless conductivity has an obvious influence on rate distribution. Since hydraulic fracture is asymmetric about wellbore, pressure wave around the short fracture wing spreads to reservoir quickly, so rate contribution of short fracture wing is higher than that of length fracture wing during while flow regime. However, as time progresses, since the control radius of short fracture wing is smaller than that of long fracture wing, rate contribution of short fracture wing decreases and rate contribution of long fracture wing increases with time increasing before linear flow regime. Rate contribution of short and long fracture wing keep an approximate constant. Therefore, the large dimensionless conductivity leads to a small rate of short fracture wing and a large of long fracture wing in the whole flow regime.

Figure 9 displays the wellbore pressure curve is affected by the off-center distance. As is shown in Figure 9, the off-center distance has an obvious effect on wellbore pressure during the arc boundary reflection regime. The larger off-center distance indicates that the well is closed to the circle closed boundary. Pressure waves propagate quickly to the boundary and boundary reaction characteristics appear in advance, which leads to the obvious “upwarping” characteristic of the derivative curve. Therefore, the larger the off-center distance is, the earlier the time of the arc boundary reflection regime and derivative curve “upwarping” is.

Figure 10 displays that the wellbore pressure curve is affected by the asymmetry factor. The asymmetry factor has an obvious effect on wellbore pressure during bilinear and linear flow regimes. It is assumed that the fracture length is equal for different asymmetry factors. The larger asymmetry factor indicates that the wellbore is closed to the fracture end-point, which leads that the fluid flows time longer in the long fracture. Therefore, the larger asymmetry factor leads to the larger seepage resistance. That is to say, the larger the asymmetry factor is, the higher of pressure and derivative curve before the radial flow regime is.

Figure 11 displays how the rate distribution curve is affected by the asymmetry factor. The solid line represents the rate distribution of the short fracture and the dots lines represent the rate distribution of the long fracture. The asymmetry factor has an obvious influence on rate distribution. Since hydraulic fracture is asymmetric about wellbore, pressure wave around the short fracture wing spread to reservoir quickly, so rate contribution of short fracture wing is higher than that of length fracture wing during while flow regime. However, as time progresses, since the control radius of the short fracture wing is smaller than that of the long fracture wing, the rate contribution of the short fracture wing decreases and the rate contribution of the long fracture wing increase with time increasing before linear flow regime. For shorter fracture wing, the larger asymmetry factor leads to a larger pressure drop between the wellbore and hydraulic fracture surface. Therefore, the larger asymmetry factor indicates a larger rate contribution of the short wing. On the contrary, the larger asymmetry factor indicates a larger rate contribution of the short wing with consideration of rate conservation.

Figure 12 displays that the wellbore pressure curve is influenced by different crossflow models. As is shown in Figure 12, a different crossflow model has an obvious influence on wellbore pressure during the crossflow regime. Transient state crossflow shows that the CBM of matrix system can crossflow into naturally fracture system quickly when the pressure of the natural fracture system decreases. On the contrary, pseudo-steady state crossflow shows that the CBM of matrix system cannot crossflow into the natural fracture system quickly when the pressure of the natural fracture system decreases. The CBM of the matrix system can flow out the surface of matrix particles by the seepage type. Therefore, the “concavity” of pressure derivative of transient state crossflow is deeper than that of pseudo-steady state crossflow during crossflow regime.

Figure 13 is the rate distribution curve affected by the storativity ratio. The solid line represents the rate distribution of the short fracture and the dotted lines represent the rate distribution of the long fracture. It is assumed that the crossflow model is transient state crossflow. A large storativity ratio indicates larger storativity volume and a more natural network, which indicates more natural fractures can connect with hydraulic fractures. In addition, since the short fracture wing is closer to the wellbore. the rate contribution of short fracture wing is higher than that of long fracture wing. Therefore, the larger storativity ratio leads to the higher rate contribution of short fracture wing and small rate contribution of long fracture wing before crossflow regime.

Figure 14 is the pressure transient response curve affected by adsorption-desorption coefficient. The adsorption-desorption coefficient mainly affects the typical curve shape after the crossflow regime. The smaller the value of adsorption-desorption coefficient is, the narrower and shallower the “concave” on the pressure derivative curve is. On the contrary, the larger the value of the adsorption-desorption coefficient is, the wider and deeper the “concave” on the pressure derivative curve is. The larger the amount of adsorbed gas in the coal matrix is, the more obvious the characteristics of crossflow after desorption are.

Figure 15 is a pressure transient response curve influenced by fracture angle. The influence of fracture angle on the pressure derivative curve is not very obvious. The fracture angle physical model is shown in Figure 16. The blue is the hydraulic and the red circle is the wellbore. Since fluid cannot flow into the hydraulic fracture and the largest rate supplement of the hydraulic fracture comes from fracture surface, the pressure wave can spread to closed boundary quickly with the increasing of the hydraulic fracture angle. Therefore, a larger fracture angle leads to higher pressure derivative curve during the arc boundary reflection regime. However, since the off-center distance is very large, the influence of a larger fracture angle on the pressure derivative curve is not very obvious.

Figure 17 is the rate distribution curve affected by fracture angle. The solid line represents the rate distribution of the short fracture and the doted lines represent the rate distribution of the long fracture. The fracture angle physical model is shown in Figure 16. The fracture angle has no influence on rate distribution curve before the arc boundary reflection regime. Rate supply around the boundary is finite, so short fracture wing rate contribution decrease with the decreasing of the fracture angle during the arc boundary reflection regime.