Issue36 coaxial pipes

[j-c-cook]

In progress

[j-c-cook]

With my most recent commit (9022abc) I have modified the function that calculates the convective heat transfer coefficient in a circular pipe with forced flow. After some discussion with J.D. Spitler, he said that for design purposes it can be important to have smooth transitions from laminar to turbulent, which can be provided via interpolation. Cengel and Ghajar (2015, pg. 476) state that considering full turbulence in piping networks for any Re > 4000 is conservative. I have incorporated a "transition" 1D interpolation between the laminar (Nu=3.66) and the Nu(Re=4000)) while 2300 < Re < 4000.

By varying the mass flow rate and holding all other values constant, I plotted the Reynolds vs. the Nusselt number. jump.pdf plots the results of what the function was, and shows that there was a sharp jump from 3.66 to about 20 for the Nusselt number when the value of the Reynolds number became greater than 2300. The new function with interpolation from 2300 < Re < 4000 creates no-jump.pdf.

[MassimoCimmino]

The range 2300 < Re < 4000 is also found elsewhere in the literature, e.g. Gnielinski (2013).

Gnielinski, V. (2013). On heat transfer in tubes. International Journal of Heat and Mass Transfer, 63, 134-140.

[j-c-cook]

Commit 843dd7e provides no functionality difference. The linear interpolation is no longer done by scipy, but uses Equations 16 and 17 from Gnielinski (2013). Documentation is enhanced for Gnielinksi() and convective_heat_transfer_coefficient_circular_pipe(); references are included.

[j-c-cook]

Hellstrom's inner and outer Nusselt numbers are used for laminar flow. The transition from 2300 < Re < 4000 makes use of the linear interpolation from Gnielinski (2013) for both inner and outer. Cengel and Ghajar (2015) state that inner and outer Nusselt numbers are equivalent for turbulent flow, and that Gnielinski (1975) can be used.

A plot of Reynolds vs. Nusselt is attached.

Re_vs_Nu_annulus.pdf

[j-c-cook]

@MassimoCimmino Could you take a look at this and provide feedback if you'd like something changed?

I still need to do some further verification before its finalized, but the methodology is there. Now would be a good time to know if you'd like to see something different. Whenever you get the chance, thanks.

[MassimoCimmino]

I like the idea of the new borehole thermal resistance function. Here are my first thoughts:

• The matrix of thermal resistances SingleCoaxialPipe._Rd needs to be evaluated in __init__ to enable all functionalities of the pipe classes, which are needed for the 'MIFT' boundary condition in gFunction and for the Network class. The effective Rb (uniform temperature) could then be evaluated from the existing borehole_thermal_resistance() (This function is a generalization of the Hellström method to N pipes).
• On that same topic, there should be a way to obtain an expression for the uniform heat flow resistance from the R_delta matrix so that the new thermal resistance calculations are compatible with any pipe object. I will look into this.
• There are implications in implementing the 'UHTR' and 'AVG' effective thermal resistances. The effective borehole wall temperature calculated using the 'MIFT' relies on the 'UBWT' resistance to find an effective borehole wall temperature (as defined in [1]). This is not a huge problem as long as the Network class is used to evaluate fluid temperatures during simulations (even for simple fields with parallel-connected boreholes). We need to make sure the default parameters do not break the implementation of the gFunction or Network classes.
• While it is not ideal, the existing pipe classes do not calculate R_f internally but instead take it as an input. We could consider opening a new issue for this change. This also provides opportunities to implement temperature dependent physical properties in the pipe and Network classes using the new media module.
• Using the same r_in and r_out arrays as in other classes may be preferable to avoid implementing a new visualize_pipes() method (I know the variable definition are from my commit).
• A new _check_geometry() method is required.

---

[1] Cimmino, M. (2019). Semi-Analytical Method for g-Function Calculation of bore fields with series-and parallel-connected boreholes. Science and Technology for the Built Environment, 25(8), 1007-1022.

[j-c-cook]

• The thermal_resistances_coaxial function takes in one pipe resistance R_p value, how will we account for the inner and outer pipe having a different resistance?
• Given that we now have pipes with different radii, how will we handle that the existing visualize_pipes() only considers an r_in and r_out?
• Do you want a new _check_geometry() method or to modify the existing?

[MassimoCimmino]

• I did not think of the different R_p values. I cannot think of an elegant way to implement it in a convenient interface. It might be a good idea to move ahead with opening a new issue and handle the calculation of film resistances at the initialization of pipe objects instead of having R_fp as an input. It will make the interface to the different pipes consistent, and v2.0.0 is the best time to change the interface in a non-backwards-compatible way.
• I got confused here and misremembered that r_in and r_out were arrays (they are not). The current visualize_pipes() method can be adapted to cover both cases (with r_in and r_out as arrays in coaxial pipes). The code below creates r_in and r_out arrays when floats or scalars are provided:

    def visualize_pipes(self):
        """
        Plot the cross-section view of the borehole.

        Returns
        -------
        fig : figure
            Figure object (matplotlib).

        """
        import matplotlib.pyplot as plt
        from matplotlib.ticker import AutoMinorLocator

        # Initialize figure
        LW = .5    # Line width
        FS = 12.    # Font size

        plt.rc('figure', figsize=(80.0/25.4, 80.0/25.4))
        fig = plt.figure()
        ax = fig.add_subplot(111)

        # Borehole wall outline
        borewall = plt.Circle((0., 0.), radius=self.b.r_b,
                              fill=False, linestyle='--', linewidth=LW)
        ax.add_patch(borewall)

        # Pipes
        r_in = np.atleast_1d(self.r_in)
        if len(r_in) == 1: r_in = np.repeat(r_in, 2*nPipes)
        r_out = np.atleast_1d(self.r_out)
        if len(r_out ) == 1: r_out = np.repeat(r_out , 2*nPipes)
        for i in range(self.nPipes):
            # Coordinates of pipes
            (x_in, y_in) = self.pos[i]
            (x_out, y_out) = self.pos[i + self.nPipes]

            # Pipe outline (inlet)
            pipe_in_in = plt.Circle((x_in, y_in), radius=r_in[i],
                                    fill=False, linestyle='-', linewidth=LW)
            pipe_in_out = plt.Circle((x_in, y_in), radius=r_out[i],
                                     fill=False, linestyle='-', linewidth=LW)
            ax.text(x_in, y_in, i + 1,
                    ha="center", va="center", size=FS)

            # Pipe outline (outlet)
            pipe_out_in = plt.Circle((x_out, y_out), radius=r_in[i + self.nPipes],
                                     fill=False, linestyle='-', linewidth=LW)
            pipe_out_out = plt.Circle((x_out, y_out), radius=r_out[i + self.nPipes],
                                      fill=False, linestyle='-', linewidth=LW)
            ax.text(x_out, y_out, i + self.nPipes + 1,
                    ha="center", va="center", size=FS)

            ax.add_patch(pipe_in_in)
            ax.add_patch(pipe_in_out)
            ax.add_patch(pipe_out_in)
            ax.add_patch(pipe_out_out)

        # Configure figure axes
        ax.set_xlabel('x (m)')
        ax.set_ylabel('y (m)')
        plt.axis('equal')
        ax.xaxis.set_minor_locator(AutoMinorLocator())
        ax.yaxis.set_minor_locator(AutoMinorLocator())
        plt.tight_layout()

        return fig

• A new _check_geometry() is required, since the current one checks for collision between pipes. It also checks the other inputs.

[j-c-cook]

@MassimoCimmino I have modified custom_borehole.py to include Double U-Tubes in series and parallel. Once we have worked out the Coaxial interface, we could add Coaxial to the example as well. The custom borehole example would then serve the purpose of seeing how to define any tube arrangement available in pygfunction all in one condensed location.

[MassimoCimmino]

@j-c-cook I started the review on my branch. I will let you know if anything comes up.

Remaining checks :

• convective_heat_transfer_coefficient_concentric_annulus()
• custom_borehole.py
• documentation
• cross-validation with the other tools

[j-c-cook]

A Coaxial_Comparison.zip has been created to help to overcome the apprehension to merge this branch due to the comparison results to other tools (GLHEDT and GLHEPRO). This apprehension was introduced due to the vast differences in effective borehole resistances shown in the plot here.

The attached file ranges through a list of volumetric flow rates and computes the effective borehole thermal resistances for pygfunction and GLHEDT. At each volumetric flow rate, the Reynolds, effective resistances and percentage difference between the two methods are printed to the console.

Reynolds	GLHEDT	pygfunction	%diff
8549	0.1024	0.1546	40.7
7694	0.1029	0.1651	46.4
6839	0.1035	0.1786	53.2
5984	0.1044	0.1967	61.3
5129	0.1055	0.2215	70.9
4274	0.1072	0.2567	82.1
3419	0.1116	0.3075	93.5
2565	0.1298	0.3807	98.3
2479	0.1350	0.3898	97.1
2394	0.1424	0.3996	94.9
2308	0.1536	0.4108	91.1
2223	0.1550	0.4253	93.2

For each volumetric flow rate, there are checks to compare differences in the local borehole thermal resistance, fluid-to-fluid thermal resistance calculation, and to see if the equation 3.68 provides any correction for short circuiting. None of the print statements are ever reached, because the boolean values are never true.

The method of GLHEDT (and GLHEPRO) are taken from Grundman (2016). The equation 3.68 is from Javed and Spitler (2016), which is first mentioned in Claesson and Hellström (2011). Claesson and Hellström reference the original derivation of the short circuiting correction, given in Eskilson and Claesson (1988). None of which discuss the validity of the method used by Grundman (2016) for concentric tube heat exchangers. Hellström (1996) presents the values for the thermal delta-circuit for a concentric tube (annular duct) heat exchanger. The same equations that apply to a single U-tube can be applied to by use of the delta network. However, the values of the delta network are different, notably the value of the resistance from the inner fluid to the outer wall is given as infinity by Hellström (1996).

---

What about GLHEDT seems off? The correction for the short circuiting (equation 3.68) is never greater than 1.0e-06. No literature broaches the subject of equation 3.68's validity for concentric tube heat exchangers. The validity of setting Ra = R12 (Grundman 2016) is not proven.

Why do the results of pygfunction make sense? We would expect the value of the short circuiting to be greater for concentric tube heat exchangers compared to u-tubes because the thermal resistance separating the two fluids flowing in opposite direction are much closer to one another.

What if GLHEDT and pygfunction are wrong? The comparison of single U-tubes in this plot validates pygfunctions implementation of Eskilson and Claesson (1998) delta resistance equations. The same functions that apply to the solution of two fluid channels with counter flow are used in the new Concentric object.

---

Grundmann, R. 2016. Improved Design Methods for Ground Heat Exchangers. M.S. Thesis, Oklahoma State University.

Javed S, Spitler JD. Calculation of borehole thermal resistance. In: Simon J. Reese, editor. Advances in ground-source heat pumps sytems. Woodhead Publishing; 2016. p. 63-95. http://dx.doi.org/10.1016/B978-0-08-100311-4.00003-0

Hellström, Göran. “Ground Heat Storage : Thermal Analyses of Duct Storage Systems.” PhD Doctoral Thesis. University of Lund; 1991.

Claesson, J., & Hellström, G. (2011). Multipole method to calculate borehole thermal resistances in a borehole heat exchanger. HVAC&R Research, 17(6), 895–911. https://doi.org/10.1080/10789669.2011.609927

Per Eskilson & Johan Claesson (1988) SIMULATION MODEL FOR THERMALLY INTERACTING HEAT EXTRACTION BOREHOLES, Numerical Heat Transfer, 13:2, 149-165, DOI: 10.1080/10407788808913609

[MassimoCimmino]

The differences in the script (Coaxial_comparison.py) are caused by a mismatch of units in the definition of the mass and volume flow rates.

Line 59:

        m_flow_borehole = V_flow_borehole / 1000. * fluid.rho

This suggests that V_flow_borehole has units of [L/s].

Lines 105-107:

        eta = H / (fluid.rho * fluid.cp * V_flow_borehole) * 1 / (2 * Rb) * \
              np.sqrt(1 + 4 * Rb / R_12)  # eq. (3.69)
        Rb_star_UBH = Rb * eta * (1 / np.tanh(eta))  # eq. (3.68)

The units of rho and cp are [kg/m3] and [J/kg-K], as given by the media module.

Introducing the following correction (divide V_flow_borehole / 1000.) makes the two methods equal:

        eta = H / (fluid.rho * fluid.cp * V_flow_borehole / 1000.) * 1 / (2 * Rb) * \
              np.sqrt(1 + 4 * Rb / R_12)  # eq. (3.69)
        Rb_star_UBH = Rb * eta * (1 / np.tanh(eta))  # eq. (3.68)

We get:

Reynolds	GLHEDT	pygfunction	%diff
Equation 3.68 is doing a short circuiting correction!!
8549	0.1546	0.1546	0.0
Equation 3.68 is doing a short circuiting correction!!
7694	0.1651	0.1651	0.0
Equation 3.68 is doing a short circuiting correction!!
6839	0.1786	0.1786	0.0
Equation 3.68 is doing a short circuiting correction!!
5984	0.1967	0.1967	0.0
Equation 3.68 is doing a short circuiting correction!!
5129	0.2215	0.2215	0.0
Equation 3.68 is doing a short circuiting correction!!
4274	0.2567	0.2567	0.0
Equation 3.68 is doing a short circuiting correction!!
3419	0.3075	0.3075	0.0
Equation 3.68 is doing a short circuiting correction!!
2565	0.3807	0.3807	0.0
Equation 3.68 is doing a short circuiting correction!!
2479	0.3898	0.3898	0.0
Equation 3.68 is doing a short circuiting correction!!
2394	0.3996	0.3996	0.0
Equation 3.68 is doing a short circuiting correction!!
2308	0.4108	0.4108	0.0
Equation 3.68 is doing a short circuiting correction!!
2223	0.4253	0.4253	0.0

[j-c-cook]

Good catch. I stand corrected. GLHEDT and pygfunction do provide the same results. I introduced a unit error.

Now, why do GLHEDT and pygfunction have a different result than GLHEPRO? Perhaps this is outside of the scope of this PR?

[MassimoCimmino]

A possible cause is differences in the calculation of the thermal resistances R_ff and R_fp. The effective borehole thermal resistance tends to decrease when R_ff increases. I will verify the implementation of the convective coefficient correlations again to make sure. Could GLHEPRO assume there is insulation between the two pipes?

[j-c-cook]

Thank you for your diligence. Progress is much appreciated.

I cannot speak for the inner workings of GLHEPRO, only the output. For example, see below:

[MassimoCimmino]

What is the length of the borehole in this case? These results match with pygfunction for H = 100.

    # -------------------------------------------------------------------------
    # Simulation parameters
    # -------------------------------------------------------------------------

    # Borehole dimensions
    D = 5.          # Borehole buried depth (m)
    H = 100.        # Borehole length (m)
    r_b = 0.075    # Borehole radius (m)

    # Pipe dimensions (all configurations)
    epsilon = 1.0e-6    # Pipe roughness (m)

    # Pipe dimensions (coaxial)
    r_in_in = 0.0221    # Inside pipe inner radius (m)
    r_in_out = 0.025    # Inside pipe outer radius (m)
    r_out_in = 0.0487   # Outer pipe inside radius (m)
    r_out_out = 0.055   # Outer pipe outside radius (m)
    # Vectors of inner and outer pipe radii
    # Note : The dimensions of the inlet pipe are the first elements of
    #        the vectors. In this example, the inlet pipe is the inside pipe.
    r_inner = np.array([r_in_in, r_out_in])     # Inner pipe radii (m)
    r_outer = np.array([r_in_out, r_out_out])   # Outer pip radii (m)

    # Ground properties
    k_s = 2.0           # Ground thermal conductivity (W/m.K)

    # Grout properties
    k_g = 1.0           # Grout thermal conductivity (W/m.K)

    # Pipe properties
    k_p = 0.4           # Pipe thermal conductivity (W/m.K)

    # Fluid properties
    fluid = gt.media.Fluid('Water', 0.)
    cp_f = fluid.cp     # Fluid specific isobaric heat capacity (J/kg.K)
    rho_f = fluid.rho   # Fluid density (kg/m3)
    mu_f = fluid.mu     # Fluid dynamic viscosity (kg/m.s)
    k_f = fluid.k       # Fluid thermal conductivity (W/m.K)
    # Total fluid mass flow rate per borehole (kg/s)
    m_flow_borehole = 0.2 / 1000. * rho_f

    # -------------------------------------------------------------------------
    # Initialize borehole model
    # -------------------------------------------------------------------------

    borehole = gt.boreholes.Borehole(H, D, r_b, x=0., y=0.)

    # -------------------------------------------------------------------------
    # Define a coaxial borehole
    # -------------------------------------------------------------------------

    # Pipe positions
    # Coaxial pipe (x, y)
    pos = (0., 0.)

    # Pipe thermal resistance
    # (the two pipes have the same thermal conductivity, k_p)
    # Inner pipe
    R_p_in = gt.pipes.conduction_thermal_resistance_circular_pipe(
        r_in_in, r_in_out, k_p)
    # Outer pipe
    R_p_out = gt.pipes.conduction_thermal_resistance_circular_pipe(
        r_out_in, r_out_out, k_p)

    # Fluid-to-fluid thermal resistance
    # Inner pipe
    h_f_in = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
        m_flow_borehole, r_in_in, mu_f, rho_f, k_f, cp_f, epsilon)
    R_f_in = 1.0 / (h_f_in * 2 * pi * r_in_in)
    # Outer pipe
    h_f_a_in, h_f_a_out = \
        gt.pipes.convective_heat_transfer_coefficient_concentric_annulus(
            m_flow_borehole, r_in_out, r_out_in, mu_f, rho_f, k_f, cp_f,
            epsilon)
    R_f_out_in = 1.0 / (h_f_a_in * 2 * pi * r_in_out)
    R_ff = R_f_in + R_p_in + R_f_out_in

    # Coaxial GHE in borehole
    R_f_out_out = 1.0 / (h_f_a_out * 2 * pi * r_out_in)
    R_fp = R_p_out + R_f_out_out
    Coaxial = gt.pipes.Coaxial(
        pos, r_inner, r_outer, borehole, k_s, k_g, R_ff, R_fp, J=2)

    # Evaluate and print the effective borehole thermal resistance
    R_b = gt.pipes.borehole_thermal_resistance(
        Coaxial, m_flow_borehole, cp_f)
    print('Coaxial tube Borehole thermal resistance: {0:.4f} m.K/W'.
          format(R_b))
    print('Convection coefficient at inside of inner pipe: {:.4f} m.K/W \n'
          'Convection coefficient at outside of inner pipe: {:.4f} m.K/W \n'
          'Convection coefficient at the outer pipe: {:.4f} m.K/W'.
          format(h_f_in, h_f_a_in, h_f_a_out))

    # Visualize the borehole geometry and save the figure
    fig_coaxial = Coaxial.visualize_pipes()
    fig_coaxial.savefig('coaxial-borehole.png')

Results:

Coaxial tube Borehole thermal resistance: 0.1929 m.K/W
Convection coefficient at inside of inner pipe: 618.7905 m.K/W 
Convection coefficient at outside of inner pipe: 71.9854 m.K/W 
Convection coefficient at the outer pipe: 57.0229 m.K/W

[j-c-cook]

Yes, the height is 100 m. That must've been the mistake all along. Here's an updated validation plot.

[MassimoCimmino]

Great!

I believe the remaining differences for double U-tubes are due to the call to borehole_thermal_resistance(pipe, m_flow_borehole, cp_f) which used m_flow_pipe instead of m_flow_borehole. I corrected this in 3a70afd.

I will wrap up the PR in the next few days.

[j-c-cook]

I believe GLHEPRO only has flow in parallel.

jcc:issue36_CoaxialPipes is now even with MassimoCimmino:issue36_CoaxialPipes.

[j-c-cook]

Yes.

Double U-tube (fixed).zip