Add Major (Darcy) and Minor (Fittings) Losses as Pressure Drop Option for Borefields

[Vbraciszewski]

@mwetter, @LoneMeertens,

Connects with Issues: #4483 and #4655.

To derisk systems and optimize pump control, it is useful to have more advanced pressure drop calculation for piping and borefield models that are aware of dynamically changing friction factor (based on changes to reynolds number especially for antifreeze mixtures).

Priority is to add friction losses (Major Losses): Add Darcy-Weisbach based option that considers u-bend diameter, length, and friction factor based on reynold's number (includes changing viscosity).

For inclusion in borefield model, we could borrow from: Buildings.Fluid.FixedResistances.HydraulicDiameter but this appears to just calculate friction factor one time, then hold it constant for simulation.

Options we are considering:

1. Extend MBL hydraulic diameter & borefield to include dynamic reynolds number calc
2. Use MSL pipe (MSL StaticPipe with WallFriction.Detailed) with MBL borefield - this appears to calculate Re# dynamically throughout simulation.

Future Steps
Minor losses (from fittings e.g., u-bend fitting): ability to sum extra minor losses based on number of elbows/tees. Unclear if there is precedent in MBL for this.

[mwetter]

I think adding a new model such as DetailedFriction (or another better name?) would be good.
I would not change the HydraulicDiameter as in early design of duct and pipe networks, one often has only limited information about bends etc that in my view does not justify a more expensive model because the uncertainty of the lack of design information is likely higher than the approximation of constant friction factor used in the underlying model Buildings.Fluid.FixedResistances.PressureDrop.
Such a model can then take the viscosity from the medium model, and the medium model determines whether viscosity is constant or a function of temperature (as is usually done in Modelica fluid and media models).

Regarding fittings and bends, this has not been implemented. For now, we often simply multiply the linear pipe resistance by a factor of 2 to account for such losses (using a rule-of-thumb that I was thought in school...).
If there are resources, we could add models that can be used for the calculation of the $\zeta$ value (ratio of static pressure loss over dynamic pressure). One however has to keep in mind that the equations or charts for $\zeta$ values are based on the assumption of undisturbed flow at the inlet and outlet. Recknagel has a figure that shows that multiple bends in series have a significantly higher $\zeta$ value as the sum of the individual bends. This makes me question whether such models of individual bends are indeed much more useful than the above rule of thumb. I guess it depends on the application, data availability and accuracy requirements.

[LoneMeertens]

Progress update — dynamic Darcy-Weisbach pressure drop for vertical borefield pipes

This update adds a dynamic Darcy-Weisbach pressure-drop calculation for the vertical borefield pipe circuits in Buildings.Fluid.Geothermal.Borefields.

Implementation branch:

• Branch: issue4656_PressureDropBorefield
• Commit: 873533fbf254e224f98088fd958c308528ebe6f3

This is a first step toward a more detailed borefield hydraulic model. The horizontal/header piping is not included yet and will be handled later in a dedicated detailed pipe component.

---

Scope

Implemented so far:

• Dynamic major pressure losses in the vertical U-tube circuits.
• Pressure drop based on instantaneous mass flow rate.
• Churchill friction factor for the Darcy-Weisbach calculation.
• Backward-compatible switch via use_DarcyPressureDrop.
• Shared pipe roughness parameter for pressure drop and friction-factor-based heat transfer.
• Validation models for constant and variable mass flow.
• Total vertical GHE fluid volume output through VTubBorFie.

Not included yet:

• Horizontal/header piping.
• Header/manifold losses.
• Temperature-dependent fluid properties in the pressure-drop calculation.
• Minor losses from U-bends or fittings.

---

Main implementation

New pressure-loss function

Added:

Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.BaseClasses.Functions.pressureLossCircularPipe

This function computes circular-pipe major pressure loss using Darcy-Weisbach with a Churchill friction factor.

The Reynolds number is computed as:

$$ Re = \frac{D |\dot m|}{A \mu} $$

The pressure-drop calculation uses the modified friction coefficient:

$$ \lambda_2 = \lambda Re^2 $$

For very small Reynolds numbers, the laminar limit is used:

$$ \lambda_2 = 64 Re $$

For larger Reynolds numbers:

$$ \lambda_2 = f_\mathrm{Churchill}(Re,\epsilon/D) Re^2 $$

The pressure drop is then computed as:

$$ \Delta p = \mathrm{sign}(\dot m) \frac{L \mu^2}{2 \rho D^3} \lambda_2 $$

This avoids using a Reynolds-number-clamped raw friction factor directly in the pressure-drop equation, which would otherwise give artificial low-flow behavior.

Validation model:

Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.BaseClasses.Functions.Validation.PressureLossCircularPipe

Figure 1 — Pressure-loss function validation

---

New pressure-drop component

Added:

Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.BaseClasses.PressureDropCircularPipe

This component is inserted into the vertical borehole fluid path and evaluates pressure drop from the instantaneous mass flow rate.

It enforces:

$$ \Delta p = p_a - p_b $$

and:

$$ \Delta p = f(\dot m) $$

using pressureLossCircularPipe.

The component has a Boolean switch:

computePressureDrop

If computePressureDrop=false, it imposes zero pressure drop and behaves as a hydraulic pass-through. This allows the component to be present without changing the previous nominal-pressure-drop behavior.

---

New configuration parameters

Updated:

Buildings.Fluid.Geothermal.Borefields.Data.Configuration.Template

Added:

• use_DarcyPressureDrop
  Enables the new Darcy-Weisbach pressure-drop calculation.

• roughness
  Absolute roughness of the vertical GHE pipes. Used by the pressure-drop calculation and by the friction-factor-based heat-transfer correlation.

• VTubBorFie
  Final derived parameter for the total fluid volume in the vertical borefield pipes.

Default behavior remains backward compatible:

use_DarcyPressureDrop = false

---

Integration in borehole models

Single U-tube

Updated:

Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.OneUTube

Added one pressure-drop component:

preDro

Hydraulic path:

port_a → preDro → intHex[1].port_a1 → U-tube circuit → port_b

For one single U-tube circuit:

$$ L = 2 h_\mathrm{Bor} $$

When use_DarcyPressureDrop=true, the old dp_nominal contribution is disabled and preDro computes the dynamic Darcy-Weisbach pressure drop.

When use_DarcyPressureDrop=false, preDro.dp = 0 and the old nominal-pressure-drop behavior is retained.

---

Double U-tube

Updated:

Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.TwoUTube

Added two pressure-drop components:

• preDro1
• preDro2

Each represents one full U-tube circuit:

$$ L = 2 h_\mathrm{Bor} $$

Parallel double U-tube

Topology:

branch 1: port_a → preDro1 → U-tube circuit 1 → port_b
branch 2: port_a → preDro2 → U-tube circuit 2 → port_b

Each branch sees about half the borehole flow:

$$ \dot m_\mathrm{branch} \approx \frac{\dot m_\mathrm{borehole}}{2} $$

Expected:

preDro1.dp ≈ preDro2.dp
borFie2UTubPar.dp ≈ preDro1.dp ≈ preDro2.dp

Series double U-tube

Topology:

port_a → preDro1 → U-tube circuit 1 → preDro2 → U-tube circuit 2 → port_b

The full flow passes through both circuits:

$$ \dot m_1 = \dot m_2 = \dot m_\mathrm{borehole} $$

The pressure drops add:

$$ \Delta p_\mathrm{series} = \Delta p_1 + \Delta p_2 $$

Expected:

borFie2UTubSer.dp ≈ preDro1.dp + preDro2.dp

---

Validation models

Constant mass-flow validation

Added:

Buildings.Fluid.Geothermal.Borefields.Examples.BorefieldsDarcyPressureDrop

This enables use_DarcyPressureDrop=true for all three configurations while keeping mass flow constant.

Expected ordering:

$$ |\Delta p_\mathrm{parallel}| < |\Delta p_\mathrm{single}| < |\Delta p_\mathrm{series}| $$

This was confirmed.

Figure 2 — Constant-flow Darcy pressure-drop validation

---

Variable mass-flow validation

Added:

Buildings.Fluid.Geothermal.Borefields.Examples.BorefieldsDarcyPressureDropVariableFlow

This enables Darcy-Weisbach pressure drop and drives the three borefields with sinusoidal mass-flow inputs. The flow varies between 50% and 150% of nominal.

Expected behavior:

m_flow increases → Re increases → dp increases
m_flow decreases → Re decreases → dp decreases

The pressure-drop ordering should remain:

$$ |\Delta p_\mathrm{parallel}| < |\Delta p_\mathrm{single}| < |\Delta p_\mathrm{series}| $$

throughout the simulation.

Figure 3 — Variable mass-flow inputs

Figure 4 — Dynamic pressure-drop response

---

Current limitations

The current implementation is dynamic with respect to mass flow:

m_flow changes → Re changes → f changes → dp changes

However, the pressure-drop calculation still uses default medium properties:

• rho_default
• mu_default

So this is not yet included:

T changes → mu(T), rho(T) change → Re changes → f changes → dp changes

This should be addressed in a later step, together with the planned temperature-dependent Rconv update.

The scope is currently limited to the vertical borefield pipe circuits. Horizontal/header piping and manifold pressure losses are not included yet and should be handled in a separate detailed horizontal piping component.

---

Open question — U-bend minor losses

Only major vertical pipe wall-friction losses are included so far. Minor losses from U-bends are currently neglected.

Since the number of U-bends is known, they could be added as an optional correction:

$$ \Delta p_\mathrm{total}= \left( K_\mathrm{major} + K_\mathrm{minor} \right) \frac{\dot m \left| \dot m \right|}{2\rho A^2}, \qquad K_\mathrm{major}=f\frac{L}{D}. $$

For the current example geometry, with $h_\mathrm{Bor}=100,\mathrm{m}$ and $D=0.036,\mathrm{m}$, one U-tube circuit has $L/D\approx5556$. For $f=0.02$ to $0.035$, this gives $K_\mathrm{major}\approx111$ to $194$.

A typical U-bend coefficient is $K_\mathrm{minor}\approx0.5$ to $2$, so the contribution is only about 0.3% to 2% of the vertical major loss. Even with $K_\mathrm{minor}=4$, it remains only about 2% to 4%.

@Vbraciszewski Decision: Include U-bend minor losses now as an optional vertical-borehole correction, or neglect them in the vertical model and add them later in the detailed horizontal/header piping component where minor losses are more relevant?

[Vbraciszewski]

Hi Lone -
Thanks for sharing the updates!

Minor Losses
My vote is to include minor losses in the vertical GHX for future flexibility (e.g., may be more significant when using different styles of GHX such as coaxial, see Dr. Spitler paper below) and provides comparison to industry tools (all which let you model minor losses). There are two ways this is done, equivalent length (add fitting length to darcy then calculate) or k factor (pull K factor from Crane technical paper 410 or other tables). There is not consensus on one approach for the vertical GHX (equivalent length is mentioned in Kavanaugh's ASHRAE book, but Dr. Spitler uses k factor which is more modern approach). More broadly, hydraulic analysis calculations use k factor so this is preferred for rest of system.

Options:

1. Equiv length only - most basic/quickest to implement (could be acceptable for quick solution for deep single and double u-bends given low overall impact)
2. k factor only - could potentially back calculate from equiv length or use tables
3. Toggle to switch between the two based on available data/user preference

Dr. Spitler's paper: https://www.sciencedirect.com/science/article/pii/S0375650522001924?ref=pdf_download&fr=RR-2&rr=a1e6f6a78966e612
"Vertical ground heat exchanger pressure loss – Experimental comparisons
and calculation procedures"

New Pressure Loss Function
What is the transition Re used (mentioned as "very small Re number") between the hard-code laminar calculation and Churchill and is it smoothed?

If there's a hard switch between the small Re laminar limit and Churchill correlation at some threshold Re (rather than a smooth blend, e.g. Modelica.Fluid style regStep/spline transition), that can introduce a small derivative discontinuity right at the switch point, which could matter for solver convergence at partial-load conditions.

Small comments
Do we need separate validation model that shows new PD calculation works when mass flow calc is inverted as function of pressure boundaries?

For clarity and support of clean documentation in the future, please consider identifying each term used in equations. For example, with reynolds number:

[LoneMeertens]

@Vbraciszewski, thanks for the detailed review and suggestions.

Good point on the “very small Re” wording. I updated the implementation and documentation to clarify that this low-Re treatment is not the laminar-to-turbulent transition. The Churchill correlation is still used for the all-regime friction behavior. The small-Re treatment is only a numerical guard to avoid evaluating the raw Darcy friction factor at zero flow.

The implementation now separates the raw Churchill friction factor from the modified pressure-drop coefficient, lambda2 = f*Re^2. The raw friction-factor function Buildings.Fluid.FixedResistances.Functions.churchillFrictionFactor no longer contains an internal Reynolds-number clamp and is documented as valid for Re > 0. For pressure-drop calculations, I added Buildings.Fluid.FixedResistances.Functions.churchillFrictionFactorRe2, which returns lambda2 directly. For very small Reynolds numbers it uses the laminar limit lambda2 = 64*Re.

The current guard is Re_min = 1e-6, so lambda2 = 64*Re for Re <= Re_min, and lambda2 = f_Churchill(Re, epsilon/D)*Re^2 for Re > Re_min. This is currently a small hard switch, not a Modelica.Fluid-style spline transition. However, it is placed at Re = 1e-6 and is only a zero-flow numerical guard. Since the Churchill correlation approaches the same laminar limit as Re -> 0, both expressions approach lambda2 = 64*Re, so the value and slope mismatch should be negligible. If preferred, I can still replace this with a smoothed regularization.

This also preserves the correct low-flow behavior, dp proportional to m_flow, rather than an artificial quadratic behavior near zero flow, which is important when the pressure-drop relation is solved inversely in a hydraulic network.

Agreed that the inverse pressure-driven case should be validated. I added component-level validations for PressureDropCircularPipe in addition to the existing function validation. The existing function validation is Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.BaseClasses.Functions.Validation.PressureLossCircularPipe, which checks the pure function direction m_flow -> dp.

I added a mass-flow-driven component validation, Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.BaseClasses.Validation.PressureDropCircularPipeMassFlowDriven, which imposes mass flow and checks m_flow -> dp. I also added a pressure-driven component validation, Buildings.Fluid.Geothermal.Borefields.BaseClasses.Boreholes.BaseClasses.Validation.PressureDropCircularPipePressureDriven, which imposes pressure difference using pressure boundaries and solves for mass flow, i.e. dp -> m_flow.

Both component validations use sinusoidal inputs that cross zero, so they check positive flow, reverse flow, and the low-flow region. They also recompute the pressure drop using pressureLossCircularPipe and compare it with the component pressure drop. The expected result is preDro.dp ≈ dpFun and errDp ≈ 0. For the pressure-driven case, the expected behavior is dpSet > 0 -> preDro.m_flow > 0, dpSet = 0 -> preDro.m_flow ≈ 0, and dpSet < 0 -> preDro.m_flow < 0.

Thanks also for the suggestion on minor losses. I agree that including them as an optional vertical GHX correction is useful for flexibility and for comparison with industry tools. I will review the Spitler paper you linked and include minor losses next.

[mwetter]

@LoneMeertens : The current hard switch at Re_min = 1e-6 needs to be make once continuously differentiable. Otherwise solvers can get stuck there. (And often, solvers iterate near zero flow rate as this region is hard for the solver due to numerical noise and possible reverse flow. So this is important for robustness.)

[LoneMeertens]

@mwetter, thanks for the clarification.

Just to make sure I understand the concern correctly: the purpose of the Re_min branch is only to avoid evaluating the raw Churchill friction factor at Re = 0, where the 64/Re behavior would be singular. It is not intended as a physical transition between laminar and turbulent flow.

Since the Churchill correlation approaches the laminar limit in this region, the two expressions

lambda2 = 64*Re

and

lambda2 = churchillFrictionFactor(Re, eps_D)*Re^2

should be nearly identical near Re_min = 1e-6. So I would expect the value jump, and probably also the derivative mismatch, to be very small.

Do you still recommend smoothing this switch to guarantee C1 continuity, even though the branch is only a zero-flow numerical guard and not a physical regime transition?

If yes, what would be your preferred implementation style for Buildings?

Option 1: use Buildings.Utilities.Math.Functions.regStep, giving a once continuously differentiable transition:

lambda2Lam := 64*Re;

lambda2 := smooth(1,
  if Re <= Re_min then
    lambda2Lam
  else
    Buildings.Utilities.Math.Functions.regStep(
      x=Re - (Re_min + Re_delta/2),
      y1=Buildings.Fluid.FixedResistances.Functions.churchillFrictionFactor(
           Re=Re,
           eps_D=eps_D)*Re^2,
      y2=lambda2Lam,
      x_small=Re_delta/2));

This would smooth over

Re_min < Re < Re_min + Re_delta

while still avoiding evaluation of churchillFrictionFactor at Re = 0.

Option 2: use Buildings.Utilities.Math.Functions.smoothHeaviside for a twice continuously differentiable transition:

lambda2Lam := 64*Re;

if Re <= Re_min then
  lambda2 := lambda2Lam;
else
  lambda2Chu :=
    Buildings.Fluid.FixedResistances.Functions.churchillFrictionFactor(
      Re=Re,
      eps_D=eps_D)*Re^2;

  w :=
    Buildings.Utilities.Math.Functions.smoothHeaviside(
      x=Re - (Re_min + Re_delta/2),
      delta=Re_delta/2);

  lambda2 := (1 - w)*lambda2Lam + w*lambda2Chu;
end if;

Do you prefer the simpler regStep approach, or would you prefer smoothHeaviside here, or is there another implementation you prefer?

[mwetter]

@LoneMeertens : I suggest using the smoothHeaviside formulation. The computing effort is similar between the two, but smoothHeaviside is C2 whereas regStep is C1.

If you must avoid any evaluation of churchillFrictionFactor outside of its validity range by the solver, even if the result is not used in the final accepted step, you also need to use noEvent.

[Vbraciszewski]

Regarding fittings and bends, this has not been implemented. For now, we often simply multiply the linear pipe resistance by a factor of 2 to account for such losses (using a rule-of-thumb that I was thought in school...). If there are resources, we could add models that can be used for the calculation of the ζ value (ratio of static pressure loss over dynamic pressure). One however has to keep in mind that the equations or charts for ζ values are based on the assumption of undisturbed flow at the inlet and outlet. Recknagel has a figure that shows that multiple bends in series have a significantly higher ζ value as the sum of the individual bends. This makes me question whether such models of individual bends are indeed much more useful than the above rule of thumb. I guess it depends on the application, data availability and accuracy requirements.

@mwetter - I realized I did not reply to this paragraph. For clarity, K factor is the same as ζ value (zeta) value (ratio of static pressure loss over dynamic pressure). At earliest stage of design, we may use rule of thumbs like what you mentioned. However, once we have a preliminary layout of the air or hydronic system we use darcy + minor losses (k factor) and accept inaccuracy inherent in 1D fluid simulation and available resources. Perhaps, designs with very close clustered fittings may add extra safety factor to final derived system pressure drop. However, with its limitations, K factor method is preferred and most common in 1D hydraulic simulation tools like Pipeflo and Datacor (formerly AFT) Fathom. Another limitation is that k factors are typically calculated under turbulent flow only. For sizing especially, that is acceptable. Future work could be to study Re# dependent k factors or some other more detailed implementation.

Link: https://docs.aft.com/fathom14/K-Factor-Loss-Model.html