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