Enh/higher derivatives jkr

Adhesion/ChangeLog.md

@@ -2,16 +2,21 @@
 Change log for Adhesion
 =======================
 
+v0.93.0 (05Jun25)
+-----------------
 
+- ENH: third derivative of stress intensity factor of a JKR contact
 
-v0.92.1 (13May25)
+v0.92.1 (05Jun25)
 -----------------
+
 - MAINT: adjust to new muFFT and NUMPI APIs
 - MAINT: drop python 3.8 support
 - DOC: correct some missing renames
 
 v0.92.0 (14Jan24)
 -----------------
+- 
 - ENH: Analytical formulas for JKR contact 
 - BUG: Fixed discover of version when in a git repository that is not the
   source directory of SurfaceTopography

Adhesion/ReferenceSolutions/JKR.py

@@ -620,6 +620,7 @@ def _stress_intensity_factor_from_penetration_radius(contact_radius, penetration
 
     dh = a ** 2 / R  # hertzian displacement
     dc = d - dh  # displacement imposed to the external crack
+    dcda = - 2 * a / R
 
     if der == "0":
         return - (Es * dc / np.sqrt(np.pi * a))
@@ -628,6 +629,8 @@ def _stress_intensity_factor_from_penetration_radius(contact_radius, penetration
             + Es / np.sqrt(np.pi * a) * 2 * a / R
     elif der == "2_a":
         return - 3 / 4 * a ** (-5 / 2) * Es * dc / np.sqrt(np.pi)
+    elif der == "3_a":
+        return 15 / 8 * a ** (-7 / 2) * Es * dc / np.sqrt(np.pi) - 3 / 4 * a ** (-5 / 2) * Es * dcda / np.sqrt(np.pi)
     elif der == "1_d":
         return - Es / np.sqrt(np.pi * a)
     elif der == "2_ad" or der == "2_da":

test/ReferenceSolutions/test_sphere_jkr.py

@@ -184,6 +184,31 @@ def test_stress_intensity_factor_second_derivative():
                                atol=1e-6, rtol=1e-4)
 
 
+def test_stress_intensity_factor_third_derivative():
+    a = np.linspace(0.5, 1, 1000)
+
+    pen = 0.5
+    am = a[2:-2]
+    da = a[1] - a[0]
+    dK_da3_num = (-1 / 2 * JKR.stress_intensity_factor(a[:-4], pen, )
+                  + JKR.stress_intensity_factor(a[1:-3], pen, )
+                  - JKR.stress_intensity_factor(a[3:-1], pen, )
+                  + 0.5 * JKR.stress_intensity_factor(a[4:], pen, )) / da ** 3
+    dK_da3_analytical = JKR.stress_intensity_factor(
+        contact_radius=am,
+        penetration=pen, der="3_a")
+    if True:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots()
+        ax.plot(am, dK_da3_num, "+", label="numerical")
+        ax.plot(am, dK_da3_analytical, "o", mfc="none",
+                label="analytical")
+        plt.show()
+
+    np.testing.assert_allclose(dK_da3_analytical, dK_da3_num,
+                               atol=1e-6, rtol=1e-4)
+
+
 @pytest.mark.parametrize("pen", (0.1, 0.5, 1))
 def test_stress_intensity_factor_second_derivative_nonadhesive(pen):
     a = JKR.contact_radius(penetration=pen, work_of_adhesion=0)