diff --git a/src/fatpy/core/damage_cumulation/damage_cumulation_palmgren_meiner.py b/src/fatpy/core/damage_cumulation/damage_cumulation_palmgren_meiner.py
new file mode 100644
index 0000000..c2330d1
--- /dev/null
+++ b/src/fatpy/core/damage_cumulation/damage_cumulation_palmgren_meiner.py
@@ -0,0 +1,88 @@
+"""Damage Accumulation Rule according to Palmgren-Miner.
+
+Resources:
+    [1] Graphical interpretation of the change in the S-N curve based on the
+        chosen version can be found e.g. in this open-access paper:
+        http://dx.doi.org/10.5545/sv-jme.2013.1348
+    [2] Miner, M. A. (1945). Cumulative damage in fatigue. Journal of Applied
+        Mechanics, 12(3), 159-164.
+"""
+
+# import numpy as np
+# from numpy.typing import NDArray
+
+
+def damage_cumulation_elementary(
+    slope_k: float,
+    constant: float,
+    sig: float,
+    number_occurrences: int,
+) -> float:
+    r"""Elementary version of Palmgren-Miner linear damage accumulation.
+
+    The same slope k of the S-N curve below and above the fatigue limit.
+
+    ??? abstract "Math Equations"
+        $$
+        D = n/N = n\,\frac{\sigma^k}{C}
+        $$
+
+    """
+    total_occurrences: float = constant / sig**slope_k
+
+    damage: float = number_occurrences / total_occurrences
+
+    return damage
+
+
+def damage_cumulation_basic(
+    slope_k: float,
+    constant: float,
+    sig_fl: float,
+    sig: float,
+    number_occurrences: int,
+) -> float:
+    """Basic version of Palmgren-Miner linear damage accumulation.
+
+    The S-N curve gets horizontal at the fatigue limit, no damage for stresses beneath.
+    Otherwise elementary damage is calculated.
+    """
+    if sig < sig_fl:
+        damage = 0.0
+    else:
+        damage = damage_cumulation_elementary(
+            slope_k, constant, sig, number_occurrences
+        )
+
+    return damage
+
+
+def damage_cumulation_haibach(
+    slope_k: float,
+    constant: float,
+    sig_fl: float,
+    sig: float,
+    number_occurrences: int,
+) -> float:
+    r"""Haibach version of Palmgren-Miner linear damage accumulation.
+
+    the original slope_k is modified below fatigue limit to 2*slope_k-1.
+
+    ??? abstract "Math Equations"
+        $$
+        D = \frac{n}{C}\,\frac{\sigma^{2k-1}}{\sigma_\mathrm{FL}^{k-1}}
+        $$
+
+    """
+    if sig < sig_fl:
+        damage: float = (
+            number_occurrences
+            * sig ** (2 * slope_k - 1)
+            / (constant * sig_fl ** (slope_k - 1))
+        )
+    else:
+        damage = damage_cumulation_elementary(
+            slope_k, constant, sig, number_occurrences
+        )
+
+    return damage
diff --git a/tests/core/test_src/fatpy/core/damage_cumulation/test_damage_cumulation.py b/tests/core/test_src/fatpy/core/damage_cumulation/test_damage_cumulation.py
new file mode 100644
index 0000000..b40e421
--- /dev/null
+++ b/tests/core/test_src/fatpy/core/damage_cumulation/test_damage_cumulation.py
@@ -0,0 +1,105 @@
+"""Test functions for damage accumulation rules."""
+
+import pytest
+import numpy as np
+# from numpy.typing import NDArray
+
+from fatpy.core.damage_cumulation import damage_cumulation_palmgren_meiner as dcpm
+
+
+@pytest.fixture
+def damage_cumulation_parameters() -> dict[str, float]:
+    """Fixture providing parameters for damage cumulation tests.
+
+    Returns:
+        dict[str, float]: Parameters including slope_k, constant, sig_fl.
+    """
+    params = {
+        "slope_k": 5.0,
+        "constant": 1e17,  # 1e15 is on the internet
+        "sig_fl": 137.97,
+    }
+    return params
+
+
+@pytest.fixture
+def fatigue_load_low() -> tuple[float, int]:
+    """Fixture providing a sample fatigue load.
+
+    Returns:
+        tuple[float, int]: Sample stress and number of occurrences.
+    """
+    return 150.0, 5000
+
+
+@pytest.fixture
+def fatigue_load_hi() -> tuple[float, int]:
+    """Fixture providing a sample fatigue load.
+
+    Returns:
+        tuple[float, int]: Sample stress and number of occurrences.
+    """
+    return 110.0, 100000
+
+
+def test_damage_cumulation_elementary(
+    damage_cumulation_parameters: dict[str, float],
+    fatigue_load_low: tuple[float, int],
+    fatigue_load_hi: tuple[float, int],
+) -> None:
+    """Elementary version
+    the same slope k of the S-N curve below and above the fatigue limit
+    """
+    slope_k = damage_cumulation_parameters["slope_k"]
+    constant = damage_cumulation_parameters["constant"]
+    sig_low, n_low = fatigue_load_low
+    sig_hi, n_hi = fatigue_load_hi
+
+    d_low = dcpm.damage_cumulation_elementary(slope_k, constant, sig_low, n_low)
+    d_hi = dcpm.damage_cumulation_elementary(slope_k, constant, sig_hi, n_hi)
+
+    assert np.around(d_low, decimals=4) == 0.0038
+    assert np.around(d_hi, decimals=4) == 0.0161
+
+
+def test_damage_cumulation_basic(
+    damage_cumulation_parameters: dict[str, float],
+    fatigue_load_low: tuple[float, int],
+    fatigue_load_hi: tuple[float, int],
+) -> None:
+    """Basic version
+    the S-N curve gets horizontal at the fatigue limit,
+    no damage for stresses beneath
+    """
+    slope_k = damage_cumulation_parameters["slope_k"]
+    constant = damage_cumulation_parameters["constant"]
+    sig_fl = damage_cumulation_parameters["sig_fl"]
+    sig_low, n_low = fatigue_load_low
+    sig_hi, n_hi = fatigue_load_hi
+
+    d_low = dcpm.damage_cumulation_basic(slope_k, constant, sig_fl, sig_low, n_low)
+    d_hi = dcpm.damage_cumulation_basic(slope_k, constant, sig_fl, sig_hi, n_hi)
+
+    assert np.around(d_low, decimals=4) == 0.0038
+    assert d_hi == 0.0
+
+
+def test_damage_cumulation_haibach(
+    damage_cumulation_parameters: dict[str, float],
+    fatigue_load_low: tuple[float, int],
+    fatigue_load_hi: tuple[float, int],
+) -> None:
+    """Haibach version
+    the original slope_k is modified below fatigue limit to 2*slope_k-1
+    """
+    slope_k = damage_cumulation_parameters["slope_k"]
+    constant = damage_cumulation_parameters["constant"]
+    sig_fl = damage_cumulation_parameters["sig_fl"]
+    sig_low, n_low = fatigue_load_low
+    sig_hi, n_hi = fatigue_load_hi
+
+    d_low = dcpm.damage_cumulation_haibach(slope_k, constant, sig_fl, sig_low, n_low)
+    d_hi = dcpm.damage_cumulation_haibach(slope_k, constant, sig_fl, sig_hi, n_hi)
+
+    assert np.around(d_low, decimals=4) == 0.0038
+    assert np.around(d_hi, decimals=5) == 0.00651