TNO Intern

Commit 2ebb5e9c authored by Hen Brett's avatar Hen Brett 🐔
Browse files

Ensuring the heat pump calculation is the same as the Java code

parent 7afe7194

src/pythermogis/workflow/utc/doublet.py

+11 −12
Original line number Diff line number Diff line
@@ -49,22 +49,14 @@ class DoubletOutput: 
    injection_temp: float



def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) -> DoubletOutput | None:

    # determine initial well distance

    well_distance = (

        (props.optim_dist_well_dist_min + props.optim_dist_well_dist_max) / 2

        if props.optim_well_dist

        else props.default_well_distance

    )





    injection_temperature = (

        max(input.temperature - props.max_cooling_temp_range,

            props.hp_minimum_injection_temperature)

        if props.use_heat_pump

        else max(input.temperature - props.max_cooling_temp_range,

                 props.dh_return_temp)

    )



    if props.use_heat_pump and props.calculate_cop and not props.hp_application_mode:

    if props.use_heat_pump:

        injection_temperature = calculate_injection_temp_with_heat_pump(

            input.temperature,

            props.hp_direct_heat_input_temp,
@@ -73,7 +65,11 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
            props.max_cooling_temp_range,

            props.hp_minimum_injection_temperature,

        )

    else:

        injection_temperature = max(input.temperature - props.max_cooling_temp_range,

            props.dh_return_temp)



    # calculate maximum pressure

    drawdown_pressure = calculate_max_pressure(

        props,

        input,
@@ -85,6 +81,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
    if drawdown_pressure == 0:

        return None



    # cooling temperature and well distance optimization

    if props.optim_well_dist:

        end_temperature_p = (

                props.optim_dist_cooling_fraction *
@@ -107,12 +104,14 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
            injection_temperature,

        )



    # stimulation capex

    stimulation_capex = (

        0.0

        if (not props.use_stimulation or input.transmissivity_with_ntg > props.stim_kh_max)

        else props.stimulation_capex

    )



    # pressure optimizer

    pressure_results = optimize_pressure(

        props=props,

        input=input,
@@ -125,6 +124,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
    if pressure_results is None:

        return None



    # everything underneath here is lightning

    heat_power_per_doublet = pressure_results.heat_power_per_doublet

    flowrate = pressure_results.flowrate

    discounted_heat_produced_p = pressure_results.discounted_heat_produced
@@ -181,7 +181,6 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
            * (1 - props.tax_rate)

    )

    opex_first_prod_year = total_opex_ts[props.drilling_time]

    hp_cop = 3.0



    return DoubletOutput(

        power=heat_power_per_doublet,
@@ -194,7 +193,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) 
        npv=npv,

        hprod=discounted_heat_produced_p,

        cop=cop,

        cophp=hp_cop,

        cophp=pressure_results.heat_pump_cop,

        pres=drawdown_pressure / 1e5,

        flow=flowrate,

        welld=well_distance,

src/pythermogis/workflow/utc/flow.py

+6 −2
Original line number Diff line number Diff line
@@ -10,6 +10,7 @@ from utils.timer import print_time 
class VolumetricFlowResults:

    hp_cop: float

    hp_added_power: float

    hp_elec_consumption: float

    heat_power_per_doublet: float

    cop: float

    flowrate: float
@@ -27,8 +28,10 @@ def calculate_volumetric_flow( 
):

    STEPS = [0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, -6]

    BAR_SI = 1e5  # Units.BAR_SI

    hp_cop = 3.0



    hp_cop = 0.0

    hp_added_power = 0.0

    hp_elec_consumption = 0.0



    d1d_results = None

    iter_index = 0
@@ -90,13 +93,14 @@ def calculate_volumetric_flow( 
            injection_temp,

            flowrate

        )



        hp_cop = hp_results.hp_cop

        hp_added_power = hp_results.hp_added_power

        hp_elec_consumption = hp_results.hp_elec_consumption



    return VolumetricFlowResults(

        hp_cop,

        hp_added_power,

        hp_elec_consumption,

        heat_power_per_doublet,

        cop,

        flowrate,

src/pythermogis/workflow/utc/heatpump.py

+18 −15
Original line number Diff line number Diff line
@@ -7,6 +7,7 @@ from pythermogis.workflow.utc.water import get_hydrostatic_pressure, get_salinit 
class HeatPumpPerformanceResults:

    hp_cop: float

    hp_added_power: float

    hp_elec_consumption: float





def calculate_heat_pump_performance(
@@ -35,28 +36,31 @@ def calculate_heat_pump_performance( 
        hydrostatic_pressure, props.surface_temperature, salinity

    )



    cp_water = heat_capacity(

        hydrostatic_pressure, props.surface_temperature, salinity

    cp_water = heat_capacity(hydrostatic_pressure, props.surface_temperature, salinity)



    # calculate the power used by the heat pump

    waste_temp = min(heat_pump_start_temp, props.dh_return_temp)

    delta_temp = max(0.0, waste_temp - injection_temp)

    hp_added_power = (flowrate / 3600.0) * delta_temp * rho_water * cp_water * 1e-6

    if hp_added_power == 0:

        return HeatPumpPerformanceResults(

            hp_cop=0,

            hp_added_power=0,

            hp_elec_consumption=0,

        )



    if props.hp_calc_cop:

    # calculate heat pump coefficient of performance

    Tout = props.hp_direct_heat_input_temp

    Tin = min(injection_temp, props.dh_return_temp)

    hp_cop = get_cop_carnot(ETA_CARNOT, Tout, Tin)

    else:

        hp_cop = props.hp_capex



    waste_temp = min(heat_pump_start_temp, props.dh_return_temp)

    delta_temp = max(0.0, waste_temp - injection_temp)



    hp_added_power = (flowrate / 3600.0) * delta_temp * rho_water * cp_water * 1e-6



    if hp_added_power == 0:

        hp_cop = 0.0

    # calculate hp_elec_consumption

    hp_elec_consumption = hp_added_power / (hp_cop -1)



    return HeatPumpPerformanceResults(

        hp_cop=hp_cop,

        hp_added_power=hp_added_power,

        hp_elec_consumption=hp_elec_consumption

    )
@@ -64,11 +68,10 @@ def calculate_heat_pump_start_temp(props, production_temp: float, injection_temp 
    """

    Python version of calculateHeatPumpStartTemp(...)

    """



    delta_temp_geothermal = production_temp - injection_temp

    delta_temp_for_hex = 0.0



    if props.hp_application_mode and production_temp > props.dh_return_temp:

    if production_temp > props.dh_return_temp:

        delta_temp_for_hex = production_temp - props.dh_return_temp

        delta_temp_for_hex = max(0.0, min(delta_temp_for_hex, delta_temp_geothermal))

src/pythermogis/workflow/utc/pressure.py

+2 −1
Original line number Diff line number Diff line
@@ -74,6 +74,7 @@ class PressureOptimizerResults: 
    drawdown_pressure: float

    flowrate: float

    heat_power_per_doublet: float

    heat_pump_cop: float

    cop: float

    utc: float

    sum_capex: float
@@ -106,7 +107,6 @@ def optimize_pressure( 
    if flow_results is None:

        return None





    if flow_results.production_temp > injection_temp:



        iter_count = 0
@@ -222,6 +222,7 @@ def optimize_pressure( 
        drawdown_pressure=drawdown_pressure,

        flowrate=flow_results.flowrate,

        heat_power_per_doublet=flow_results.heat_power_per_doublet,

        heat_pump_cop=flow_results.hp_cop,

        cop=cop,

        utc=econ.utc.utc,

        sum_capex=econ.capex.sum_capex,

tests/utc/test_doublet.py

+2 −2
Original line number Diff line number Diff line
@@ -36,7 +36,7 @@ def test_calculate_doublet_performance_precise(): 




    # Assert

    rtol = 0.002

    rtol = 0.002 # accurate to 0.2%

    assert np.isclose(result.flow, 227.2757568359375, rtol=rtol)

    assert np.isclose(result.pres, 60, rtol=rtol)

    assert np.isclose(result.utc, 6.616096470753937, rtol=rtol)
@@ -81,7 +81,7 @@ def test_calculate_doublet_performance_approximate(): 




    # Assert

    rtol = 0.005

    rtol = 0.005 # accurate to 0.5%

    assert np.isclose(result.flow, 227.2757568359375, rtol=rtol)

    assert np.isclose(result.pres, 60, rtol=rtol)

    assert np.isclose(result.utc, 6.616096470753937, rtol=rtol)