Using pygfunction in Ground Source Heat Pump simulation

[TobiasEnergyMachines]

Hi,
Thank you for building a great package.

I am attempting to use pygfunction to simulate a Ground Source Heat Pump (GSHP) system with a combined heat pump and chiller connected to a borehole field (BTES).

There seems to be a circular dependency when combining a BTES simulation based on g-functions with load aggregation together with a heat pump simulation based on (compressor) polynomials. Because the BTES model needs a load (Q) to compute the temperature at outlet and the heat pump model needs the temperature at inlet to compute the source load.

Do you have any ideas or have any good reference of how to overcome this circular dependency?

Thank you,
Tobias

[j-c-cook]

Dr. Jeffrey Spitler might have recently discussed this topic in a self-funded paper, possibly recognizing the need for addressing it in the context of incorporating a model into ghedt (his rebranded version is GHEDesigner). The abstract highlights the necessity for a validated approach using typical manufacturers' data sheets. This paper could be a valuable starting point for further exploration, and I've added it to my reading list due to my own interest in the topic.

I'm uncertain if @MassimoCimmino considers integrating this work within the package's scope. If he does, implementing an STS g-function could be a prerequisite step.

[MassimoCimmino]

@TobiasEnergyMachines You are correct that the heat pump simulation and the bore field simulation are coupled. This can be solved iteratively at each time step. Here is a sample simulation block that can replace the one in the fluid_temperature_multiple_boreholes example (in https://github.com/MassimoCimmino/pygfunction/blob/78fb638264b27ef84569bbdff638a7f7e80d4336/examples/fluid_temperature_multiple_boreholes.py). The code can be adapted for other heat pump modelling approaches.

    # -------------------------------------------------------------------------
    # Simulation
    # -------------------------------------------------------------------------

    # COP of the heat pump
    # These are not recommended or typical values and are simply here to
    # demonstrate how to integrate temperature-dependent COPs into the
    # simulation. The COP functions can be functions of more than one variable
    # by writing `COP = lambda var1, var2, var3: [...]`.

    # Cooling COP (-)
    COP_C = lambda T: 3.86
    # Heating COP (-)
    COP_H = lambda T: 4.00 + 0.065 * T

    # Evaluate building heating demand (positive for heating)
    Q_building = nBoreholes*synthetic_load(time/3600.)

    # Threshold for the load to be considered non-zero (this avoids errors
    # such as division by zero), in Watts
    Q_small = 10.0

    T_b = np.zeros(Nt)
    T_f_in = np.zeros(Nt)
    T_f_out = np.zeros(Nt)
    Q_ground = np.zeros(Nt)
    for i, (t, Q_i) in enumerate(zip(time, Q_building)):
        # Increment time step by (1)
        LoadAgg.next_time_step(t)

        # Iterate on the outlet fluid temperature
        T_f_out[i] = T_g
        error = np.inf
        while error > 0.005:
            T_f_out_last = T_f_out[i]

            # Apply the corresponding COP for cooling or heating
            if Q_i > 0.:
                COP = COP_H(T_f_out_last)
                Q_ground[i] = Q_i * (COP - 1) / COP
            else:
                COP = COP_C(T_f_out_last)
                Q_ground[i] = Q_i * (COP + 1) / COP

            # Apply current load (in watts per meter of borehole)
            Q_b = Q_ground[i]/nBoreholes
            LoadAgg.set_current_load(Q_b/H)
    
            # Evaluate borehole wall temperature
            deltaT_b = LoadAgg.temporal_superposition()
            T_b[i] = T_g - deltaT_b
    
            # Evaluate inlet fluid temperature (all boreholes are the same)
            T_f_in[i] = network.get_network_inlet_temperature(
                    Q_ground[i], T_b[i], m_flow_network, cp_f, nSegments=1)
    
            # Evaluate outlet fluid temperature
            T_f_out[i] = network.get_network_outlet_temperature(
                    T_f_in[i],  T_b[i], m_flow_network, cp_f, nSegments=1)

            error = np.abs(T_f_out[i] - T_f_out_last)

[wouterpeere]

Hi @TobiasEnergyMachines

I experienced kind of the same challenge with some of my thesis students a couple of years ago in the context of active/passive cooling with geothermal boreholes. They solved it by doing some iterations, like @MassimoCimmino suggested. It is implemented in GHEtool in this example.

Note that in contrast with GHEDesigner, GHEtool does not (yet) account for the short-term effects, which can lead to overestimation of the fluid temperatures, when using hourly simulated values.