Merge f3bfad9 into 6110f12

README.md

@@ -146,6 +146,72 @@ cable_impedance = calculate_impedance(ConcentricNeutralCarsonsEquations(Cable())
 For examples of how to use the model, see the [concentric cable
 tests](https://github.com/opusonesolutions/carsons/blob/master/tests/test_concentric_neutral_cable.py).
 
+### Multi-Conductor Cable
+
+`carsons` also supports modelling of phased duplex, triplex, quadruplex cables and triplex secondary.
+It only requires a few more parameters to describe cable's geometry.
+
+```python
+from carsons import (MultiConductorCarsonsEquations,
+                     calculate_impedance)
+
+class Cable:
+    resistance: {
+        'A': per-length resistance of conductor A in ohm/meters
+        ...
+    }
+    geometric_mean_radius: {
+        'A': geometric mean radius of conductor A in meters
+        ...
+    }
+    wire_positions: {
+        'A': (x, y) cross-sectional position of conductor A in meters
+        ...
+    }
+    radius: {
+        'A': radius of conductor A
+        ...
+    }
+    insulation_thickness: {
+        'A': insulation thickness of conductor A
+        ...
+    }
+    phases: {'A', ... }
+
+cable_impedance = calculate_impedance(MultiConductorCarsonsEquations(Cable()))
+```
+
+To model a triplex secondary cable, the inputs should be keyed on secondary conductors `S1` and `S2`. The impedance result
+is a 2 x 2 matrix.
+
+```python
+class Cable:
+    resistance: {
+        'S1': per-length resistance of conductor A in ohm/meters
+        ...
+    }
+    geometric_mean_radius: {
+        'S1': geometric mean radius of conductor A in meters
+        ...
+    }
+    wire_positions: {
+        'S1': (x, y) cross-sectional position of conductor A in meters
+        ...
+    }
+    radius: {
+        'S1': radius of conductor A
+        ...
+    }
+    insulation_thickness: {
+        'S1': insulation thickness of conductor A
+        ...
+    }
+    phases: {'S1', ... }
+```
+
+For examples of how to use the model, see the [multi-conductor cable
+tests](https://github.com/opusonesolutions/carsons/blob/master/tests/test_multi_conductor.py).
+
 Problem Description
 -------------------
 

carsons/__init__.py

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

carsons/carsons.py

@@ -18,11 +18,12 @@ def convert_geometric_model(geometric_model) -> ndarray:
 
 def calculate_impedance(model) -> ndarray:
     z_primitive = model.build_z_primitive()
-    z_abc = perform_kron_reduction(z_primitive)
+    z_abc = perform_kron_reduction(z_primitive, dimension=model.dimension)
+
     return z_abc
 
 
-def perform_kron_reduction(z_primitive: ndarray) -> ndarray:
+def perform_kron_reduction(z_primitive: ndarray, dimension=3) -> ndarray:
     """ Reduces the primitive impedance matrix to an equivalent impedance
         matrix.
 
@@ -56,8 +57,10 @@ def perform_kron_reduction(z_primitive: ndarray) -> ndarray:
                             [Zba, Zbb, Zbc]
                             [Zca, Zcb, Zcc]
     """
-    Ẑpp, Ẑpn = z_primitive[0:3, 0:3], z_primitive[0:3, 3:]
-    Ẑnp, Ẑnn = z_primitive[3:,  0:3], z_primitive[3:,  3:]
+    Ẑpp, Ẑpn = (z_primitive[0:dimension, 0:dimension],
+                z_primitive[0:dimension, dimension:])
+    Ẑnp, Ẑnn = (z_primitive[dimension:,  0:dimension],
+                z_primitive[dimension:,  dimension:])
     Z_abc = Ẑpp - Ẑpn @ inv(Ẑnn) @ Ẑnp
     return Z_abc
 
@@ -78,17 +81,11 @@ def __init__(self, model):
         self.ω = 2.0 * π * self.ƒ  # angular frequency radians / second
 
     def build_z_primitive(self) -> ndarray:
-        neutral_conductors = sorted([
-            ph for ph in self.phases
-            if ph.startswith("N")
-        ])
-        conductors = ["A", "B", "C"] + neutral_conductors
-
-        dimension = len(conductors)
+        dimension = len(self.conductors)
         z_primitive = zeros(shape=(dimension, dimension), dtype=complex)
 
-        for index_i, phase_i in enumerate(conductors):
-            for index_j, phase_j in enumerate(conductors):
+        for index_i, phase_i in enumerate(self.conductors):
+            for index_j, phase_j in enumerate(self.conductors):
                 if phase_i not in self.phases or phase_j not in self.phases:
                     continue
                 R = self.compute_R(phase_i, phase_j)
@@ -190,8 +187,46 @@ def get_h(self, i):
         _, yᵢ = self.phase_positions[i]
         return yᵢ
 
+    @property
+    def dimension(self):
+        return 2 if getattr(self, 'is_secondary', False) else 3
+
+    @property
+    def conductors(self):
+        neutral_conductors = sorted([
+            ph for ph in self.phases
+            if ph.startswith("N")
+        ])
+
+        return ["A", "B", "C"] + neutral_conductors
+
+
+class ModifiedCarsonsEquations(CarsonsEquations):
+    """
+    Modified Carson's Equation. Two approximations are made:
+    only the first term of P and the first two terms of Q are considered.
+    """
+    number_of_P_terms = 1
+
+    def compute_P(self, i, j, number_of_terms=1) -> float:
+        return super().compute_P(i, j, self.number_of_P_terms)
+
+    def compute_X(self, i, j) -> float:
+        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 * π)
+
 
-class ConcentricNeutralCarsonsEquations(CarsonsEquations):
+class ConcentricNeutralCarsonsEquations(ModifiedCarsonsEquations):
     def __init__(self, model, *args, **kwargs):
         super().__init__(model)
         self.neutral_strand_gmr: Dict[str, float] = model.neutral_strand_gmr
@@ -246,22 +281,45 @@ def compute_d(self, i, j) -> float:
             # Distance between two neutral/phase conductors
             return distance_ij
 
-    def compute_X(self, i, j) -> float:
-        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) -> float:
         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)
+
+
+class MultiConductorCarsonsEquations(ModifiedCarsonsEquations):
+    def __init__(self, model):
+        super().__init__(model)
+        self.radius: Dict[str, float] = model.radius
+        self.insulation_thickness: Dict[str, float] = \
+            model.insulation_thickness
+
+    def compute_d(self, i, j) -> float:
+        # Assumptions:
+        # 1. All conductors in the cable are touching each other and
+        #    therefore equidistant.
+        # 2. In case of quadruplex cables, the space between conductors
+        #    which are diagonally positioned is neglected.
+        return (self.radius[i] + self.radius[j]
+                + self.insulation_thickness[i] + self.insulation_thickness[j])
+
+    @property
+    def conductors(self):
+        neutral_conductors = sorted([
+            ph for ph in self.phases if ph.startswith("N")
+        ])
+        if self.is_secondary:
+            conductors = ["S1", "S2"] + neutral_conductors
+        else:
+            conductors = ["A", "B", "C"] + neutral_conductors
+
+        return conductors
+
+    @property
+    def is_secondary(self):
+        phase_conductors = [ph for ph in self.phases if not ph.startswith('N')]
+        if phase_conductors == ["S1", "S2"]:
+            return True
+        else:
+            return False

tests/helpers.py

@@ -88,3 +88,45 @@ def diameter_over_neutral(self):
     @property
     def neutral_strand_count(self):
         return self._neutral_strand_count
+
+
+class MultiLineModel:
+    def __init__(self, conductors, is_secondary=False):
+        self._resistance = {}
+        self._geometric_mean_radius = {}
+        self._wire_positions = {}
+        self._radius = {}
+        self._insulation_thickness = {}
+
+        for phase, val in conductors.items():
+            self._resistance[phase] = val['resistance']
+            self._geometric_mean_radius[phase] = val['gmr']
+            self._wire_positions[phase] = val['wire_positions']
+            self._radius[phase] = val['radius']
+            self._insulation_thickness[phase] = val['insulation_thickness']
+
+        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 radius(self):
+        return self._radius
+
+    @property
+    def insulation_thickness(self):
+        return self._insulation_thickness