Fix VRF FluidTCtrl negative outdoor unit fan power

[yujiex]

Pull request overview

• Fixes VRF FluidTCtrl model OU Fan power has negative values #10647
• Fixes VRFFluidTCtrl HR model division by zero when hard sized #10800
• Also fixes an issue where CompSpdActual is casted from real to int in function VRFOU_CompCap: the CompSpdActual (an int) that's fed into VRFOU_CompCap gets its value from the output of function VRFOU_CompSpd (which is a real). As a result, when VRFOU_CompCap gets its input CompSpdActual, it will cast it from real to int. When the value is small (less than 1), it will be zeroed out. The PR fixes the problem by changing the type of CompSpdActual from int to real.

Description of the fix

When running VariableRefrigerantFlow_FluidTCtrl_HR_5Zone.idf with Miami airport epw, using detailed timestep, with runperiod from 1/1 to 1/7, there are some negative fan power values like this

Tracing the code, we can see these values are calculated in VRFHR_OU_HR_Mode and stored in N_fan_OU

The following are related code snippets

 N_fan_OU = N_fan_OU_evap + N_fan_OU_cond;
...
 N_fan_OU_cond = this->RatedOUFanPower * m_air_cond / m_air_rated;
... 
m_air_cond = this->VRFOU_FlowRate(
            state, HXOpMode::CondMode, Tdischarge, this->SC, Q_h_OU, state.dataEnvrn->OutDryBulbTemp, state.dataEnvrn->OutHumRat);
...
// Inside VRFOU_FlowRate
    if (OperationMode == HXOpMode::CondMode) {
        // IU Cooling: OperationMode 0

        BF = this->RateBFOUCond; // 0.219;
        deltaT = this->C3Tc * pow_2(SHSC) + this->C2Tc * SHSC + this->C1Tc;
        T_coil_surf = TeTc - deltaT;
        T_coil_out = T_coil_in + (T_coil_surf - T_coil_in) * (1 - BF);
        m_air = Q_coil / (T_coil_out - T_coil_in) / 1005.0;   // <------------- here is the problem

    } else if (OperationMode == HXOpMode::EvapMode) {
    ...

The m_air calculation is related to the operation mode, which seems problematic.
The operation mode is related to the following two rps terms (rps1_evap <= rps2_cond indicates it's in mode 5 of heating dominant mode)

// Calculate compressor speed satisfying IU loads: rps1_evap & rps2_cond
this->VRFOU_CompSpd(state, Q_c_TU_PL, HXOpMode::EvapMode, temp_Tsuction, Tdischarge, h_IU_evap_in, h_IU_PLc_out, rps1_evap);
this->VRFOU_CompSpd(state, Q_h_TU_PL, HXOpMode::CondMode, temp_Tsuction, Tdischarge, h_IU_evap_in, h_IU_PLc_out, rps2_cond);

When the load is very small, rps2_cond could become negative when Q_cond_req < CompEvaporatingPWRSpd(1), as in VRFOU_CompSpd function else branch (Q_type != HXOpMode::EvapMode)

Q_cond_req = Q_req;
...
Q_evap_req = Q_cond_req - CompEvaporatingPWRSpd(CounterCompSpdTemp); // <-- Q_evap_req is negative, note here CounterCompSpdTemp = 1
...
CompSpdActual = this->CompressorSpeed(1) * (Q_evap_req * C_cap_operation) / CompEvaporatingCAPSpd(1); <-- negative as well

with this calculation, even when heating load is larger than cooling load, it might still report as cooling dominant

This PR fixes the problem by adding the following CompSpdActual calculation when Q_cond_req < CompEvaporatingPWRSpd(1)

if (Q_cond_req - CompEvaporatingPWRSpd(CounterCompSpdTemp) < 0.0) { // use compressor power to meet condenser required load
    CompSpdActual = this->CompressorSpeed(1) * (Q_cond_req * C_cap_operation) / CompEvaporatingPWRSpd(1);
}

Using the 3rd timestamp in the 8th hour (7:35) as an example, before the fix, the cooling load (Q_c_TU_PL) is about 74 and the heating load (Q_h_TU_PL) is about 125, but the HRmode is 2 (cooling dominant); after the fix, the HRmode became mode 5 (heating dominant) submode 2.

before vs after the fix (updated)

The negative fan power is removed after the fix. The following plot shows how the negative fan power changes after the fix (blue is before the fix, orange is after the fix), using Miami and Atlanta weather inputs.

Attached is the simulation output comparison table.
compare_fan_power_dev_feature.xlsx

NOTE: ENHANCEMENTS MUST FOLLOW A SUBMISSION PROCESS INCLUDING A FEATURE PROPOSAL AND DESIGN DOCUMENT PRIOR TO SUBMITTING CODE

Pull Request Author

Add to this list or remove from it as applicable. This is a simple templated set of guidelines.

• Title of PR should be user-synopsis style (clearly understandable in a standalone changelog context)
• Label the PR with at least one of: Defect, Refactoring, NewFeature, Performance, and/or DoNoPublish
• Pull requests that impact EnergyPlus code must also include unit tests to cover enhancement or defect repair
• Author should provide a "walkthrough" of relevant code changes using a GitHub code review comment process
• If any diffs are expected, author must demonstrate they are justified using plots and descriptions
• If changes fix a defect, the fix should be demonstrated in plots and descriptions
• If any defect files are updated to a more recent version, upload new versions here or on DevSupport
• If IDD requires transition, transition source, rules, ExpandObjects, and IDFs must be updated, and add IDDChange label
• If structural output changes, add to output rules file and add OutputChange label
• If adding/removing any LaTeX docs or figures, update that document's CMakeLists file dependencies

Reviewer

This will not be exhaustively relevant to every PR.

• Perform a Code Review on GitHub
• If branch is behind develop, merge develop and build locally to check for side effects of the merge
• If defect, verify by running develop branch and reproducing defect, then running PR and reproducing fix
• If feature, test running new feature, try creative ways to break it
• CI status: all green or justified
• Check that performance is not impacted (CI Linux results include performance check)
• Run Unit Test(s) locally
• Check any new function arguments for performance impacts
• Verify IDF naming conventions and styles, memos and notes and defaults
• If new idf included, locally check the err file and other outputs

[mjwitte]

If I'm understanding Q_coil correctly, it should never be <0. While this fixes the issue (which occurs only a few timesteps in the year) it masks a problem higher up. The Q that is passed in to this function shouldn't be <0, ever.

[yujiex]

So Q_coil is negative when it's in cooling mode. It is positive when it's heating

[mjwitte]

So Q_coil is negative when it's in cooling mode. It is positive when it's heating

I'm not sure that's the case. Q_coil here is the heat rejection when the unit is in cooling mode. So it's heating the outdoor air, so (Tcoilout-Tcoilin) should be positive, so you would expect Q_coil to be positive.

In the next block where it's in heating mode (so it's cooling the outdoor air) delta is (h_coil_in - h_coil_out so, again deltaH is positive and so Q_coil should be positive.

I ran in the debugger and saw both positive and negative values of Q_coil in the OperationMode == HXOpMode::CondMode block. Tried to watch the values in the OperationMode == HXOpMode::EvapMode block, but it never goes there. Even with Atlanta weather. That doesn't seem correct. Is it possible that the value passed in for OperationMode is not correct here?

[yujiex]

Now I'm not certain whether Q_coil should be always positive or not.
Although, from the definition of this function, it seems they do need to be absolute value of Q_coil.

Real64 VRFCondenserEquipment::VRFOU_FlowRate(EnergyPlusData &state,
                                             HXOpMode const OperationMode, // Mode 0 for running as condenser, 1 for evaporator
                                             Real64 const TeTc,            // VRF Tc at cooling mode, or Te at heating mode [C]
                                             Real64 const SHSC,            // SC for OU condenser or SH for OU evaporator [C]
                                             Real64 const Q_coil,          // absolute value of OU coil heat release or heat extract [W] <----- here they require absolute value
                                             Real64 const T_coil_in,       // Temperature of air at OU coil inlet [C]
                                             Real64 const W_coil_in        // Humidity ratio of air at OU coil inlet [kg/kg]
) const
{

Maybe I should just add a Q_coil = abs(Q_coil) at the beginning of the call to this function? @mjwitte

[mjwitte]

Maybe I should just add a Q_coil = abs(Q_coil) at the beginning of the call to this function?

No, I think that will mask other problems. To me, the "absolute value" comment implies that the value passed in will always be positive. We need to understand if HXOpMode is being passed in correctly and better understand the variables that are being passed in for Q_coil.

[yujiex]

I see. I will look into it. Currently there are definitely both positive and negatives for the same mode

[yujiex]

I've done some investigations and had a few findings. Although I'm not sure if this is enough to figure out a fix. @mjwitte

about the operation mode

According to the documentation, there are 6 heat recovery modes

• Mode 1: Cooling load only. No heating load. Both OU heat exchangers operate as condensers.
• Mode 2: Simultaneous heating and cooling. The sum of the zone cooling loads and compressor heat is much larger than the sum of the zone heating loads. Both OU heat exchangers operate as condensers.
• Mode 3: Simultaneous heating and cooling. The sum of the zone cooling loads and compressor heat is slightly larger than the sum of the zone heating loads. One OU heat exchanger operates as a condenser while the other as an evaporator.
• Mode 4: Simultaneous heating and cooling. The sum of the zone cooling loads and compressor heat is slightly smaller than the sum of the zone heating loads. One OU heat exchanger operates as a condenser while the other as an evaporator.
• Mode 5: Simultaneous heating and cooling. The sum of the zone cooling loads and compressor heat is much smaller than the sum of the zone heating loads. Both OU heat exchangers operate as evaporators.
• Mode 6: Heating load only. No cooling load. Both OU heat exchangers operate as evaporators.

Here is how the heat recovery mode is determined

// Determine FlagToLower
// this->DiffOUTeTo is "Difference between Outdoor Unit Evaporating Temperature and Outdoor Air Temperature in Heat Recovery Mode" from user input
if (state.dataEnvrn->OutDryBulbTemp - this->DiffOUTeTo < Tsuction) {
    temp_Tsuction = state.dataEnvrn->OutDryBulbTemp - this->DiffOUTeTo;
    FlagToLower = true;
} else {
    temp_Tsuction = Tsuction;
    FlagToLower = false;
}

// Calculate compressor speed satisfying IU loads: rps1_evap & rps2_cond
this->VRFOU_CompSpd(state, Q_c_TU_PL, HXOpMode::EvapMode, temp_Tsuction, Tdischarge, h_IU_evap_in, h_IU_PLc_out, rps1_evap);
this->VRFOU_CompSpd(state, Q_h_TU_PL, HXOpMode::CondMode, temp_Tsuction, Tdischarge, h_IU_evap_in, h_IU_PLc_out, rps2_cond);

// Determine FlagMode5
if (rps1_evap <= rps2_cond) {
    FlagMode5 = true;
} else {
    FlagMode5 = false;
}

// Determine HR Mode
if (FlagMode5) {
    HRMode = 5;
    if (FlagToLower)
        HRMode_sub = 1;
    else
        HRMode_sub = 2;
} else {

    if (FlagToLower)
        HRMode = 3; // Mode 3&4 share the same logics below
    else
        HRMode = 2;
}

One issue I saw is that rps2_cond is negative sometimes. Here's some debugging print. I used the
VariableRefrigerantFlow_FluidTCtrl_HR_5Zone.idf and USA_FL_Miami.Intl.AP.722020_TMY3.epw combination. For the 01/02 07, 4th time step, rps2_cond is negative.

01/02 07-4 HRMode=2, flag5=false, rps1_evap=14.796895201613392,rps2_cond=-103.04599155624868: Q_c_TU_PL=279.7809733457191, Q_h_TU_PL=300.59323456715885, temp_Tsuction=14.000000000000002,Tdischarge=36

This is caused by the following in VRFCondenserEquipment::VRFOU_CompSpd, where Q_req is smaller than the CompEvaporatingPWRSpd(CounterCompSpdTemp) when CounterCompSpdTemp = 1.

Q_cond_req = Q_req;
...
Q_evap_req = Q_cond_req - CompEvaporatingPWRSpd(CounterCompSpdTemp);

I was considering, maybe Q_evap_req = Q_cond_req - CompEvaporatingPWRSpd(CounterCompSpdTemp); should be Q_evap_req = max(0.0, Q_cond_req - CompEvaporatingPWRSpd(CounterCompSpdTemp));

Even with this, there still remains a question of how this "much larger" or "much smaller" in mode 2 and 5 is reflected in the code. I feel the current code seems to be just testing for "larger" or "smaller".

about the Q_coil sign

the input argument Q_h_OU in m_air_cond = this->VRFOU_FlowRate( state, HXOpMode::CondMode, Tdischarge, this->SC, Q_h_OU, state.dataEnvrn->OutDryBulbTemp, state.dataEnvrn->OutHumRat); corresponds to the Q_coil referred to in the last comment.

Here is how it's computed. As Q_c_tot + Ncomp is smaller than Q_h_TU_PL, the Q_h_OU is negative.

// compressor: Ncomp & Q_c_tot
this->VRFOU_CompCap(state, CompSpdActual, Tsuction, Tdischarge, h_IU_evap_in, h_comp_in, Q_c_tot, Ncomp);

Real64 Q_h_tot;                 // Total condenser capacity [W]
Real64 Q_h_TU_PL;                // Heating load to be met at heating mode, including the piping loss (W)
Real64 &Q_h_OU,        // OU condenser load [W]
Q_h_tot = Q_c_tot + Ncomp; 
Q_h_OU = Q_h_tot - Q_h_TU_PL;

Following is the debug print

01/02 07-4 Q_h_OU=-14.889397577117904=285.70383699004094 - 300.59323456715885

For this specific instance, the Q_h_TU_PL is mostly on the piping loss term.
Q_h_TU_PL = TU_HeatingLoad + Pipe_Q_h;

01/02 07-4 Q_h_TU_PL=300.59323456715885=0.8205014721511905 + 299.77273309500765

Maybe being in HR mode 2 should have ensured Q_c_tot + Ncomp much larger than Q_h_TU_PL? (as Q_c_tot and Ncomp is the compressor performance at the compressor speed rps1_evap that meets the cooling load Q_c_TU_PL). But currently it's not.

about the high piping loss

There's also some strange behavior about the very large piping loss.

Piping loss is calculated here, in the iterations (label230 region), other input arguments stay the same, the h_IU_cond_in updates, leading to updates in piping loss (Pipe_Q_h).

 this->VRFOU_PipeLossH(
            state, m_ref_IU_cond, max(min(Pcond, RefPHigh), RefPLow), h_IU_cond_in, OutdoorDryBulb, Pipe_Q_h, Pipe_DeltP_h, h_comp_out);

The convergence for the label230 checks whether h_comp_out and h_comp_out_new is close and whether h_IU_cond_in is maxed out. So before it's converged h_IU_cond_in keeps increasing till it reaches the upper bound h_IU_cond_in_up.

((std::abs(h_comp_out - h_comp_out_new) > Tolerance * h_comp_out) && (h_IU_cond_in < h_IU_cond_in_up))

I've added some debug print after this line.

        //* Update h_comp_out in iteration (Label230)
        h_comp_out_new = Ncomp / (m_ref_IU_evap + m_ref_OU_evap) + h_comp_in;

The following is the output for the iteration. h_comp_out_new stays the same while h_comp_out gets updated in function VRFOU_PipeLossH(). However, the h_comp_out is not going in the right direction. It's much larger than h_comp_out_new at the beginning and keeps increasing each round till h_IU_cond_in is maxed out. Large h_IU_cond_in will lead to high piping loss. This potentially contributes to the fact that Q_h_tot - Q_h_TU_PL is negative (Q_h_TU_PL gets too large as the piping loss term is large)

-------------before iteration !!! ----------------------------------
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405876 + 0) + 435658.9325403223, h_comp_out=25374199.615089763
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405876 + 0) + 435658.9325403223, h_comp_out=29927864.42684107
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405876 + 0) + 435658.9325403223, h_comp_out=34929581.8246451
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405876 + 0) + 435658.9325403223, h_comp_out=40339405.37992493
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405878 + 0) + 435658.9325403223, h_comp_out=46137187.431805536
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405878 + 0) + 435658.9325403223, h_comp_out=52302549.03700723
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405878 + 0) + 435658.9325403223, h_comp_out=58808375.97752049
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405878 + 0) + 435658.9325403223, h_comp_out=65654625.1845655
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405878 + 0) + 435658.9325403223, h_comp_out=72837775.9464217
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405876 + 0) + 435658.9325403223, h_comp_out=80330393.31041604
01/02 07-4 h_comp_out_new=451462.4133081403=20.990626990893254 / (0.0013282280846405878 + 0) + 435658.9325403223, h_comp_out=88762784.65486081
-------------after iteration !!! ----------------------------------

[rraustad]

So during these iterations h_comp_out_new did not change but h_comp_out did change? That makes no sense. It feels like this iteration loop isn't updating something that SHOULD change in the h_comp_out_new equation otherwise the loop should terminate.

[yujiex]

Yeah, for the specific time step when the fan power became negative, this is what happens. Maybe that could be another criteria (checking whether h_comp_out_new changes) to add to the convergence check

[rraustad]

Could we change the review process on this model to require a graph that shows how this fix changes the result? This model is so hard to understand that reviewing code is not productive. Let's look at a result and then discuss code. Adding new report variables that show certain data does not go negative might be a good path forward. I would hope that an iteration loop that hits a rail would be easily identified in a graph.

[yujiex]

Could we change the review process on this model to require a graph that shows how this fix changes the result? This model is so hard to understand that reviewing code is not productive. Let's look at a result and then discuss code. Adding new report variables that show certain data does not go negative might be a good path forward. I would hope that an iteration loop that hits a rail would be easily identified in a graph.

Good point. I'll produce plots illustrating before vs after the fix how fan power changes

[yujiex]

@rraustad I added a plot comparing the before and after OU fan power difference.

[rraustad]

That plot looks good. And there are no more negative condenser fan power results? I assume no so this should be ready.

[yujiex]

That plot looks good. And there are no more negative condenser fan power results? I assume no so this should be ready.

Right, just one time step with negative fan power.
I also just added that additional check that confirms whether h_comp_out_new changed or not

[dareumnam]

@yujiex, I think if you can at least explain why there's one time step with negative fan power, this PR will be good to go. Do you have any further comments, @rraustad?

[yujiex]

@yujiex, I think if you can at least explain why there's one time step with negative fan power, this PR will be good to go. Do you have any further comments, @rraustad?

Hi @dareumnam and @rraustad , some more details are in this comment #10649 (comment).

The high-level understanding I have is:

• The negative value is directly related to the Q_coil (OU heat release or extract) in function VRFOU_FlowRate.
• In the examples I've tested, this negative Q_coil is related to some high piping loss.
• The high piping loss is caused by the convergence logic in the iterations around label230. The loop converges if h_comp_out and h_comp_out_new are close enough. Currently, the iterations continues even if the distance between the two quantities gets larger in each followup iteration

label230:
// stuff
if (converge) {
    goto label230:
}

• The variable update of h_IU_cond_in (used in VRFOU_PipeLossH to compute h_comp_out) is modified so that it can also decrease in the next iteration depending on the sign of h_comp_out - h_comp_out_new (previously it can only increase). Also an early termination is added if h_comp_out_new doesn't change any more.

[mjwitte]

So the model has other issues that need to be fixed but we should never let the oa fan power be negative.

[dareumnam]

@mjwitte Thanks for your insight. I think we’ll need more time to fix this. I’ll move the milestone to 25.1.

[yujiex]

Thanks @mjwitte and @dareumnam. Even if the current change in convergence criteria fixed the negative fan power, maybe more extensive testing is needed

[yujiex]

Hi @mjwitte and @rraustad I might have found the culprit for the negative fan power. The negative rps2_cond might have made some heating-dominant cases be misclassified as cooling-dominant. The negative value seems to be from an un-handled edge case where Q_cond_req - CompEvaporatingPWRSpd(CounterCompSpdTemp). I have included some details in the PR description. Please see if it make sense. Thanks

[mjwitte]

@yujiex I think I understand the intent here. This could be simpler.

}  else {
    if (Q_evap_req < 0.0) {  // use compressor power to meet condenser required load
          CompSpdActual = this->CompressorSpeed(1) * (Q_cond_req * C_cap_operation) /CompEvaporatingPWRSpd(1);
    } else {
          CompSpdActual = this->CompressorSpeed(1) * (Q_evap_req * C_cap_operation) / CompEvaporatingCAPSpd(1);
    }

[yujiex]

Indeed, this is more clear.

[rraustad]

What I don't understand is what looks like a COP term = Qc / PWR. I just don't see how that gets back to a speed 1 fraction. What am I missing?

CompSpdActual = this->CompressorSpeed(1) * (Q_cond_req * C_cap_operation) /CompEvaporatingPWRSpd(1);

[yujiex]

@rraustad The way I think about this is when Q_cond_req - CompEvaporatingPWRSpd(1) <= 0.0, then what CompEvaporatingPWRSpd(1) can supply is too much to meet Q_cond_req. The speed should be somewhere between 0 and the 1st speed level. So in this case Q_cond_req should be met by CompEvaporatingPWRSpd(1) * ratio. This ratio is calculated as Q_cond_req / CompEvaporatingPWRSpd(1)

[rraustad]

But CompEvaporatingPWRSpd is power, not capacity. So Q_cond_req / CompEvaporatingPWRSpd(1) looks strange to me. Shouldn't both equations have CompEvaporatingCAPSpd in the denominator (line 13729)?

Array1D<Real64> CompEvaporatingPWRSpd; // Array for the compressor power at certain speed [W]
Array1D<Real64> CompEvaporatingCAPSpd; // Array for the evaporating capacity at certain speed [W]

[yujiex]

I can also try to change the VRF evaporative capacity to make it reach higher speed levels.

[yujiex]

This is how they look when they can reach higher speed levels (I reduced the capacity to 10000, the autosized heating and cooing capacity are around 70000W)

[rraustad]

Can you make it span speeds? like from speed 1 to speed 3? I want to see what Q_evap_req * C_cap_operating / CompEvaporatingCAPSpd(x) looks like (or your derivation). I suspect there are discontinuities because C_cap_operatiion = 0.6. And make a scatter plot like before. Q_evap_req versus compressor speed.

[rraustad]

Oh I see, your plot is Q_evap_req versus comp speed. I was confused because I didn't see a 45 degree line like I was expecting. Why is Q_evap_req = 2000 at comp speeds 0 - 10000 ??

[yujiex]

Somehow there's some weird thing happening when I plot with excel. CompSpdActual is always below 6000 for both the feature and develop branch output, but it somehow plotted it out to over 10000. I have attached the data here
plot data.zip
I used R to plot, this is how they look

[mjwitte]

Unit tests pass, defect file has no negative fan power, but the overall result changes very slightly (annual simulation). This is good once the code is updated as suggested.

[yujiex]

Thanks, @mjwitte . I just updated the code according to suggestions here: #10649 (comment)

[Myoldmopar]

Well that's weird...cppcheck couldn't be installed on the action runner? It's always something... :) I'll look into it.

[Myoldmopar]

If you feel like this needs more commits, pulling develop in will get the code integrity check passing. Don't waste CI cycles purely for that though, only if you need more relevant commits here.

[yujiex]

If you feel like this needs more commits, pulling develop in will get the code integrity check passing. Don't waste CI cycles purely for that though, only if you need more relevant commits here.

Hi @Myoldmopar, I don't think I have more commits relevant to this bug fix right now. I guess I will not pull develop just to fix the integrity check.