[PL] BHE; Fix calculateTMFLowPropsAnnulus() impl.

This is used in the CXA/CXC implementations.

ProcessLib/HeatTransportBHE/BHE/Physics.h

@@ -59,6 +59,43 @@ inline double nusseltNumber(double const reynolds_number,
                       (std::pow(prandtl_number, 2.0 / 3.0) - 1.0)) *
            (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0));
 }
+
+// Pipe aspect ratio is the pipe's diameter equivalent over pipe's length.
+inline double nusseltNumberAnnulus(double const reynolds_number,
+                                   double const prandtl_number,
+                                   double const diameter_ratio,
+                                   double const pipe_aspect_ratio)
+{
+    if (reynolds_number < 2300.0)
+    {
+        return 3.66 + (4.0 - 0.102 / (diameter_ratio + 0.02)) *
+                          std::pow(diameter_ratio, 0.04);
+    }
+    if (reynolds_number < 10000.0)
+    {
+        double const gamma = (reynolds_number - 2300) / (10000 - 2300);
+
+        return (1.0 - gamma) *
+                   (3.66 + (4.0 - 0.102 / (diameter_ratio + 0.02))) *
+                   std::pow(diameter_ratio, 0.04) +
+               gamma *
+                   ((0.0308 / 8.0 * 1.0e4 * prandtl_number) /
+                    (1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
+                               (std::pow(prandtl_number, 2.0 / 3.0) - 1.0)) *
+                    (1.0 + std::pow(pipe_aspect_ratio, 2.0 / 3.0)) *
+                    ((0.86 * std::pow(diameter_ratio, 0.84) + 1.0 -
+                      0.14 * std::pow(diameter_ratio, 0.6)) /
+                     (1.0 + diameter_ratio)));
+    }
+    double const xi = std::pow(1.8 * std::log10(reynolds_number) - 1.5, -2.0);
+    return (xi / 8.0 * reynolds_number * prandtl_number) /
+           (1.0 + 12.7 * std::sqrt(xi / 8.0) *
+                      (std::pow(prandtl_number, 2.0 / 3.0) - 1.0)) *
+           (1.0 + std::pow(pipe_aspect_ratio, 2.0 / 3.0)) *
+           ((0.86 * std::pow(diameter_ratio, 0.84) + 1.0 -
+             0.14 * std::pow(diameter_ratio, 0.6)) /
+            (1.0 + diameter_ratio));
+}
 }  // end of namespace BHE
 }  // end of namespace HeatTransportBHE
 }  // end of namespace ProcessLib

ProcessLib/HeatTransportBHE/BHE/Pipe.h

@@ -29,11 +29,20 @@ struct Pipe
     double const wall_thickness;
     double const wall_thermal_conductivity;
 
+    /// Area of the pipe's inside without the wall.
     double area() const
     {
         constexpr double pi = boost::math::constants::pi<double>();
         return pi * diameter * diameter / 4;
     }
+
+    /// Area of the pipe's outside including the wall thickness.
+    double outsideArea() const
+    {
+        constexpr double pi = boost::math::constants::pi<double>();
+        double const d = diameter + 2 * wall_thickness;
+        return pi * d * d / 4;
+    }
 };
 
 Pipe createPipe(BaseLib::ConfigTree const& config);

ProcessLib/HeatTransportBHE/BHE/ThermoMechanicalFlowProperties.h

@@ -46,22 +46,28 @@ calculateThermoMechanicalFlowPropertiesPipe(Pipe const& pipe,
 
 inline ThermoMechanicalFlowProperties
 calculateThermoMechanicalFlowPropertiesAnnulus(
-    PipeParameters const& pipe, double const length,
+    Pipe const& inner_pipe, Pipe const& outer_pipe, double const length,
     RefrigerantProperties const& fluid, double const flow_rate)
 {
     double const Pr =
-        prandtlNumber(fluid.mu_r, fluid.specific_heat_capacity, fluid.lambda_r);
+        prandtlNumber(fluid.dynamic_viscosity, fluid.specific_heat_capacity,
+                      fluid.thermal_conductivity);
+
+    double const inner_pipe_outside_diameter =
+        inner_pipe.diameter + 2. * inner_pipe.wall_thickness;
 
+    // Velocity between the outer pipe and inner pipe.
     double const velocity =
-        annulusFlowVelocity(flow_rate, pipe.r_outer, pipe.r_inner + pipe.b_in);
+        flow_rate / (outer_pipe.area() - inner_pipe.outsideArea());
 
     double const Re = reynoldsNumber(
-        velocity, 2.0 * (pipe.r_outer - (pipe.r_inner + pipe.b_in)), fluid.mu_r,
-        fluid.density);
+        velocity, outer_pipe.diameter - inner_pipe_outside_diameter,
+        fluid.dynamic_viscosity, fluid.density);
 
-    double const diameter_ratio = (pipe.r_inner + pipe.b_in) / pipe.r_outer;
+    double const diameter_ratio =
+        inner_pipe_outside_diameter / outer_pipe.diameter;
     double const pipe_aspect_ratio =
-        (2 * pipe.r_outer - 2 * (pipe.r_inner + pipe.b_in)) / length;
+        (outer_pipe.diameter - inner_pipe_outside_diameter) / length;
     double const nusselt_number =
         nusseltNumberAnnulus(Re, Pr, diameter_ratio, pipe_aspect_ratio);
     return {velocity, nusselt_number};