Difference between BHResist and pygfunction

[wouterpeere]

Hi @MassimoCimmino,

I don't know if you are familiar with the GitHub package BHResist (https://github.com/NREL/BHResist), but they have implemented the closed-form solution for the multipole method (first order) for single, double and coaxial probes and they have compared it to multipole method in pygfunction (https://github.com/NREL/BHResist/blob/main/notebooks/Comparisons.ipynb). They found a good overlap between their first order model and your, by default, second order multipole.

However, when you do the comparison between their explicit method and a first order multipole in pygfunction, there is still a difference, although I would think it should be identical.

Do you have an idea what might be the cause of this difference?

Best,
Wouter

[MassimoCimmino]

I was not aware of BHResist.

I see a few possible causes for the (small) differences, ordered here based on chronology of computations to evaluate Rb:

1. Convection resistances: The correlations used for the evaluation of fluid resistances in circular pipes may not be the same. We see fifth decimal differences on Rf in the notebook you linked. The annular pipe coefficients are based on the same reference so they should be the same. Both should be compared in more detail to verify the causses for the differences.
2. Multipole method: Differences persist when using J=1 (first order multipole) in pygfunction. The explicit solutions used in BHResist should be exact but pygfunction uses the iteratice scheme from Hellström, which is expected to cause some (hopefully negligible) errors on the delta-circuit resistances. The evaluation of the delta-circuit resistances (for known Rf) should be compared. In my other project geothermsim, a direct solver (rather) than an iterative solver is used to evaluate delta-circuit resistances, so it could be added to the comparison.
3. Effective borehole thermal resistance: I see they used the Hellström correction Rb* = Rb * n * coth(n) to evaluate the effective thermal resistance. Both this and the method implemented in pygfunction should be equivalent and exact. The evaluation of the local borehole thermal resistance should be compared first between the two tools and then the evaluation of the effective borehole thermal resistance.

[wouterpeere]

Hi @MassimoCimmino , thank you for your response. I investigated somewhat further.

1. There are some differences here, for example BHResist has a transition zone between Re=2000-4000 (pyfunction 2300-4000), has a Nu_lam of 4.01 (pygfunction 3.66) and they both have different smoothing curve (Sigmoid function for BHResist and linear interpolation for pygfunction). The difference seems rather small (indeed fifth decimal), but for the sake of comparison, I will use constant values for the comparison of the other two points.
2. When I manually overwrite the Rf in BHResist to match the one from pygfunction, the local Rb is identical (up to 5 decimal places) (at least for the U-pipes, the coaxial I was not yet able to check).
3. With the same identical local resistances, the effective borehole resistance still differs between BHResist and pygfunction. I also added geothermsim (with J=1), but since I'm not that familiar with it, I have not changed the Rf manually. The results from pygfunction and geothermsim are very alike.

To me, it looks like there are some differences related to the calculation of the effective borehole thermal resistance.

[MassimoCimmino]

I was also working on this on my side. My notes below repeat your observations.

Important notes:

• BHResist uses UNIFORM_HEAT_FLUX as a default value to evaluate the effective thermal resistance. It should be overwritten with the parameter boundary_condition="UNIFORM_BOREHOLE_WALL_TEMP" when instanciating the object.
• pygfunction uses J=2 as the default multipole order. It should be overwritten with the parameter J=1 when instanciating the object.

Convective heat transfer coefficient (circular pipes)

In BHResist, I also found that:

1. In the laminar regime, the convection heat transfer coefficient is the mean of the uniform heat flux and uniform wall temperature (so Nu = 4.01, as you pointed out). In pygfunction : Nu = 3.66 (uniform wall temperature).
2. The smoothing function between uses a sigmoid from Re=2000 to Re=4000. In pygfunction : linear interpolation from Re=2300 to Re=4000.

Figure A. Convective heat transfer coefficient in circular pipes

Table A. Convective heat transfer coefficient in circular pipes

These two differences explain the differences below Re=4000. Above Re=4000, the small remaining differences are due to choices of correlations.

Internal thermal resistances (single U-tube borehole)

Here, I used R_f from BHResist in pygfunction (opposite to what you did) and evaluated R_b, R_a and R_b*. If I read the BHResist code correctly, R_f is evaluated every call so it cannot be overwritten.

Figure B. Local borehole thermal resistance R_b in a single U-tube borehole

Figure C. Local short-circuit thermal resistance R_a in a single U-tube borehole

Table B. Effective borehole thermal resistance R_b* in a single U-tube borehole

The differences are small when using boundary_condition="UNIFORM_BOREHOLE_WALL_TEMP". The differences disappear when using the same R_f for single U-tube boreholes.

I will do coaxial boreholes next.

[MassimoCimmino]

I may have found an issue with the evaluation of the Reynolds number in BHResist, which I reported in NREL/BHResist#9.

There might also be a mishandling of the hydraulic diameter in pygfunction:

pygfunction/pygfunction/pipes.py

Lines 3257 to 3274 in 443d866

	# Hydraulic diameter and radius for concentric tube annulus region
	D_h = 2 * (r_a_out - r_a_in)
	r_h = D_h / 2
	# Cross-sectional area of the annulus region
	A_c = pi * ((r_a_out ** 2) - (r_a_in ** 2))
	# Volume flow rate
	V_dot = np.abs(m_flow_pipe) / rho_f
	# Average velocity
	V = V_dot / A_c
	# Reynolds number
	Re = rho_f * V * D_h / mu_f
	# Prandtl number
	Pr = cp_f * mu_f / k_f
	# Ratio of radii (Grundmann, 2016)
	r_star = r_a_in / r_a_out
	# Darcy-Wiesbach friction factor
	fDarcy = fluid_friction_factor_circular_pipe(
	m_flow_pipe, r_h, mu_f, rho_f, epsilon)

The code above calls fluid_friction_factor_circular_pipe using r_h instead of r_in as an attempt to avoid the implementation of a dedicated function for annular pipes. However, this may lead to a mistake due to the use of r_h to evaluate the fluid velocity.

Additional differences come from different correlations in the turbulent regime (BHResist: Dittus-Boelter vs pygfunction: Gnielinski).

After fixing both issues and using Dittus-Boelter in pygfunction:

[MassimoCimmino]

Convective heat transfer coefficient (coaxial pipes)

Figure D. Convective heat transfer coefficients in annular channels

Table C. Convective heat transfer coefficients in coaxial boreholes

[wouterpeere]

@MassimoCimmino

Thank you for the evaluation! I missed the difference between the uniform borehole wall/uniform heat extraction setting.
Can you perhaps share your code for the evaluation? I want to compare some explicit second order degree multipole methods as well.

[MassimoCimmino]

Here is for single U-tube boreholes.

SingleUTube.ipynb

[wouterpeere]

Thanks!
I tried running the same for the DoubleU comparison and there seems to be quite a lot of differences there, even though pyfunction and pygfunction with the R_f from BHResist is (again) more or less identical.

test_double.py

[MassimoCimmino]

I think the pipe arrangements are different. In your script (lines 89-103):

double_bhr = Borehole()
double_bhr.init_double_u_borehole(
    borehole_diameter=r_b * 2,
    pipe_outer_diameter=r_out * 2,
    pipe_dimension_ratio=(2 * r_out) / (r_out - r_in),
    length=H,
    shank_space=D_s,
    pipe_conductivity=k_p,
    grout_conductivity=k_g,
    soil_conductivity=k_s,
    fluid_type="PROPYLENEGLYCOL",
    fluid_concentration=0.2,
    boundary_condition="UNIFORM_BOREHOLE_WALL_TEMP",
    pipe_inlet_arrangement="ADJACENT" # <-------- CHANGE HERE (from pipe_inlet_arrangement="DIAGONAL").
)

This gives:

------  -------  --------  -----------  ------------  -------------------------------  ------------
m_flow  Re       BHResist  pygfunction  % difference  pygfunction (R_f from BHResist)  % difference
0.050   533.32   1.05575   1.05533      0.04          1.05575                          0.00
0.056   598.40   0.94094   0.94057      0.04          0.94094                          0.00
0.063   671.41   0.83863   0.83830      0.04          0.83863                          0.00
0.071   753.34   0.74746   0.74718      0.04          0.74746                          0.00
0.079   845.26   0.66627   0.66604      0.03          0.66627                          0.00
0.089   948.40   0.59401   0.59387      0.02          0.59401                          0.00
0.100   1064.12  0.52983   0.52979      0.01          0.52983                          0.00
0.112   1193.96  0.47298   0.47308      -0.02         0.47298                          0.00
0.126   1339.64  0.42283   0.42312      -0.07         0.42283                          0.00
0.141   1503.11  0.37884   0.37938      -0.14         0.37884                          0.00
0.158   1686.51  0.34053   0.34136      -0.24         0.34053                          0.00
0.177   1892.30  0.30745   0.30859      -0.37         0.30745                          0.00
0.199   2123.19  0.27756   0.28060      -1.09         0.27756                          0.00
0.223   2382.26  0.24990   0.24988      0.01          0.24990                          0.00
0.251   2672.94  0.22275   0.22150      0.57          0.22275                          0.00
0.281   2999.09  0.19944   0.19993      -0.25         0.19944                          0.00
0.315   3365.03  0.18129   0.18204      -0.41         0.18129                          0.00
0.354   3775.63  0.16676   0.16707      -0.19         0.16676                          -0.00
0.397   4236.33  0.15464   0.15469      -0.03         0.15464                          0.00
0.446   4753.24  0.14445   0.14449      -0.03         0.14445                          0.00
0.500   5333.22  0.13601   0.13605      -0.03         0.13601                          0.00
------  -------  --------  -----------  ------------  -------------------------------  ------------

[wouterpeere]

Interesting. Is there any way in pygfunction to differentiate between diagonal and adjacent configuration?

[wouterpeere]

By trial-and-error, I found out this is determined by the position array, but I couldn't find anything in the documentation about this. Perhaps it could be a good idea to clarify this somewhere.