Rock Mech Rock Eng DOI 10.1007/s00603-017-1202-6

ORIGINAL PAPER

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between Flowing Fluid and Rock Fracture Walls Bing Bai1 • Yuanyuan He1 • Shaobin Hu1 • Xiaochun Li1

Received: 26 August 2016 / Accepted: 8 March 2017  Springer-Verlag Wien 2017

Abstract The convective heat transfer coefficient (HTC) is a useful indicator that characterizes the convective heat transfer properties between flowing fluid and hot dry rock. An analytical method is developed to explore a more realistic formula for the HTC. First, a heat transfer model is described that can be used to determine the general expression of the HTC. As one of the novel elements, the new model can consider an arbitrary function of temperature distribution on the fracture wall along the direction of the rock radius. The resulting Dirichlet problem of the Laplace equation on a semi-disk is successfully solved with the Green’s function method. Four specific formulas for the HTC are derived and compared by assuming the temperature distributions along the radius of the fracture wall to be zeroth-, first-, second-, and third-order polynomials. Comparative verification of the four specific formulas based on the test data shows that the formula A corresponding to the zeroth-order polynomial always predicts stable HTC values. At low flow rates, the four formulas predict similar values of HTC, but at higher flow rates, formulas B and D, respectively, corresponding to the firstand third-order polynomials, predict either too large or too small values of the HTC, while formula C, corresponding to the second-order polynomial, predicts relatively acceptable HTC values. However, we cannot tell which one is the more rational formula between formulas A and C due to the limited information measured. One of the clear advantages of formula C is that it can avoid the drawbacks

& Bing Bai [email protected] 1

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China

of the discontinuity of temperature and the singular integral of HTC at the points (±R, 0). Further experimental work to measure the actual temperature distribution of water in the fracture will be of great value. It is also found that the absorbed heat of the fluid, Q, has a significant impact on the prediction results of the HTC. The temperatures at the inlet and the outlet used for Q should be consistent with the assumptions adopted in the derivation of its corresponding HTC formula. A mismatched value of Q might be the reason that some existing HTC formulas predict negative or extremely large HTCs at high flow rates. Keywords Convection heat transfer coefficient  Rock fracture  Granite  Hot dry rock (HDR)  Water flow List of symbols Ti (x) Inner surface temperature distribution of the rock sample (K) Tw (x) Temperature distribution of the fracture water (K) Tw1 (x) Temperature distribution at the fracture inlet Tw2 (x) Temperature distribution at the fracture outlet Ti0 Inner surface temperature at the center of the rock sample (K) T1 The temperature of water at the fracture inlet (K) T2 The temperature of water at the fracture outlet (K) L Length of the rock sample (m) k Thermal conductivity of rock (W/(m K)) d Width of aperture (m) HTC Heat transfer coefficient n The order of polynomial

123

B. Bai et al.

f1 (x) f2 (x) h R Q Tc Cp q u G x y u (x, y) b, bw1, bw2 Z?

The boundary condition at L2[-R, 0], in ‘‘Appendix’’ The boundary condition at L3[0, R], in ‘‘Appendix’’ Heat transfer coefficient (W/(m2 K)) Radius of the rock specimen (m) Heat quantity transferred to the water (J) The temperature at the outer wall surface of specimen, i.e., the oil temperature (K) Specific heat capacity at constant pressure (J/(kg K)) Density of water (kg/m3) Flow velocity (m/s) Green’s function Horizontal coordinate Vertical coordinate Solution to Dirichlet problem of Laplace equation on a semi-disk, in ‘‘Appendix’’ Coefficients in polynomial The set of positive integers

1 Introduction The convection heat transfer coefficient (HTC) is a key indicator that characterizes the convective heat transfer properties between flowing fluid and hot dry rock (HDR) (Mohais et al. 2011; Zhao 2014; Xu et al. 2015). It is required in both numerical simulations of enhanced geothermal system (EGS) recovery and the design of geothermal power plants (Ogino et al. 1999). Flow and heat transfer tests in rock fractures are crucial to understand the physical process and to verify the prediction ability of the heat transfer coefficient (Xu et al. 2015; Lu and Xiang 2012). As a typical implementation, conducting the tests based on the conventional triaxial test platform with the capabilities of high temperature and fluid flow components is considered a convenient and effective approach (Zhao 1987, 1992; Zhao and Brown 1992; Zhao and Tso 1993). In the tests, fluid flows through and exchanges heat with the fracture wall of a cylindrical rock specimen, and the fracture aperture can be controlled by varying the confining pressure. The convective heat transfer coefficient can be obtained with the measured parameters (e.g., temperature and geometry) and the prediction models. To derive the formula of the convective heat transfer coefficient for such tests, one needs to first obtain the temperature distribution on the fracture wall. However, so far, existing temperature solutions have been obtained under highly simplified conditions. In an early formula for the HTC by Zhao (1992), the semi-disk was simplified into an equivalent rectangle, and the temperature along the

123

rectangular wall was assumed to be constant. A recent comparison of the existing formulas shows that the resulting prediction formula can produce abnormal results at higher flow rates (Bai et al. 2017). Subsequently, a new formula was developed for which the shape approximation of the rock was eliminated, but the assumption of constant temperature distribution on the fracture wall was maintained (Zhao 1999). The same approach was recently investigated by Zhang et al. (2015) and Huang et al. (2016). Comparisons of the existing formulas of the heat transfer coefficient made by Bai et al. (2017) demonstrate the stability of this formula. However, the reason behind this stability is still not well understood, and the assumption of constant temperature distribution leads to temperature discontinuity at the points (-R, 0) and (R, 0), which further produces singular integrals at the points. Recently, Zhao (2014) presented a new prediction formula based on a two-dimensional analytical solution that they developed. However, the formula predicted some negative values that cannot be explained. The actual temperature distribution on the fracture wall along the direction of radius is not constant. When the diameter of the specimen increases, the assumption of a constant temperature distribution can increase the deviation of HTC. However, little is known about the realistic distribution of temperature along the fracture wall, as no measured data are currently available, and therefore the effect of a non-constant distribution of temperature on the HTC is not understood. An analytical solution that can accommodate an arbitrary temperature distribution of temperature on the boundary will not only advance the understanding of the temperature distribution and its effect on the HTC but also help to develop a more reasonable formula for the HTC. The aim of this paper is to develop an analytical method for the HTC that can consider the arbitrary distribution of temperature along the fracture wall. For this purpose, in Sect. 2, a heat transfer model is first described, and a general formula for the convection HTC is presented. Then, in Sect. 3, some specific analytical solutions for the heat transfer coefficient are developed corresponding to different distribution functions of temperature. To derive these formulas, accurate solutions of temperature of the heat conduction equation in the semi-disk region are required. The heat conduction problem is characterized as a Dirichlet problem of the Laplace equation on the semidisk, which is solved with the method of Green’s functions. The detailed derivation process is shown in ‘‘Appendix.’’ In Sect. 4, we present the analysis and discussion conducted on the new solution by comparison and verification. Finally, in Sect. 5, we summarize our primary findings and present the study’s conclusions.

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between…

2 Description of the Heat Transfer Model 2.1 Experimental Principle The heat transfer coefficient model is based on a flow and heat transfer experiment conducted on a cylindrical specimen with a single fracture. Figure 1 shows the experimental principle and how the sensors are set up. This type of test can be implemented on a function-enhanced triaxial test system. The fractured specimen is placed inside a pressurization chamber with heating capability. The specimen is encapsulated in a heat-shrinkable sleeve, and hydraulic oil is heated to a preset temperature. During the process of fluid flowing in from one side and out the other side, the fluid exchanges are heated with the fracture surface. The heat transfer coefficient can be calculated with a model and the parameters measured in the experiment. For this study, a series of water flow and heat transfer tests were conducted at three confining pressures (0, 3, and 6 MPa) at the confining temperature of 90 C. The fractured specimen was obtained by splitting off a U50 mm 9 100 mm cylindrical granite rock using the Brazil splitting test method (ISRM 1978). More technical details and experimental results can be found in recent work by Bai et al. (2017). The width of the fracture aperture under the confining pressure can be obtained by subtracting the lateral deformation of the specimen from the initial width (Bai et al. 2017). The image measurement approach was used to obtain the value of the initial width of the fracture aperture. The lateral deformation recorded by the lateral strain gauge includes the deformations of the fracture and the two halves of the granite specimen. However, the deformation of the granite specimen can be ignored because the

confining pressures are not high and the stiffness of the granite specimen is considerably larger than that of the fracture. The uncertainty of the size measurement is ±0.01 mm. 2.2 Description of the Heat Transfer Model The actual problem is, of course, three-dimensional; in this study, however, we are investigating the average heat transfer coefficient of the whole fracture because it is much more difficult to obtain the analytical solution of the threedimensional model. Therefore, the three-dimensional specimen is simplified into an equivalent two-dimensional model, as shown in Fig. 2. Curvature of the fracture wall will be ignored to make the heat conduction equation solvable analytically. The seepage of fluid into the rock is neglected because the permeability and the porosity of the granite specimen are extremely low. In addition, no radiation heat transfer or no heat loss will be considered as heatinsulating jacket is wrapped around the test cell (Bai et al. 2017). Tc is the temperature at the outer wall surface of the specimen, i.e., the oil temperature (K). Ti (x) is the temperature distribution of the inner surface of the rock sample (K). Tw (x) is the temperature distribution of the fracture water. Ti0 is the inner surface temperature at the center of the rock sample (K). R is the radius of the rock specimen (m). For a differential segment dx, the heat transferred from the fracture wall to the water is dQ ¼ hðTi ðxÞ  Tw ðxÞÞdA ¼ hðTi ðxÞ  Tw ðxÞÞLdx

ð1Þ

where h is the heat transfer coefficient (W/(m2 K)). Therefore, h¼

4L

RR 0

Q

ð2Þ

ðTi ðxÞ  Tw ðxÞÞdx

It can be observed from this expression that the temperature functions Ti (x) and Tw (x) and the heat quantity Q are y

Tc

Confining pressure Semi-disk rock specimen

Fracture wall Fracture -R Tw= Tw(x)

Fig. 1 Experimental schematic diagram (Bai et al. 2016, 2017)

Ti0= Ti(x=0)

R

x

Ti= Ti(x)

Fig. 2 Conceptual model of the flow and heat transfer in a half-disk rock specimen (Bai et al. 2017)

123

B. Bai et al.

required to calculate the heat transfer coefficient h. Firstly, let us discuss how to determine Ti (x) and Tw (x). The specific forms of these two functions have to be assumed because we have no knowledge about them. Although theoretically any functional form can be used in our model, here an nth-order polynomial will be employed for both Ti (x) and Tw (x) to derive the model without loss of generality. In addition, the candidate functions should be even functions because of the symmetry with respect to the y-axis. The temperature distribution function on the fracture wall of the half-disk is assumed to be Ti ðxÞ ¼ Ti0 þ bj xjn

Tc  Ti0 n j xj Rn

ð4Þ

In Fig. 2, Tw (x) is the temperature distribution on the middle line of the 2-D fracture, which is a compact equivalent of the actual 3-D fracture. The middle line is actually a compact equivalent of the actual middle section of the fracture. Thus, Tw (x) should be an equivalent approximation of the temperature on this middle section. As we do not know the actual temperature distribution on the entire specimen section, Tw (x) can be approximated as the average of the temperature distributions Tw1 (x) and Tw2 (x) at the fracture inlet and the fracture outlet, respectively, i.e., Tw ð xÞ ¼ ðTw1 ðxÞ þ Tw2 ðxÞÞ=2

ð5Þ

Similar to the expression of the fracture temperature in (3), Tw1 (x) and Tw2 (x) are, respectively, Tw1 ðxÞ ¼ aw1 þ bw1 j xj

n

ð6Þ

Tw2 ðxÞ ¼ aw2 þ bw2 j xjn

ð7Þ

According to the experimental principle, we can express the boundary condition at the fracture inlet as  Tw1 ðx ¼ 0Þ ¼ T1 ¼ aw1 ð8Þ Tw1 ðx ¼ RÞ ¼ Tc ¼ Tw1 ðxÞ ¼ aw1 þ bw1 Rn where T1 is the temperature of water at the fracture inlet (from the ‘‘2-inlet thermometer’’ in Fig. 1). Substituting the expressions of aw1 and bw1 solved from (8) to (6) yields Tc  T1 n ð9Þ j xj : Rn At the fracture outlet, because the outlet thermometer is placed in the hole of the upper pad downstream of the fracture outlet, the temperature from this thermometer T2 is treated as the temperature of water at (R/2, 0). Moreover,

123

Tw2 ðxÞ ¼ Tc 

Tc  T2 Tc  T2 þ n j xjn n 12 R ð1  2n Þ

ð11Þ

Equation (5) can be updated by substituting (9) and (11). Substituting the updated Eq. (5) together with Eq. (4) into Eq. (2), we have ð2n  1Þðn þ 1ÞQ ; 2nLRðT1 þ Tc  2n ðT1 þ T2  2Ti0 Þ  2Ti0 Þ n 2 Zþ:

h¼

ð12Þ

We now discuss the determination approach for Ti0. As far as the distribution function of the temperature of the fracture wall is assumed [i.e., Eq. (4)], we can obtain the analytic solution of the heat conduction equation on the semi-disk region (to be presented in Sect. 3 and ‘‘Appendix’’). We denote the solution as T(x, y; Ti0). According to the principle of conservation of energy, the heat from the fracture wall should be equal to that absorbed by the water, i.e.,    Z R   oT ðx; y; Ti0 Þ   dx ¼ Q 4kL  oy 0 y¼0 

ð13Þ

This equation can be used to determine the parameter Ti0, which should be a function of the heat quantity Q: Ti0 ¼ Ti0 ðQÞ

and

Tw1 ðxÞ ¼ T1 þ

Substituting the expressions of aw2 and bw2 solved from (10) to (7) yields

ð3Þ

The continuity of temperature (R, 0) requires Ti (x = R) = Tc, which leads to b = (Tc - Ti0)/Rn. Therefore, Eq. (3) can be restated as Ti ðxÞ ¼ Ti0 þ

the temperature of water is considered to be continuous at (R, 0) with the hydraulic oil. Thus, we have the following boundary conditions of water temperature at the outlet: 8    n R R < ¼ T2 ¼ aw2 þ bw2 Tw2 x ¼ ð10Þ 2 2 : Tw2 ðx ¼ RÞ ¼ Tc ¼ aw2 þ bw2 Rn

ð14Þ

Here we only show the approach of how to obtain Ti0; its explicit expression will be derived later in Sect. 3 when the order of the polynomial, n, that the explicit solution T(x, y; Ti0) depends on is assigned a particular value. We now address the calculation of transferred heat quantity Q. In previous studies (Zhao 1999; Zhang et al. 2015), Q was obtained by Q ¼ cp qm DTw ¼ cp qm ðT2  T1 Þ;

ð15Þ

where qm ¼ 2qudR: In Eq. (15), constant temperature distributions, T1 and T2, are assumed at the inlet and the outlet of the fracture, respectively. However, in our model, the temperature distributions are no longer constant. Therefore, Q can be calculated with the following formula instead:

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between…

Q¼

Z

An equation on Ti0 can be obtained by substituting Eq. (20) into the left-hand term of Eq. (13) as follows:

R

cp ðTw2 ðxÞ  Tw1 ðxÞÞ2quddx

0

¼ 2cp qudR

nðð2n  1ÞT1  2n T2 þ Tc Þ ð 2n  1Þ ð n þ 1Þ

ð16Þ

where the right-hand term is obtained by substituting Eqs. (9) and (11) into the integral of the left-hand term.

16kLðTc  Ti0 Þ ¼Q p

ð21Þ

We can solve for Ti0 by solving Eq. (21): Ti0 ¼ Tc 

pQ 16kL

ð22Þ

3 Specific Analytical Solution for Heat Transfer Coefficient

Finally, the solution of h can be obtained by substituting Q and Ti0 into Eq. (12):

The target of this section is to derive the specific analytical solution for the heat transfer coefficient. As noted previously, the analytic solution of the heat conduction equation on the semi-disk region T(x, y; Ti0) should be obtained first. The heat conduction equation is a special case of the Dirichlet problem of the Laplace equation on a semi-disk, which is solved in ‘‘Appendix.’’ Here, we present derivation details of the HTC by using n = 2 as an example. When n = 2, a quadratic polynomial is used as the temperature distribution function of the fracture wall and the water. The heat conduction equation and its boundary conditions are

h¼

DT ¼ 0; in X; oX ¼ L1 þ L2 þ L3 8 on L1 < Tc ; T¼ Tc  Ti0 2 : Ti0 þ x ; on L2 þ L3 R2

where the absorbed heat Q can be obtained by letting n = 2 in Eq. (16): 4 Q ¼ cp qudR½4T2  3T1  Tc : 9

ð24Þ

Using the same process, the solutions corresponding to n = 1 or n = 2 are also obtained. All the solutions are listed in Table 1. It should be noted that the solution corresponding to n = 0 is also listed for later comparison although it is not covered in the general formula (12).

4 Analysis and Discussion 4.1 Verification of the Solution of the Heat Conduction Equation

ð18Þ

Equation (17) can be transformed to the following by substituting (18), 8 > < Du ¼(0; in X; oX ¼ L1 þ L2 þ L3 0; on L1 Tc  T0i 2 u ¼ > : Ti0 þ x  Tc ; on L2 þ L3 R2

ð23Þ

ð17Þ

where Ti0 is a parameter to be determined representing the temperature at the center of the fracture wall. Let u ¼ T  Tc

9Q 4LRð6Ti0 þ Tc  3T1  4T2 Þ 18kQ ¼ R½8kLð7Tc  3T1  3T2 Þ  3pQ

ð19Þ

This equation is the Dirichlet problem of Laplace equation on a semi-disk in ‘‘Appendix’’ whose solution is (40). Therefore, the solution of Eq. (19) can be obtained by using (40) as, 1 T ¼ Tc þ yðTc aÞ p      2 3 2ðR2 x2 y2 Þ R2 x2 þy2 Rx Rþx þArctan Arctan 6 Rðx2 þy2 Þ  7 2 y y R y 6 7 6     7 2 2 2 2 2 2 2 2 2 2 6 R ðx y Þ ðx þy Þ 7 Rxþx þy Rxþx þy 6 7 Arctan þArctan 6 7 Ry Ry yðx2 þy2 Þ2 6 7 6 x  7    6 7 2 2 2 2 6  2 ln ðRþxÞ þy ln ðRxÞ þy 7 6 R 7 6 7     4 5 R2 x   2  2 2 2 2 2 þ ln R ðRþxÞ þy ln R ðRxÞ þy 2 ðx2 þy2 Þ

ð20Þ

As emphasized in previous sections, the correctness of the analytical solution of the heat conduction Eq. (17) is a prerequisite to ensure a rational formula for the heat transfer coefficient. In this section, we verify the solution by comparing it with a numerical solution based on the finite element method. As we cannot verify all the solutions from the general solution suitable for arbitrary boundary functions in ‘‘Appendix,’’ we only verify the solution (20) as an example, where the boundary function on the diameter is in the form of a ? bx2. In the verification example, Ti0 = 30 C, Tc = 80 C, and R = 0.025 m. Figure 3 shows the contour lines of the analytical solution of the example. Figure 4 shows a comparison of the analytical and the numerical solutions at selected points. The numerical solutions were obtained with the finite element method on the COMSOL Multiphysics Platform by directly inputting the Dirichlet problem of the Laplace Eq. (17) with the same parameters. In the numerical simulation, the half-disk region was discretized with 274 domain elements of triangle and 40

123

123

Order of polynomial

0

1

2

3

No.

A

B

C

D

8T2 Tc 7

2 þ 87 TcRT j xj 3 2

1 T1 þ TcRT j xj3 3

R

2 2 þ 43 TcRT x 2   Tc Ti0  3  Ti0 þ 3 x

4T2 Tc 3

1 2 T1 þ TcRT x 2

i0 2 x Ti0 þ Tc RT 2

2T2  Tc þ 2ðTcRT2 Þ j xj

1 T1 þ Tc T R j xj

i0 Ti0 þ Tc T R j xj

T2

T1

Ti0

Temp. distribution function on y = 0 Ti ¼ Tw1 ¼ Tw2 ¼

Tc  8kLðpQ 1þln 4Þ

pQ Tc  16kL

 ð8T2  7T1  Tc Þ

 ð4T2  3T1  Tc Þ

3 14 Cp dRqu

4 9 Cp dRqu

Cp dRqu  ð2T2  T1  Tc Þ

2Cp qudRðT2  T1 Þ

Tc  0:074Q kL

Tc  4kLpQ ln 16

Quantity of heat absorbed by water Q=

Temperature at the center of the semi-disk specimen Ti0 ¼

Table 1 HTC solutions corresponding to different temperature distribution functions on the fracture wall

14Q 3LR½Tc þ14Ti0 7T1 8T2 

9Q 4LRð6Ti0 þTc 3T1 4T2 Þ

ð1þln 4ÞkQ ¼ 3R½4kLð1þln 56 4Þð15Tc 7T1 8T2 Þ7pQ

18kQ ¼ R½8kLð7Tc 3T 1 4T2 Þ3pQ

ln 16 ¼ R½2 ln 16kLð2kQ 3Tc T1 2T2 ÞpQ

¼ R½0:5kLð2Tc 0:25kQ T1 T2 Þ0:074Q

Q LRðTc þ2Ti0 T1 2T2 Þ

Q 2LR½2Ti0 T1 T2 

Convection heat transfer coefficient h¼

B. Bai et al.

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between…

Vertical position (m)

80 °C 74.4 °C 68.2 °C 62.0 °C 55.8 °C 49.6 °C 43.4 °C 37.2 °C

Horizontal position (m)

Fig. 3 Contour maps of the analytical solution of the example

Fig. 5 The fractured specimen used in the experiment

Temperature (°C)

corresponding to the confining temperature of 90 C will be used as working data for the analytical formulas. These test results are listed in Table 2. 4.2.2 Comparisons of the Different HTC Solutions

Position to the center of the specimen (m)

Fig. 4 Comparison of the analytical and the numerical solutions at selected points

boundary elements. Additionally, various meshes from the extremely coarse element size (16 domain elements and 10 boundary elements) to extremely fine element size (11,218 domain elements and 258 boundary elements) were simulated. The results show that with the increase of the number of the elements, the numerical results rapidly converge to the corresponding analytical solutions. Overall, the numerical and the analytical solutions are highly consistent. Therefore, we conclude the correctness of our analytical solution has been thoroughly confirmed. 4.2 Comparative Verification of the HTC Formulas 4.2.1 Experimental Results Following the experimental principle presented previously in Sect. 2.1, a series of water flow and heat transfer tests were conducted at various confining pressures (0, 3, and 6 MPa) and various confining temperatures in the very recent work by Bai et al. (2017). The fractured specimen was obtained by splitting off a U50 mm 9 100 mm cylindrical granite rock using the Brazil splitting test method. Figure 5 shows the fractured specimen used in the experiment. In this paper, only the test results

Figure 6 shows the comparison of the predicted HTC solutions using test data under the confining pressures of 0, 3, and 6 MPa. During the calculations of the HTCs, the specific heat capacity at constant pressure of water 4200 (J/ (kg K)) is used to calculated the heat quantity Q transferred to the water. Other parameters involved in the formulas of h in Table 1 are shown in Table 2. Overall, the plots of the HTCs based on the three sets of test data under the confining pressures show similar characteristics. For each confining pressure, some solutions predict abnormal (i.e., negative or extremely large) HTC values. This finding suggests that these solutions are not suitable for this test. Interestingly, of these solutions, this solution and the one by Zhao (1992) predict similar trends with increasing flow rates, although the HTC values are different. The other solutions show a roughly positive correlation between the HTC and the flow rates. They predict close HTC values at low flow rates; at higher flow rates, different solutions predict significantly different results. The lower the order of the polynomials, the smaller the HTC they will predict at higher flow rates. One exception is the solution A (in Table 1), which is a zeroth-order polynomial, that predicts larger values of HTC than that predicted by solution C based on second-order polynomial. Actually, solution A predicts the most stable HTC for all the test cases. To explain this phenomenon, we investigate the temperature distribution on the straight edge when different polynomials are assumed. Without loss of generality, we choose only one test case at the confining pressure 3 MPa and injection flow rate = 15 mL/min as an illustration example. Other parameters used include Tc = 90 C, L = 0.1 m, R = 0.025 m, and k = 2.78 W/(m K).

123

B. Bai et al. Table 2 Temperature outside the specimen: 90 C (Bai et al. 2017) Flow rate (mL/ min)

Confining pressure = 0 MPa

Confining pressure = 3 MPa Pressure drop (kPa)

Average aperture (lm)

Average aperture (lm)

Inlet temp (C)

Outlet temp (C)

5

125.5

78.8

89.2

71.0

95.0

83.2

89.4

463.3

80.5

83.1

89.6

749.5

8

125.5

76.4

88.3

102.0

95.0

78.8

89.2

658.7

80.5

80.2

89.4

1304.5

10

125.5

74.1

87.7

121.0

95.0

75.1

88.9

818.7

80.5

76.6

88.9

2000.5

13

125.5

70.4

86.6

150.0

95.0

71.0

87.9

1053.7

80.5

73.1

87.9

2834.5

15

125.5

67.5

85.1

169.0

95.0

67.8

86.4

1196.5

80.5

69.2

87.6

3439.5

18

125.5

64.3

83.2

200.0

95.0

63.3

84.9

1413.3

80.5

65.1

86.3

4329.0

20

125.5

62.4

81.4

223.0

95.0

60.9

83.3

1548.3

80.5

62.8

83.5

4852.0

25

125.5

58.6

79.5

280.0

95.0

56.7

80.8

1918.7

80.5

59

81.8

6070.0

30

125.5

55.5

77.1

340.0

95.0

52.8

77.6

2276.9

80.5

55.2

79.1

6316.5

35

125.5

53.4

75.0

402.0

95.0

50.6

75.3

2589.6

80.5

52.1

76.4

6385.0

Fig. 6 Comparisons of the HTC solutions with the test data

Inlet temp (C)

Confining pressure = 6 MPa

Outlet temp (C)

Pressure drop (kPa)

h

heat transfer coefficient (W / (m2. K))

D:x3

C:x2

Confining pressure = 0 MPa

5

Average aperture (lm)

Inlet temp (C)

5

Pressure drop (kPa)

A:x0

B:x

Confining pressure = 6 MPa

Confining pressure = 3 MPa

35

Outlet temp (C)

35

5

35

Flow rate (mL/min)

Figure 7 shows the temperature distributions along the radius [0, R] of all the functions we investigated previously. The intercept on the vertical axis is Ti0. The constant function has the highest value of Ti0. With the increase of the order of the polynomials, the curves of the functions gradually approach the constant function and the predicted value of Ti0 also gradually increases. These curves can be observed as the predicted temperature distributions on the rock walls. Therefore, different orders of polynomials will induce significant differences in the temperature distributions. When the radius of the specimen increases, the difference will also increase. 4.3 The Effect of Absorbed Heat Quantity Q In the formulas of h listed in Table 1, comparing to the L, R, and Tc, which are known constants, Q relies on additional formulas that are not unique. For example, Eqs. (15) and (16) are two options, and more options are possible as

123

Fig. 7 Temperature distributions along the radius of the rock wall corresponding to different polynomial assumptions

far as they are rational. In this section, we investigate the effect of the absorbed heat quantity Q on heat transfer coefficient h.

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between…

Q Heat quantity transferred to water

Figure 8 shows the four groups of heat quantity Q corresponding to the four formulas in Table 1 using the test results under a confining pressure of 0 MPa in Table 2. The four formulas clearly predict remarkably different values of Q. Formula A, corresponding to the constant temperature on the fracture wall, predicts the highest values. In the left three formulas, the formulas with higher-order polynomials predict larger values of Q. We note that heat quantity indeed has crucial effect on h. During the initial development process of our model, the Q in formula A was used in all four formulas of h. As a result, however, all three formulas produced extremely abnormal and unphysical (i.e., extremely large or negative) values of h, which can be verified from Fig. 9. When their own formulas of Q were used, the three formulas of B–D produced much more reasonable values of h, which is shown in Fig. 6. Therefore, the formula for h and the formula for Q should be matched. The reason that Q has a significant effect is that it exists in the denominator of the formula of h. From the last

column of Table 1, it can be observed that the denominators of all the formulas of h include a subtraction operation with Q as the subtrahend. Simultaneously, for a given test record, the minuend is constant. Therefore, the variation of Q can induce qualitative changes in the relationship between subtrahend and minuend. The first possibility is that when Q grows, the subtrahend will be close to the minuend. This will produce a very small denominator and thus an extremely large value of h, but it possibly includes numerical errors and is not believable. A second possibility is that when Q made the subtrahend greater than the minuend, a negative h is produced, which is unacceptable. According to Fig. 8, formula A, corresponding to a constant temperature on the fracture wall, predicts much higher values of Q than other formulas; therefore, if this formula of Q is used in formulas B–D of h, abnormal values of h might be produced, in agreement with the results shown in Fig. 9. In addition, comparing to the Q that appears in the denominator, the Q that appears in the numerator of the formulas of h has no contribution to producing abnormal values of h. It has been found that some existing formulas of h also predict extremely large or negative values of h (Bai et al. 2017), especially at high flow rates that imply larger values of Q. Therefore, this result indicates that the calculation methods used there might not match with the formulas of h. Direct use of Eq. (15) in various formulas of h might overestimate the values of Q corresponding to their model of h. 4.4 On the Assumptions of the Model

Flow rate (mL/min)

Fig. 8 Q values of different solutions

Assumptions, which cannot be avoided in any mathematical model, represent in a certain sense the spirit and the limitation of a model. Here we emphasize the main assumptions of the

h

heat transfer coefficient (W / (m2. K))

Fig. 9 HTC solutions with fixed heat flow

Flow rate (mL/min)

123

B. Bai et al.

heat transfer model, which are different from that in the previous studies. As one of the novel elements of our model, variable temperature distributions of fracture wall and water along the direction of radius are used to update the assumption of constant temperatures in existing studies. Although we consider the polynomial function to be simple and sufficient to describe the practical temperature distributions during the derivation of the model in Sect. 2, theoretically the methodology described there is suitable to any analytical function that is needed. Other notable assumptions are the boundary conditions of the water temperatures at the inlet and the outlets, i.e., Equations (8) and (10). These assumptions are required because temperatures T1 and T2 from the probes are actually the temperatures of the water in the holes of the down pad and the up pad, respectively, according to the experimental principle in Fig. 1. Clearly, T1 is not the temperature of the water in the inlet fracture, and T2 is not the temperature of the water in the outlet fracture. As the position of the inlet probe of the water is close to the center of the fracture, T1 is approximated to be the temperature of the water at the center of the fracture, and therefore we have the first boundary condition Tw1 (x = 0) = T1 = aw1 in Eq. (8). The second boundary in Eq. (8) involves the temperature of the water at (R, 0). As shown in Fig. 1, there is a thin heat-shrinkable film between the water and the outside hydraulic oil at (R, 0), and there might be a temperature drop across the film. However, this temperature difference was ignored in the second boundary Tw1 (x = R) = Tc = aw1 in Eq. (8) because no measured water temperature is available there. The same approximation was used in the second boundary condition at the outlet in Eq. (10). However, a different process was used in the first boundary condition in Eq. (10). Different from T1, T2 is the temperature of the water downstream of the fracture outlet. Physically, we considered T2 to be the temperature of the mixed water flowing out of the fracture, and T2 was approximated as the temperature at x = R/2 rather than x = 0, i.e., Tw2 (x = R/2) = T2 in the first boundary condition in Eq. (10). We note that different definitions of T2 will lead to entirely different expressions of heat quantity Q and thereby have a significant effect on the results of the heat transfer coefficient. When T2 was treated as Tw2 (x = R/2) = T2, overestimated heat produced abnormal (extremely large or negative) values of heat transfer coefficient. In any case, a more reasonable determination of the boundary conditions of the water relies on additional measurements of the water temperature at the inlet and the outlet fractures.

5 Conclusions This study presents an analytical method developed for evaluating the convection heat transfer coefficient between flowing fluid and rock fracture walls. Based on this

123

method, derived formulas of the heat transfer coefficient were compared, and further discussion was presented on the reason why certain formulas produced abnormal values for the heat transfer coefficient. Some specific valuable conclusions are drawn as follows: 1.

2.

3.

The heat transfer model describing how to determine the HTC was detailed in this paper. Although the methodology framework used here in our model is the same as that in the existing studies, the constant distribution of temperature on the fracture wall was replaced with an arbitrary function along the direction of the rock radius. Thus, the new solution, on the one hand, contains the existing formula; on the other hand, it also accommodates a variety of new solutions. If a continuous function of temperature is used, the drawbacks of the discontinuity of temperature and the singular integral at the points (±R, 0) of HTC can be overcome, which will also provide the possibility of presenting a realistic distribution of temperature on the rock fractures. Generalized boundary conditions of temperature on the fracture wall require the analytical solution of the steady heat conduction equation in a semi-disk region, which is actually a special case of the Dirichlet problem of the Laplace equation on a semi-disk. The analytical solution of this problem was developed with the Green’s function method, which can consider a general Dirichlet boundary condition of arbitrary functions. The accuracy of the analytical solution of temperature has been verified by the almost identical numerical results based on the finite element method. The temperature solution was then used as the core to derive the specific formula of the heat transfer coefficient. We derived four specific formulas of HTC by assuming the temperature distributions along the radius of fracture wall to be zeroth-, first-, second-, and thirdorder polynomials. Without loss of generality, only the derivation process of the formulas corresponding to the second-order polynomial was detailed. Comparative verification of the four specific HTC formulas based on the test data shows that formula A corresponding to the zeroth-order polynomial always predicted stable HTC values, and at low flow rates, the four formulas predicted similar values for the HTC. When the flow rates increased, formulas B and D corresponding to the first- and third-order polynomials, respectively, predicted either too large or too small values of HTC, but formula C corresponding to the second-order polynomial presents relatively acceptable HTC values. However, this is only a subjective judgement because we cannot tell which formula gives more reasonable

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between…

4.

5.

values of the HTC between formulas A and C due to the limited information measured. However, as one of the advantages, formula C simultaneously presents a continuous distribution of temperature along the fracture wall. For a given test record, as the inlet and outlet temperatures are fixed, the predicted temperature on the fracture wall and the HTC of the specimen are primarily dependent on the absorbed heat Q. Therefore, how Q is obtained will have a significant impact on the prediction results. The main factors that affect the heat Q are the temperatures of the water at the inlet and outlet ends of the fracture. When the measured temperatures are directly treated as the average values at the inlet and the outlet, then the Q obtained with these temperatures is found suitable for formula A, but not suitable for formulas B, C, and D. This finding indicates that the temperatures at the inlet and the outlet used for Q should be consistent with the assumptions adopted in the derivation of the corresponding HTC formula. Our analytical method provides a tool to explore more realistic formulas for both the temperature of fracture water and the HTC. However, the existing test data are not sufficient to judge which formula is more reasonable. This is because both the new parameters using higherorder polynomials and the final verifications depend on additional experimental results. Further experimental work to measure the actual temperature distribution of water in the fracture will be of great value.

Acknowledgements The authors gratefully acknowledge the support of this work by the National Natural Science Foundation of China (Grant No. 41672252). Thanks a lot for the help by Prof. Guji Tian from the Institute of Physics and Mathematics of Chinese Academy of Sciences.

Appendix: Solution to Dirichlet Problem of Laplace Equation on a Semi-disk Green Function for Semi-disk in a Plane The Dirichlet problem of Laplace equation on a semi-disk reads, 8 Du ¼80 in X > < on L
u ¼ f1 ðxÞ > : : f2 ðxÞ

1

on L2 on L3

The materials below originated from Chapter 12 (Asmar 2005). We would like to apply the so-called method of images, which uses basic facts from plane geometry about

the circle and reflection, to derive the Green functions. Generally, the Green’s function for X  R2 is of the form,  1  Gðx; y; x0 ; y0 Þ ¼ ln ðx  x0 Þ2 þðy  y0 Þ2 þhðx; y; x0 ; y0 Þ 2 ð25Þ where h satisfies ( Dm ¼ 0   1 m ¼  ln ðx  x0 Þ2 þðy  y0 Þ2 2

in X on oX

ð26Þ

and then Gðx; y; x0 ; y0 Þ ¼ 0;

ð27Þ

on oX

Now we introduce a well-known transformation called the Steiner inversion, see Section 12.4 (Asmar 2005). When A is not the origin, the point A* is the inverse of the point A through the circle {(x, y) 2 R2: |x|2 ? |y|2 = R2}. By definition of the Steiner inversion, A* is the point on the ray from the origin O through A at a distance such that OA  OA ¼ R2 . In particular, if A = (x0,y0) [ {(x, y) 2 R2: 2 2 2 |x| ? |y| B R }, then 

x0 ; y0 ¼

R2 2

jx0 ; y0 j

ðx 0 ; y0 Þ ¼

R2 2

jx0 j þjy0 j2

ðx 0 ; y0 Þ

ð28Þ

Let " # 1 ðx  x0 Þ2 þðy  y0 Þ2 G1 ðx; y; x0 ; y0 Þ ¼ ln  2  2 2 x  x0 þ y  y0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ln x20 þ y20  R

ð29Þ

Then G1 ðx; y; x0 ; n y0 Þ;

on L1 ¼ ðx; yÞ 2 R2 : j xj2 þj yj2 ¼ R2 ; y [ 0

o

ð30Þ

The reflection of ðx0 ; y0 Þ about x-axis (y = 0) is ðx0 ; y0 Þ; let Gðx; y; x0 ; y0 Þ ¼ G1 ðx; y; x0 ; y0 Þ  G1 ðx; y; x0 ; y0 Þ " # " # 1 ðx  x 0 Þ2 þ ðy  y 0 Þ2 1 ðx  x 0 Þ2 þ ðy þ y 0 Þ 2 ¼ ln  2  2  ln  2  2 2 2 x  x þ y  y x  x þ y þ y 0

0

0

0

ð31Þ then Gðx; y; x0 ; y0 Þjy¼0 ¼ 0 which, by virtue of (29), implies Gðx; y; x0 ; y0 Þ ¼ 0; on oX ð32Þ n o 2 2 and with X ¼ ðx; yÞ 2 R2 : j xj þj yj \R2 ; y [ 0 ; therefore Gðx; y; x0 ; y0 Þ is the unique Green’s function

123

B. Bai et al.

determined by the region X. The uniqueness of Green’s function is from Theorem 3, Section 12.3, in the work by Asmar (2005). We separate the boundary of X into two parts: One is a n o semicircle ðx; yÞ 2 R2 : j xj2 þj yj2 ¼ R2 ; y [ 0 ; and the n o other is ðx; yÞ 2 R2 : j xj2 þj yj2  R2 ; y ¼ 0 : When we deal with the problem in a semicircle, it is convenient to express Green’s function in polar coordinates. Let 8 A ¼ ðx0 ; y0 Þ ¼ ðr0 cos /;  sin /Þ > > < 2 2 R R ð33Þ A ¼ cos /; sin / > r0 r0 > : P ¼ ðr cos h; r sin hÞ By (29) G1 ðx; y; x0 ; y0 Þ ¼ G1 ðr; h; r0 ; /Þ  1 r 2 þ r02  2rr0 cosðh  /Þ ¼ ln R2 2 2 2 r r0 þ R4  2rr0 R2 cosðh  /Þ Similarly G1 ðx; y; x0 ; y0 Þ ¼ G1 ðr; h; r0 ; /Þ  1 r 2 þ r02  2rr0 cosðh þ /Þ 2 ¼ ln R 2 2 2 r r0 þ R4  2rr0 R2 cosðh þ /Þ Therefore Gðr; h; r0 ; /Þ ¼ G1 ðr; h; r0 ; /Þ  G1 ðr; h; r0 ; /Þ  2 1 r þ r02  2rr0 cosðh  /Þ ¼ ln 2 2 2 r r0 þ R4  2rr0 R2 cosðh  /Þ  2 1 r þ r02  2rr0 cosðh þ /Þ  ln 2 2 2 r r0 þ R4  2rr0 R2 cosðh þ /Þ 

Now we turn to the case ðx; yÞ 2 R2 : j xj  R; y ¼ 0 ; and it follows from (28) that 8 h  i2 h  i2 2 2 2 2 2 2 2 2 > > x x þ y R x y x þ y R y 0 0   0 0 0 0 > >  2  2 > þ  2 2  2 2 > < x  x0 þ y  y0 ¼ x0 þ y20 x0 þ y20 h  i h i2 2 2  2 > > x x20 þ y20 R2 x0 y x20 þ y20 þR2 y0 >   > 2 2 > x  x þ y þ y ¼ > þ  2 2  2 2 : 0 0 x0 þ y20 x0 þ y20

Inserting the equalities above into (31), we have " # 1 ðx  x0 Þ2 þðy  y0 Þ2 Gðx; y;x0 ;y0 Þ ¼ ln   2   2 2 x x20 þ y20  R2 x0 þ y x20 þ y20  R2 y0 " # 1 ðx  x0 Þ2 þðy þ y0 Þ2  ln   :     2 2 2 x x2 þ y2  R2 x0 þ y x2 þ y2 þ R2 y0 0

0

0

We restate Theorem 2, Section 12.3 (Asmar 2005), as follows. Theorem 1 Suppose that u is harmonic in X and continuous on its boundary C. Then for ðx; yÞ 2 X, we have Z 1 oG uðx0 ; y0 ; x; yÞ ðx0 ; y0 ; x; yÞdsðx0 ; y0 Þ uðx; yÞ ¼ 2p C on ð36Þ where G is the Green’s function for X: Here we can have a different version from the original one of Asmar (2005) because the fact Gðx; y; x0 ; y0 Þ ¼ Gðx0 ; y0 ; x; yÞ is employed. In general, if we write uðx; yÞ ¼ uðr0 cos /; r0 sin /Þ; the outward normal vector n ¼ n o ðcos /; sin /Þ on the circle ðx; yÞ 2 R2 : j xj2 þj yj2 ¼ R2 . Therefore ou   ux ; uy  n ¼ ux cos / þ uy sin / on o ¼ ½uðr0 cos /; r0 sin /Þ or0 Now we calculate ou 0 ; y0 ; x; yÞ on ðxo n 2 2 2 2 ðx; yÞ 2 R : j xj þj yj ¼ R as

the

circle

where ds ¼ Rd/ is used. Calculations show that  1 o r 2 þ r02  2rr0 cosðh  /Þ ln R2 2 2 j 2 or0 r r0 þ R4  2rr0 R2 cosðh  /Þ r0 ¼R 1 R2  r 2 ¼ 2 R r þ R2  2rR cosðh  /Þ and similarly  1 o r 2 þ r02  2rr0 cosðh þ /Þ 2 ln R 2 2 j 2 or0 r r0 þ R4  2rr0 R2 cosðh þ /Þ r0 ¼R 2 2 1 R r ¼ 2 2 R r þ R  2rR cosðh þ /Þ By Theorem 1, we have Z jx0 j2 þjy0 j2 ¼R2 ;y [ 0  Z p

¼

uð/Þ

0

uðx0 ; y0 Þ

oG ðx0 ; y0 ; x; yÞdsðx0 ; y0 Þ on

R2  r2 R2  r 2  d/ r 2 þ R2  2rR cosðh  /Þ r2 þ R2  2rR cosðh þ /Þ

0

ð37Þ

ð35Þ Because 0  h; /  p; we see that

123

on

ou o ðx0 ; y0 ; x; yÞdsðx0 ; y0 Þ ¼ ½Gðr; h; r0 ; /Þjr0 ¼R Rd/ on or0

ð34Þ

Analytical Solution for Arbitrary Boundary Functions

An Analytical Method for Determining the Convection Heat Transfer Coefficient Between…

References

R2  r 2 R2  r 2  0 r 2 þ R2  2rR cosðh  /Þ r 2 þ R2  2rR cosðh þ /Þ

ð38Þ Since the outward normal vector on

n

ðx; yÞ 2 R2 : j xj2 þ

j yj2  R2 ; y ¼ 0g is ð0; 1Þ; then oG oG ¼ j j on y0 ¼0 oy0 y0 ¼0 and calculation shows by (35) that oG oG jy0 ¼0 ¼  j on oy0 y0 ¼0 2y 2yR2 ¼  ðx  x0 Þ2 þy2 ðxx0  R2 Þ2 þðyx0 Þ2 Notice that  x2    R2  2xx0 þ x2 þ y2 02  x20  2xx0 þ x2 þ y2 R    2 x2 þ y2 2 ¼ R  x0 1  0 R2 We see that oG j 0 on y0 ¼0 which, together with (38), yields oG j 0 on oX and, from (36) uðx; yÞ 0

if ujoX 0

By Theorem 1, we have Z

oG ðx0 ; y0 ; x; yÞdsðx0 ; y0 Þ on " # Z 0 2y 2yR2 ¼ f 1 ð xÞ  dx0 ðx  x0 Þ2 þy2 ðxx0  R2 Þ2 þðyx0 Þ2 R " # Z R 2y 2yR2 þ f 2 ð xÞ  dx0 ðx  x0 Þ2 þy2 ðxx0  R2 Þ2 þðyx0 Þ2 0

fjx0 j2 þjy0 j2 ¼R2 ;y¼0g

u ðx 0 ; y 0 Þ

ð39Þ

Asmar NH (2005) Partial Differential equations with fourier series and boundary value problems, 2nd edn. ISBN:0-13-148096-0 Bai B, He Y, Li X, Hu S (2016) Local heat transfer characteristics of water flowing through a single fracture within a cylindrical granite specimen. Environ Earth Sci 75:1460. doi:10.1007/ s12665-016-6249-2 Bai B, He Y, Li X, Li J, Huang X, Zhu J (2017) Experimental and analytical study of the overall heat transfer coefficient of water flowing through a single fracture in a granite core. Appl Therm Eng 116:79–90 Huang X, Zhu J, Li J, Bai B, Zhang G (2016) Fluid friction and heat transfer through a single rough fracture in granitic rock under confining pressure. Int Commun Heat Mass Transf 1(1):111–123 ISRM (1978) Suggested methods for determining tensile strength of rock materials part 2: suggested Method for determining indirect tensile strength by the Brazil Test. Int J Rock Mech Min Sci 15:99–103 Lu W, Xiang YY (2012) Experiments and sensitivity analyses for heat transfer in a meter-scale regularly fractured granite model with water flow. J Zhejiang Univ Sci A 13(12):958–968 Mohais R, Xu C, Dowd PA (2011) Fluid flow and heat transfer within a single horizontal fracture in an enhanced geothermal system. ASME J Heat Transf 133:112603-1–112603-8 Ogino F, Yamamura M, Fukuda T (1999) Heat transfer from hot dry rock to water flowing through a circular fracture. Geothermics 28(1):21–44 Xu RN, Zhang L, Zhang FZ, Jiang PX (2015) A review on heat transfer and energy conversion in the enhanced geothermal systems with water/CO2 as working fluid. Int J Energy Res 39(13):1722–1741 Zhang G, Zhu J, Li J, Wan Q (2015) The analytical solution of the water-rock heat transfer coefficient and sensitivity analyses of parameters. Proceedings World Geothermal Congress, Melbourne, Australia Zhao J (1987) Experimental studies of the hydro-thermo-mechanical behavior of joints in granite. Imperial College, London Zhao J (1992) Analytical and experimental studies of heat convection by water flow in rock fractures. In: Tillerson JRWW (ed) Proceeding of 33rd US Symposium on Rock Mechanics. Balkema, Rotterdam Zhao J (1999) Experimental study of flow-rock heat transfer in rock fractures. Chin J Rock Mech Eng 18(02):1–5 Zhao ZH (2014) On the heat transfer coefficient between rock fracture walls and flowing fluid. Comput Geotech 59:105–111 Zhao J, Brown ET (1992) Hydro-thermo-mechanical properties of joints in the carnmenellis granite. Q J Eng Geol 25(4):279–290 Zhao J, Tso CP (1993) Heat-transfer by water-flow in rock fractures and the application to hot dry rock geothermal systems. Int J Rock Mech Min Sci Geomech Abstr 30(6):633–641

From (36), (37), and (39), we have the explicit solution to the Dirichlet problem of Laplace equation on a semi-disk in plane uðx; yÞ ¼

 Z 1 p R2  r 2 R2  r 2  d/ uð/Þ 2 r þ R2  2rR cosðh  /Þ r 2 þ R2  2rR cosðh þ /Þ 2p 0 " # Z 0 2y 2yR2 þ f1 ðxÞ  dx0 ðx  x0 Þ2 þy2 ðxx0  R2 Þ2 þðyx0 Þ2 R " # Z R 2y 2yR2 þ f2 ð xÞ  dx0 : ðx  x0 Þ2 þy2 ðxx0  R2 Þ2 þðyx0 Þ2 0

ð40Þ

123