diff --git a/process/models/engineering/materials.py b/process/models/engineering/materials.py
index d42f673b93..fef43d7e1d 100644
--- a/process/models/engineering/materials.py
+++ b/process/models/engineering/materials.py
@@ -2,6 +2,7 @@
 
 import logging
 
+import numba
 import numpy as np
 
 logger = logging.getLogger(__name__)
@@ -48,6 +49,7 @@ def eurofer97_thermal_conductivity(temp: float, fw_th_conductivity: float) -> fl
     )
 
 
+@numba.njit(cache=True)
 def calculate_tresca_stress(
     stress_x: float | np.ndarray,
     stress_y: float | np.ndarray,
@@ -71,13 +73,16 @@ def calculate_tresca_stress(
         Tresca stress (maximum shear stress criterion) in Pa, defined as the
         maximum of |stress_x - stress_y|, |stress_y - stress_z|, |stress_x - stress_z|.
     """
-    return max(
-        abs(stress_x - stress_y),
-        abs(stress_y - stress_z),
-        abs(stress_x - stress_z),
+    return np.maximum(
+        np.maximum(
+            np.abs(stress_x - stress_y),
+            np.abs(stress_y - stress_z),
+        ),
+        np.abs(stress_x - stress_z),
     )
 
 
+@numba.njit(cache=True)
 def calculate_von_mises_stress(
     stress_x: float | np.ndarray,
     stress_y: float | np.ndarray,
diff --git a/process/models/tfcoil/base.py b/process/models/tfcoil/base.py
index cba4e39506..51458e7c54 100644
--- a/process/models/tfcoil/base.py
+++ b/process/models/tfcoil/base.py
@@ -19,6 +19,10 @@
 from process.data_structure.build_variables import TFCSRadialConfiguration
 from process.data_structure.pfcoil_variables import PFConductorModel
 from process.data_structure.physics_variables import DivertorNumberModels
+from process.models.engineering.materials import (
+    calculate_tresca_stress,
+    calculate_von_mises_stress,
+)
 
 logger = logging.getLogger(__name__)
 
@@ -3135,19 +3139,20 @@ def stresscl(
         # Tresca / Von Mises yield criteria calculations
         # -----------------------------
         # Array equation
-        s_shear_tf = np.maximum(
-            np.absolute(sig_tf_r - sig_tf_z), np.absolute(sig_tf_z - sig_tf_t)
+
+        s_shear_tf = calculate_tresca_stress(
+            stress_x=sig_tf_r, stress_y=sig_tf_t, stress_z=sig_tf_z
         )
 
         # Array equation
 
-        sig_tf_vmises = np.sqrt(
-            0.5e0
-            * (
-                (sig_tf_r - sig_tf_t) ** 2
-                + (sig_tf_r - sig_tf_z) ** 2
-                + (sig_tf_z - sig_tf_t) ** 2
-            )
+        sig_tf_vmises = calculate_von_mises_stress(
+            stress_x=sig_tf_r,
+            stress_y=sig_tf_t,
+            stress_z=sig_tf_z,
+            stress_shear_xy=0.0,
+            stress_shear_yz=0.0,
+            stress_shear_zx=0.0,
         )
 
         # Array equation
@@ -3163,9 +3168,23 @@ def stresscl(
             ):
                 # Addaped Von-mises stress calculation to WP strucure [Pa]
 
-                svmxz = sigvm(0.0e0, sig_tf_t[ii], sig_tf_z[ii], 0.0e0, 0.0e0, 0.0e0)
+                svmxz = calculate_von_mises_stress(
+                    stress_x=0.0e0,
+                    stress_y=sig_tf_t[ii],
+                    stress_z=sig_tf_z[ii],
+                    stress_shear_xy=0.0e0,
+                    stress_shear_yz=0.0e0,
+                    stress_shear_zx=0.0e0,
+                )
 
-                svmyz = sigvm(sig_tf_r[ii], 0.0e0, sig_tf_z[ii], 0.0e0, 0.0e0, 0.0e0)
+                svmyz = calculate_von_mises_stress(
+                    stress_x=sig_tf_r[ii],
+                    stress_y=0.0e0,
+                    stress_z=sig_tf_z[ii],
+                    stress_shear_xy=0.0e0,
+                    stress_shear_yz=0.0e0,
+                    stress_shear_zx=0.0e0,
+                )
                 sig_tf_vmises[ii] = max(svmxz, svmyz)
 
                 # Maximum shear stress for the Tresca yield criterion using CEA
@@ -3702,44 +3721,6 @@ def eyoung_parallel(
     return eyoung_j_3, a_3, poisson_j_perp_3
 
 
-@numba.njit(cache=True)
-def sigvm(sx: float, sy: float, sz: float, txy: float, txz: float, tyz: float) -> float:
-    """Calculates Von Mises stress in a TF coil
-
-    This routine calculates the Von Mises combination of
-    stresses (Pa) in a TF coil.
-
-    Parameters
-    ----------
-    sx :
-        In-plane stress in X direction [Pa]
-    sy :
-        In-plane stress in Y direction [Pa]
-    sz :
-        In-plane stress in Z direction [Pa]
-    txy :
-        Out-of-plane stress in X-Y plane [Pa]
-    txz :
-        Out-of-plane stress in X-Z plane [Pa]
-    tyz :
-        Out-of-plane stress in Y-Z plane [Pa]
-
-    Returns
-    -------
-    :
-        Von Mises combination of stresses (Pa) in a TF coil.
-    """
-    return np.sqrt(
-        0.5
-        * (
-            (sx - sy) ** 2
-            + (sx - sz) ** 2
-            + (sz - sy) ** 2
-            + 6 * (txy**2 + txz**2 + tyz**2)
-        )
-    )
-
-
 @numba.njit(cache=True, error_model="numpy")
 def extended_plane_strain(
     nu_t,
diff --git a/tests/unit/models/engineering/test_materials.py b/tests/unit/models/engineering/test_materials.py
index 6b05f78652..ec3f576900 100644
--- a/tests/unit/models/engineering/test_materials.py
+++ b/tests/unit/models/engineering/test_materials.py
@@ -1,6 +1,11 @@
+import numpy as np
 import pytest
 
-from process.models.engineering.materials import eurofer97_thermal_conductivity
+from process.models.engineering.materials import (
+    calculate_tresca_stress,
+    calculate_von_mises_stress,
+    eurofer97_thermal_conductivity,
+)
 
 
 def test_eurofer97_thermal_conductivity(monkeypatch, process_models):
@@ -9,3 +14,147 @@ def test_eurofer97_thermal_conductivity(monkeypatch, process_models):
     assert eurofer97_thermal_conductivity(
         1900.0, process_models.data.fwbs.fw_th_conductivity
     ) == pytest.approx(326.70406785462256)
+
+
+def test_calculate_tresca_stress_floats():
+    """Test Tresca stress calculation with float inputs."""
+    result = calculate_tresca_stress(100.0, 50.0, 25.0)
+    assert result == pytest.approx(75.0)
+
+
+def test_calculate_tresca_stress_ndarrays():
+    """Test Tresca stress calculation with ndarray inputs."""
+    stress_x = np.array([100.0, 200.0])
+    stress_y = np.array([50.0, 100.0])
+    stress_z = np.array([25.0, 50.0])
+
+    result = calculate_tresca_stress(stress_x, stress_y, stress_z)
+    expected = np.array([75.0, 150.0])
+
+    np.testing.assert_array_almost_equal(result, expected)
+
+
+def test_calculate_tresca_stress_mixed_float_ndarray():
+    """Test Tresca stress calculation with mixed float and ndarray inputs."""
+    stress_x = 100.0
+    stress_y = np.array([50.0, 75.0])
+    stress_z = 25.0
+
+    result = calculate_tresca_stress(stress_x, stress_y, stress_z)
+    expected = np.array([75.0, 75.0])
+
+    np.testing.assert_array_almost_equal(result, expected)
+
+
+def test_calculate_tresca_stress_zeros():
+    """Test Tresca stress calculation with zero values."""
+    result = calculate_tresca_stress(0.0, 0.0, 0.0)
+    assert result == pytest.approx(0.0)
+
+
+def test_calculate_tresca_stress_zeros_ndarray():
+    """Test Tresca stress calculation with zero ndarray values."""
+    stress_x = np.array([0.0, 100.0])
+    stress_y = np.array([0.0, 50.0])
+    stress_z = np.array([0.0, 25.0])
+
+    result = calculate_tresca_stress(stress_x, stress_y, stress_z)
+    expected = np.array([0.0, 75.0])
+
+    np.testing.assert_array_almost_equal(result, expected)
+
+
+def test_calculate_von_mises_stress_floats():
+    """Test von Mises stress calculation with float inputs."""
+    result = calculate_von_mises_stress(100.0, 50.0, 25.0, 10.0, 5.0, 2.0)
+    expected = np.sqrt(
+        0.5
+        * (
+            (100.0 - 50.0) ** 2
+            + (50.0 - 25.0) ** 2
+            + (25.0 - 100.0) ** 2
+            + 6 * (10.0**2 + 5.0**2 + 2.0**2)
+        )
+    )
+    assert result == pytest.approx(expected)
+
+
+def test_calculate_von_mises_stress_ndarrays():
+    """Test von Mises stress calculation with ndarray inputs."""
+    stress_x = np.array([100.0, 200.0])
+    stress_y = np.array([50.0, 100.0])
+    stress_z = np.array([25.0, 50.0])
+    stress_shear_xy = np.array([10.0, 20.0])
+    stress_shear_yz = np.array([5.0, 10.0])
+    stress_shear_zx = np.array([2.0, 4.0])
+
+    result = calculate_von_mises_stress(
+        stress_x, stress_y, stress_z, stress_shear_xy, stress_shear_yz, stress_shear_zx
+    )
+    expected = np.sqrt(
+        0.5
+        * (
+            (stress_x - stress_y) ** 2
+            + (stress_y - stress_z) ** 2
+            + (stress_z - stress_x) ** 2
+            + 6 * (stress_shear_xy**2 + stress_shear_yz**2 + stress_shear_zx**2)
+        )
+    )
+
+    np.testing.assert_array_almost_equal(result, expected)
+
+
+def test_calculate_von_mises_stress_mixed_float_ndarray():
+    """Test von Mises stress calculation with mixed float and ndarray inputs."""
+    stress_x = 100.0
+    stress_y = np.array([50.0, 75.0])
+    stress_z = 25.0
+    stress_shear_xy = 10.0
+    stress_shear_yz = np.array([5.0, 7.5])
+    stress_shear_zx = 2.0
+
+    result = calculate_von_mises_stress(
+        stress_x, stress_y, stress_z, stress_shear_xy, stress_shear_yz, stress_shear_zx
+    )
+    expected = np.sqrt(
+        0.5
+        * (
+            (stress_x - stress_y) ** 2
+            + (stress_y - stress_z) ** 2
+            + (stress_z - stress_x) ** 2
+            + 6 * (stress_shear_xy**2 + stress_shear_yz**2 + stress_shear_zx**2)
+        )
+    )
+
+    np.testing.assert_array_almost_equal(result, expected)
+
+
+def test_calculate_von_mises_stress_zeros():
+    """Test von Mises stress calculation with zero values."""
+    result = calculate_von_mises_stress(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+    assert result == pytest.approx(0.0)
+
+
+def test_calculate_von_mises_stress_zeros_ndarray():
+    """Test von Mises stress calculation with zero ndarray values."""
+    stress_x = np.array([0.0, 100.0])
+    stress_y = np.array([0.0, 50.0])
+    stress_z = np.array([0.0, 25.0])
+    stress_shear_xy = np.array([0.0, 10.0])
+    stress_shear_yz = np.array([0.0, 5.0])
+    stress_shear_zx = np.array([0.0, 2.0])
+
+    result = calculate_von_mises_stress(
+        stress_x, stress_y, stress_z, stress_shear_xy, stress_shear_yz, stress_shear_zx
+    )
+    expected = np.sqrt(
+        0.5
+        * (
+            (stress_x - stress_y) ** 2
+            + (stress_y - stress_z) ** 2
+            + (stress_z - stress_x) ** 2
+            + 6 * (stress_shear_xy**2 + stress_shear_yz**2 + stress_shear_zx**2)
+        )
+    )
+
+    np.testing.assert_array_almost_equal(result, expected)
diff --git a/tests/unit/models/tfcoil/test_tfcoil.py b/tests/unit/models/tfcoil/test_tfcoil.py
index da98774ee6..7e314c5e55 100644
--- a/tests/unit/models/tfcoil/test_tfcoil.py
+++ b/tests/unit/models/tfcoil/test_tfcoil.py
@@ -1950,21 +1950,6 @@ def test_eyoung_parallel_array(eyoungparallelarrayparam, monkeypatch):
     )
 
 
-@pytest.mark.parametrize(
-    ("sx", "sy", "sz", "expected"),
-    [
-        (0, -3.2e8, 2.4e8, 486621002.42385757),
-        (-2.8e8, 0, 2.4e8, 450777106.7833858),
-    ],
-)
-def test_sigvm(sx, sy, sz, expected):
-    # could not find an example of a use in PROCESS where
-    # tx, ty, or tz were anything other than 0
-    ret = tfcoil_module.sigvm(sx, sy, sz, 0, 0, 0)
-
-    assert ret == pytest.approx(expected)
-
-
 @pytest.mark.parametrize(
     (
         "ind_tf_coil",