1ph short circuit: calculate ip/ith, and honour the neutral earthing impedance of trafo3w

CHANGELOG.rst

@@ -5,6 +5,8 @@ Change Log
 -------------------------------
 - [FIXED] restored ``OpenDSSDirect.py`` to the ``all``/``dev`` extras so the OpenDSS converter is exercised (and its coverage reported) in CI again; it was dropped in #3062 because installing it alongside ``pytest~=9.1`` crashed the pytest process on Windows. Root cause (see `dss-extensions/OpenDSSDirect.py#148 <https://github.com/dss-extensions/OpenDSSDirect.py/issues/148>`_): pytest enables Python's ``faulthandler`` by default, which intercepts a first-chance Windows structured exception that OpenDSSDirect.py's native backend raises -- and normally handles itself -- during import, and misreports it as fatal. Bracketing the import with ``faulthandler.disable()``/``.enable()`` avoids the false crash while leaving ``faulthandler`` protecting the rest of the test run; the converter now runs on Windows instead of merely skipping there.
 - [ADDED] OpenDSS converter: series (bus-to-bus) ``Reactor`` elements are now imported as a fixed-impedance ``line``, the pattern some feeder libraries (e.g. EPRI's Ckt5/Ckt7) use to model the substation's Thevenin-equivalent source impedance instead of a ``Transformer``.
+- [FIXED] short circuit: ``calc_sc(..., fault="1ph", ip=True, ith=True)`` accepted ``ip``/``ith`` but never calculated them, so ``res_bus_sc.ip_ka`` and ``ith_ka`` came back all ``NaN`` without any warning
+- [ADDED] short circuit: neutral earthing impedance of a three winding transformer in the zero sequence, per winding via the optional columns ``xn_hv_ohm``/``rn_hv_ohm``, ``xn_mv_ohm``/``rn_mv_ohm``, ``xn_lv_ohm``/``rn_lv_ohm`` (the ``trafo3w`` counterpart of ``net.trafo.xn_ohm``/``rn_ohm``)
 
 [3.5.4] - 2026-07-08
 -------------------------------

pandapower/pd2ppc_zero.py

@@ -612,6 +612,68 @@ def _add_impedance_sc_impedance_zero(net, ppc):
     branch[f:t, BR_STATUS] = net["impedance"]["in_service"].values.astype(np.int64)
 
 
+def _split_trafo3w_vector_group(vector_group):
+    """
+    Split the vector group of a three winding transformer into the connection
+    code of its hv, mv and lv winding, e.g. "YNyd" -> ["yn", "y", "d"].
+    """
+    vg = vector_group.lower()
+    windings, ix = [], 0
+    while ix < len(vg):
+        if vg[ix + 1:ix + 2] == "n":
+            windings.append(vg[ix:ix + 2])
+            ix += 2
+        else:
+            windings.append(vg[ix])
+            ix += 1
+    return windings
+
+
+def _get_trafo3w_zn(t3, column):
+    value = t3.get(column, 0.)
+    if value is None:
+        return 0.
+    value = float(value)
+    return 0. if math.isnan(value) else value
+
+
+def _add_trafo3w_neutral_earthing_impedance(net, ppc, r, x, n_t3):
+    """
+    Add the neutral earthing impedance of a three winding transformer to the
+    zero sequence branch of the winding it belongs to. A star point earthed
+    through Z_N appears as 3*Z_N in series with the zero sequence impedance of
+    that winding (IEC 60909-0), the same way net.trafo.xn_ohm / rn_ohm is
+    handled for two winding transformers. The columns are given per winding:
+    xn_hv_ohm / rn_hv_ohm, xn_mv_ohm / rn_mv_ohm, xn_lv_ohm / rn_lv_ohm.
+
+    The branches of the star model are hv_bus -> star point, star point ->
+    mv_bus and star point -> lv_bus. Their impedance is referred to the to
+    side, and the star point is on the hv bus base voltage, so in all three
+    cases that is the base voltage of the winding's own terminal bus.
+    """
+    sides = ("hv", "mv", "lv")
+    columns = [f"{p}n_{side}_ohm" for side in sides for p in ("r", "x")]
+    if not any(col in net.trafo3w.columns for col in columns):
+        return
+    bus_lookup = net["_pd2ppc_lookups"]["bus"]
+    for t3_ix in range(n_t3):
+        t3 = net.trafo3w.iloc[t3_ix, :]
+        windings = _split_trafo3w_vector_group(t3.vector_group)
+        for side_ix, (side, winding) in enumerate(zip(sides, windings)):
+            # only an earthed star point carries a neutral earthing impedance
+            if not winding.endswith("n"):
+                continue
+            zn_ohm = complex(_get_trafo3w_zn(t3, f"rn_{side}_ohm"),
+                             _get_trafo3w_zn(t3, f"xn_{side}_ohm"))
+            if zn_ohm == 0:
+                continue
+            vn_kv = ppc["bus"][bus_lookup[int(t3[f"{side}_bus"])], BASE_KV]
+            zn_pu = 3 * zn_ohm / (vn_kv ** 2 / ppc["baseMVA"])
+            br_ix = t3_ix + side_ix * n_t3
+            r[br_ix] += zn_pu.real
+            x[br_ix] += zn_pu.imag
+
+
 def _add_trafo3w_sc_impedance_zero(net, ppc):
     # TODO Roman: check this/expand this
     branch_lookup = net["_pd2ppc_lookups"]["branch"]
@@ -766,6 +828,8 @@ def _add_trafo3w_sc_impedance_zero(net, ppc):
         else:
             raise UserWarning(f"{t3.vector_group} not supported yet for trafo3w!")
 
+    _add_trafo3w_neutral_earthing_impedance(net, ppc, r, x, n_t3)
+
     branch[f:t, BR_R] = r
     branch[f:t, BR_X] = x
     branch[f:t, BR_B] = 0

pandapower/shortcircuit/calc_sc.py

@@ -251,6 +251,13 @@ def _calc_sc_1ph(net, bus):
     _calc_rx(net, ppci_0, ppci_bus)
     _calc_ikss_1ph(net, ppci, ppci_0, ppci_bus)
 
+    # kappa is taken from the positive sequence network, as for the symmetrical
+    # fault: ip = kappa * sqrt(2) * ikss (IEC 60909-0, 4.6.1)
+    if net._options["ip"]:
+        _calc_ip(net, ppci)
+    if net._options["ith"]:
+        _calc_ith(net, ppci)
+
     if net._options["branch_results"]:
         if net._options["fault"] == "3ph":
             _calc_branch_currents_complex(net, ppci, ppci_bus)

pandapower/test/shortcircuit/test_1ph.py

@@ -592,5 +592,92 @@
     assert build("ZNyn", rn_ohm=5.0)[0] < znyn_mv
 
 
+@pytest.mark.parametrize("inverse_y", (True, False))
+@pytest.mark.parametrize("side, vector_group, faulted_bus", [
+    ("hv", "YNyd", 1), ("hv", "YNynd", 1),
+    ("mv", "YNynd", 2), ("mv", "YYnd", 2),
+    ("lv", "Yndyn", 3), ("lv", "Ydyn", 3)])
+def test_trafo3w_neutral_earthing_impedance(side, vector_group, faulted_bus, inverse_y):
+    # Same as test_trafo_neutral_earthing_impedance, for the three winding
+    # transformer: a star point earthed through Z_N adds 3*Z_N to the zero
+    # sequence impedance of that winding. The columns are per winding, since
+    # a trafo3w can have more than one earthed star point.
+    # The hv bus is fed through a delta so that the only zero-sequence earth of
+    # the transformer's terminals is the trafo3w itself -- otherwise the ext_grid
+    # would sit in parallel and the added Z_N would not show up undiluted.
+    # The rated voltages differ from the bus voltages on purpose: Z_N is an
+    # impedance at the terminal bus, so it must be referred with the bus base
+    # voltage and not with the rated voltage of the winding.
+    def build(rn_ohm=0., xn_ohm=0.):
+        net = create_empty_network(sn_mva=1.)
+        b_src = create_bus(net, vn_kv=220.)
+        b_hv = create_bus(net, vn_kv=110.)
+        b_mv = create_bus(net, vn_kv=20.)
+        b_lv = create_bus(net, vn_kv=10.)
+        create_ext_grid(net, b_src, s_sc_max_mva=5000., s_sc_min_mva=5000., rx_max=0.1, rx_min=0.1,
+                        x0x_max=1., x0x_min=1., r0x0_max=0.1, r0x0_min=0.1)
+        create_transformer_from_parameters(
+            net, b_src, b_hv, sn_mva=100., vn_hv_kv=220., vn_lv_kv=110., vk_percent=12., vkr_percent=0.4,
+            pfe_kw=0., i0_percent=0., vector_group="YNd", shift_degree=330., vk0_percent=12.,
+            vkr0_percent=0.4, mag0_percent=100., mag0_rx=0., si0_hv_partial=0.9)
+        create_transformer3w_from_parameters(
+            net, hv_bus=b_hv, mv_bus=b_mv, lv_bus=b_lv, vn_hv_kv=115., vn_mv_kv=21., vn_lv_kv=10.5,
+            sn_hv_mva=40., sn_mv_mva=20., sn_lv_mva=20., vk_hv_percent=10., vk_mv_percent=11.,
+            vk_lv_percent=12., vkr_hv_percent=0.3, vkr_mv_percent=0.31, vkr_lv_percent=0.32, pfe_kw=0.,
+            i0_percent=0., vk0_hv_percent=10., vk0_mv_percent=11., vk0_lv_percent=12., vkr0_hv_percent=0.3,
+            vkr0_mv_percent=0.31, vkr0_lv_percent=0.32, vector_group=vector_group,
+            **{f"rn_{side}_ohm": rn_ohm, f"xn_{side}_ohm": xn_ohm})
+        calc_sc(net, fault="1ph", case="max", inverse_y=inverse_y)
+        return net.res_bus_sc.loc[faulted_bus]
+
+    base = build()
+    with_rn = build(rn_ohm=5.)
+    with_xn = build(xn_ohm=7.)
+
+    assert np.isfinite(base.xk0_ohm)
+    assert np.isclose(with_rn.rk0_ohm - base.rk0_ohm, 3 * 5., rtol=0, atol=1e-6)
+    assert np.isclose(with_rn.xk0_ohm, base.xk0_ohm, rtol=0, atol=1e-6)
+    assert np.isclose(with_xn.xk0_ohm - base.xk0_ohm, 3 * 7., rtol=0, atol=1e-6)
+    assert np.isclose(with_xn.rk0_ohm, base.rk0_ohm, rtol=0, atol=1e-6)
+    # the positive sequence impedance is unaffected by neutral earthing
+    assert np.isclose(with_rn.rk_ohm, base.rk_ohm, rtol=0, atol=1e-6)
+    assert np.isclose(with_xn.xk_ohm, base.xk_ohm, rtol=0, atol=1e-6)
+    # and the earth-fault current is correspondingly reduced
+    assert with_rn.ikss_ka < base.ikss_ka
+    assert with_xn.ikss_ka < base.ikss_ka
+
+
+@pytest.mark.parametrize("inverse_y", (True, False))
+def test_1ph_ip_ith(inverse_y):
+    # ip and ith were accepted for fault="1ph" but never calculated, so
+    # res_bus_sc.ip_ka / ith_ka came back all NaN without any warning.
+    net = create_empty_network(sn_mva=1.)
+    b_hv = create_bus(net, vn_kv=110.)
+    b_lv = create_bus(net, vn_kv=20.)
+    create_ext_grid(net, b_hv, s_sc_max_mva=2000., s_sc_min_mva=2000., rx_max=0.1, rx_min=0.1,
+                    x0x_max=1., x0x_min=1., r0x0_max=0.1, r0x0_min=0.1)
+    create_transformer_from_parameters(
+        net, b_hv, b_lv, sn_mva=40., vn_hv_kv=110., vn_lv_kv=20., vk_percent=12., vkr_percent=0.4,
+        pfe_kw=0., i0_percent=0., vector_group="Dyn", shift_degree=330., vk0_percent=12.,
+        vkr0_percent=0.4, mag0_percent=100., mag0_rx=0., si0_hv_partial=0.9)
+    create_line_from_parameters(net, b_lv, create_bus(net, vn_kv=20.), length_km=5., r_ohm_per_km=0.16,
+                                x_ohm_per_km=0.12, c_nf_per_km=0., max_i_ka=0.5, r0_ohm_per_km=0.5,
+                                x0_ohm_per_km=0.5, c0_nf_per_km=0., endtemp_degree=80.)
+
+    calc_sc(net, fault="3ph", case="max", ip=True, ith=True, tk_s=0.1, kappa_method="C", inverse_y=inverse_y)
+    kappa = net.res_bus_sc.ip_ka.values / (np.sqrt(2) * net.res_bus_sc.ikss_ka.values)
+
+    calc_sc(net, fault="1ph", case="max", ip=True, ith=True, tk_s=0.1, kappa_method="C", inverse_y=inverse_y)
+    ikss, ip, ith = (net.res_bus_sc[c].values for c in ("ikss_ka", "ip_ka", "ith_ka"))
+    assert np.all(np.isfinite(ip)) and np.all(np.isfinite(ith))
+    # kappa only depends on the positive sequence R/X, so it is the same as for
+    # the symmetrical fault at the same bus (IEC 60909-0, 4.6.1)
+    assert np.allclose(ip, np.sqrt(2) * kappa * ikss, rtol=0, atol=1e-9)
+    assert np.all(ith >= ikss)
+    # ip and ith stay unset when they are not requested
+    calc_sc(net, fault="1ph", case="max", inverse_y=inverse_y)
+    assert "ip_ka" not in net.res_bus_sc.columns or np.all(np.isnan(net.res_bus_sc.ip_ka.values))
+
+
 if __name__ == "__main__":
     pytest.main([__file__])

pandapower/test/shortcircuit/test_iec60909_4.py

@@ -88,8 +88,9 @@ def iec_60909_4():
                                          vkr_mv_percent=.5, vk_lv_percent=12, vkr_lv_percent=.5, vk0_hv_percent=12,
                                          vkr0_hv_percent=0.5, vk0_mv_percent=12, vkr0_mv_percent=0.5, vk0_lv_percent=12,
                                          vkr0_lv_percent=0.5, vector_group="Yynd", tap_changer_type="Ratio", tap_max=10,
-                                         tap_min=-10, tap_pos=0, tap_neutral=0, tap_side="hv", tap_step_percent=0.1)
-                                         # reactor is 100 Ohm
+                                         tap_min=-10, tap_pos=0, tap_neutral=0, tap_side="hv", tap_step_percent=0.1,
+                                         # T6 is the only earth of the 10 kV network, through a 100 Ohm reactor
+                                         xn_mv_ohm=100.)
 
     create_motor(net, b7, pn_mech_mw=5.0, cos_phi=0.88, cos_phi_n=0.88, efficiency_n_percent=97.5, vn_kv=10, rx=0.1,
                  lrc_pu=5)
@@ -386,17 +387,23 @@ def test_iec_60909_4_2ph():
     assert np.allclose(net.res_bus_sc.skss_mw.values[:10], np.array(skss), atol=1e-1)
 
 
-@pytest.mark.skip("1ph gen-close sc calculation still under develop")
 def test_iec_60909_4_1ph():
     net = iec_60909_4()
-    calc_sc(net, fault="1ph", case="max", ip=True, tk_s=0.1, kappa_method="C")
+    calc_sc(net, fault="1ph", case="max", ip=True, ith=True, tk_s=0.1, kappa_method="C")
 
     ikss = [24.6526, 15.9722, 10.4106, 9.0498, 17.0452, 0.06337, 0.0633, 0, 0.0001, 0.0001]
     ip = [60.9982, 40.5086, 24.2424, 20.5464, 42.8337, 0.1656, 0.1279, 0.0, 0.00025, 0.00033]
     # No ib for 1ph sc calculation
 
-    assert np.allclose(net.res_bus_sc.ikss_ka.values[:10], np.array(ikss),
-                       atol=1e-4)  # assert np.allclose(net.res_bus_sc.ip.values[:8], np.array(ip), rtol=1e-4)
+    assert np.allclose(net.res_bus_sc.ikss_ka.values[:10], np.array(ikss), atol=1e-4)
+    # ip = kappa * sqrt(2) * ikss, kappa from the positive sequence network; the
+    # reference values are given to four decimals, hence the same atol as for 2ph
+    assert np.allclose(net.res_bus_sc.ip_ka.values[:10], np.array(ip), atol=1e-3)
+    # ith is not part of the reference set, but must at least be finite and,
+    # with m >= 0 and n = 1, never below ikss
+    ith = net.res_bus_sc.ith_ka.values[:10]
+    assert np.all(np.isfinite(ith))
+    assert np.all(ith >= net.res_bus_sc.ikss_ka.values[:10] - 1e-12)
 
 
 def test_detect_power_station_units():