Support concentric neutral cables

README.rst

@@ -64,19 +64,19 @@ of the conductor for that phase.
 .. code:: python
 
 
-   from carsons import CarsonsEquations, perform_kron_reduction
+    from carsons import CarsonsEquations, perform_kron_reduction, impedance
 
     class Line:
        gmr: {
            'A': geometric_mean_radius_A
            ...
        }
        r: {
-            'A' => per-length resistance of conductor A in ohms
+            'A': per-length resistance of conductor A in ohms
             ...
        }
        phase_positions: {
-            'A' => (x, y) cross-sectional position of the conductor in meters
+            'A': (x, y) cross-sectional position of the conductor in meters
             ...
        }
        phases: {'A', ... }
@@ -86,6 +86,7 @@ of the conductor for that phase.
 
     z_primitive = CarsonsEquations(Line()).build_z_primitive()
     z_abc = perform_kron_reduction(z_primitive)
+    line_impedance = impedance(CarsonsEquations(Line()))
 
 
 The model supports any combination of ABC phasings (for example BC, BCN etc...)
@@ -101,6 +102,63 @@ For examples of how to use the model, see the `tests <https://github.com/opusone
 ``carsons`` is tested against several cable configurations from the
 `IEEE 4-bus test network <http://sites.ieee.org/pes-testfeeders/resources/>`_.
 
+
+### Concentric Neutral Cable
+
+``carsons`` also supports modelling of concentric neutral cables of any phasings.
+Its usage is very similar to the example above, only requires a few more
+parameters about the neutral conductors in the line model object.
+
+.. code:: python
+
+
+    from carsons import (ConcentricNeutralCarsonsEquations,
+                         perform_kron_reduction,
+                         impedance)
+
+    class Line:
+       resistance: {
+           'A': per-length resistance of conductor A in ohms
+           ...
+       }
+       geometric_mean_radius: {
+           'A': geometric_mean_radius_A
+           ...
+       }
+       phase_positions: {
+            'A' => (x, y) cross-sectional position of the conductor in meters
+            ...
+       }
+       phases: {'A', 'NA', ... }
+       neutral_strand_gmr: {
+           'NA': neutral_strand_gmr_A
+           ...
+       }
+       neutral_strand_resistance: {
+           'NA': neutral_strand_resistance_A
+           ...
+       }
+       neutral_strand_diameter: {
+           'NA': neutral_strand_diameter_A
+           ...
+       }
+       diameter_over_neutral: {
+           'NA': diameter_over_neutral_A
+           ...
+       }
+       neutral_strand_count: {
+           'NA': neutral_strand_count_A
+           ...
+       }
+
+
+    z_primitive = ConcentricNeutralCarsonsEquations(Line()).build_z_primitive()
+    z_abc = perform_kron_reduction(z_primitive)
+    line_impedance = impedance(ConcentricNeutralCarsonsEquations(Line()))
+
+For examples of how to use the model, see the `tests <https://github.com/opusonesolutions/carsons/blob/master/tests/test_concentric_neutral_cable.py>`_.
+
+
 Problem Description
 -------------------
 

carsons/__init__.py

@@ -1,3 +1,5 @@
-from carsons.carsons import convert_geometric_model  # noqa 401
+from carsons.carsons import (convert_geometric_model,               # noqa 401
+                             impedance,                             # noqa 401
+                             ConcentricNeutralCarsonsEquations)     # noqa 401
 
 name = "carsons"

carsons/carsons.py

@@ -1,13 +1,9 @@
-from numpy import pi as π
+from collections import defaultdict
+from itertools import islice
 
-from numpy import zeros
+from numpy import arctan, cos, log, sin, sqrt, zeros
+from numpy import pi as π
 from numpy.linalg import inv
-from numpy import sqrt
-from numpy import log
-from numpy import cos
-from numpy import sin
-from numpy import arctan
-from itertools import islice
 
 
 def convert_geometric_model(geometric_model):
@@ -18,6 +14,12 @@ def convert_geometric_model(geometric_model):
     return z_abc
 
 
+def impedance(model):
+    z_primitive = model.build_z_primitive()
+    z_abc = perform_kron_reduction(z_primitive)
+    return z_abc
+
+
 def perform_kron_reduction(z_primitive):
     """ Reduces the primitive impedance matrix to an equivalent impedance
         matrix.
@@ -107,16 +109,17 @@ def compute_X(self, i, j):
         Qᵢⱼ = self.compute_Q(i, j)
         ΔX = self.μ * self.ω / π * Qᵢⱼ
 
+        # calculate geometry ratio 𝛥G
         if i != j:
             Dᵢⱼ = self.compute_D(i, j)
             dᵢⱼ = self.compute_d(i, j)
-            geometry_ratio = Dᵢⱼ / dᵢⱼ
+            𝛥G = Dᵢⱼ / dᵢⱼ
         else:
             hᵢ = self.get_h(i)
             gmrⱼ = self.gmr[j]
-            geometry_ratio = 2.0 * hᵢ / gmrⱼ
+            𝛥G = 2.0 * hᵢ / gmrⱼ
 
-        X_o = self.ω * self.μ / (2 * π) * log(geometry_ratio)
+        X_o = self.ω * self.μ / (2 * π) * log(𝛥G)
 
         return X_o + ΔX
 
@@ -172,6 +175,7 @@ def compute_d(self, i, j):
 
     def compute_D(self, i, j):
         xⱼ, yⱼ = self.phase_positions[j]
+
         return self.calculate_distance(self.phase_positions[i], (xⱼ, -yⱼ))
 
     @staticmethod
@@ -183,3 +187,74 @@ def calculate_distance(positionᵢ, positionⱼ):
     def get_h(self, i):
         _, yᵢ = self.phase_positions[i]
         return yᵢ
+
+
+class ConcentricNeutralCarsonsEquations(CarsonsEquations):
+    def __init__(self, model, *args, **kwargs):
+        super().__init__(model)
+        self.neutral_strand_gmr = model.neutral_strand_gmr
+        self.neutral_strand_count = defaultdict(
+            lambda: None, model.neutral_strand_count)
+        self.neutral_strand_resistance = model.neutral_strand_resistance
+        self.radius = defaultdict(lambda: None, {
+            phase: (diameter_over_neutral -
+                    model.neutral_strand_diameter[phase]) / 2
+            for phase, diameter_over_neutral
+            in model.diameter_over_neutral.items()
+        })
+        self.phase_positions.update({
+            f"N{phase}": self.phase_positions[phase]
+            for phase in self.phase_positions.keys()
+        })
+        self.gmr.update({
+            phase: self.GMR_cn(phase)
+            for phase in model.diameter_over_neutral.keys()
+        })
+        self.r.update({
+            phase: resistance / model.neutral_strand_count[phase]
+            for phase, resistance in model.neutral_strand_resistance.items()
+        })
+        return
+
+    def compute_d(self, i, j):
+        I, J = set(i), set(j)
+        r = self.radius[i] or self.radius[j]
+
+        one_neutral_same_phase = I ^ J == set('N')
+        different_phase = not I & J
+        one_neutral = 'N' in I ^ J
+
+        if one_neutral_same_phase:
+            # Distance between a neutral/phase conductor of same phase
+            return r
+
+        distance_ij = self.calculate_distance(self.phase_positions[i],
+                                              self.phase_positions[j])
+        if different_phase and one_neutral:
+            # Distance between a neutral/phase conductor of different phase
+            # approximate by modelling the concentric neutral cables as one
+            # equivalent conductor directly above the phase conductor
+            return (distance_ij**2 + r**2) ** 0.5
+        else:
+            # Distance between two neutral/phase conductors
+            return distance_ij
+
+    def compute_X(self, i, j):
+        Q_first_term = super().compute_Q(i, j, 1)
+
+        # Simplify equations and don't compute Dᵢⱼ explicitly
+        kᵢⱼ_Dᵢⱼ_ratio = sqrt(self.ω * self.μ / self.ρ)
+        ΔX = Q_first_term * 2 + log(2)
+
+        if i == j:
+            X_o = -log(self.gmr[i]) - log(kᵢⱼ_Dᵢⱼ_ratio)
+        else:
+            X_o = -log(self.compute_d(i, j)) - log(kᵢⱼ_Dᵢⱼ_ratio)
+
+        return (X_o + ΔX) * self.ω * self.μ / (2 * π)
+
+    def GMR_cn(self, phase):
+        GMR_s = self.neutral_strand_gmr[phase]
+        k = self.neutral_strand_count[phase]
+        R = self.radius[phase]
+        return (GMR_s * k * R**(k-1))**(1/k)

setup.py

@@ -38,6 +38,6 @@ def readme():
         'numpy>=1.13.1',
     ],
     extras_require={
-        "test": ["pytest>=3.6", "pytest-cov"],
+        "test": ["pytest>=3.6", "pytest-cov", "pint"],
     },
 )

tests/helpers.py

@@ -26,3 +26,65 @@ def wire_positions(self):
     @property
     def phases(self):
         return self._phases
+
+
+class ConcentricLineModel:
+    def __init__(self, conductors):
+        self._resistance = {}
+        self._geometric_mean_radius = {}
+        self._wire_positions = {}
+        self._phases = {}
+        self._neutral_strand_gmr = {}
+        self._neutral_strand_resistance = {}
+        self._neutral_strand_diameter = {}
+        self._diameter_over_neutral = {}
+        self._neutral_strand_count = {}
+
+        for phase, val in conductors.items():
+            if 'N' in phase:
+                self._neutral_strand_gmr[phase] = val['neutral_strand_gmr']
+                self._neutral_strand_resistance[phase] = val['neutral_strand_resistance']   # noqa 401
+                self._neutral_strand_diameter[phase] = val['neutral_strand_diameter']       # noqa 401
+                self._diameter_over_neutral[phase] = val['diameter_over_neutral']           # noqa 401
+                self._neutral_strand_count[phase] = val['neutral_strand_count']
+            else:
+                self._resistance[phase] = val['resistance']
+                self._geometric_mean_radius[phase] = val['gmr']
+                self._wire_positions[phase] = val['wire_positions']
+            self._phases = sorted(list(conductors.keys()))
+
+    @property
+    def resistance(self):
+        return self._resistance
+
+    @property
+    def geometric_mean_radius(self):
+        return self._geometric_mean_radius
+
+    @property
+    def wire_positions(self):
+        return self._wire_positions
+
+    @property
+    def phases(self):
+        return self._phases
+
+    @property
+    def neutral_strand_gmr(self):
+        return self._neutral_strand_gmr
+
+    @property
+    def neutral_strand_resistance(self):
+        return self._neutral_strand_resistance
+
+    @property
+    def neutral_strand_diameter(self):
+        return self._neutral_strand_diameter
+
+    @property
+    def diameter_over_neutral(self):
+        return self._diameter_over_neutral
+
+    @property
+    def neutral_strand_count(self):
+        return self._neutral_strand_count