Documented Hertz Contact test

Tests/Data/Mechanics/Linear/Python/post.py

@@ -117,8 +117,8 @@ def get_y_top(t):
 fig, ax = plt.subplots()
 
 fig.subplots_adjust(right=0.75)
-ax.set_xlabel("$r$ / m")
-ax.set_ylabel(r"$\sigma_{yy}$ / Pa")
+ax.set_xlabel(r"$\xi$ / m")
+ax.set_ylabel(r"$\bar p$ / Pa")
 add_leg = True
 
 
@@ -268,11 +268,11 @@ def computeStrain(grid):
 
     def average_stress(rs, stress):
         r_contact = max(rs)
-        rs_int = np.linspace(min(rs), max(rs), max(len(rs), 200))
+        rs_int = np.linspace(min(rs), max(rs)-1e-8, max(len(rs), 200))
         stress_int = interp1d(rs, stress, bounds_error=False, fill_value=0.0)
         avg_stress = np.empty_like(rs_int)
 
-        for i, r in enumerate(rs_int[:-1]):
+        for i, r in enumerate(rs_int):
             rho_max = np.sqrt(r_contact**2 - r**2)
             rhos = np.linspace(0, rho_max, 100)
             xis = np.sqrt(rhos**2 + r**2)
@@ -282,7 +282,7 @@ def average_stress(rs, stress):
                 print(max(xis), r_contact)
                 raise
             avg_stress[i] = 1.0 / rho_max * np.trapz(x=rhos, y=stress_int(xis))
-        avg_stress[-1] = 0.0
+        # avg_stress[-1] = 0.0
 
         return rs_int, avg_stress
 
@@ -304,19 +304,31 @@ def stress_at_contact_area():
             print(min(xs), max(xs), r_contact, max(xs) - r_contact)
             raise
 
+        w_0 = 2.0 * (1.0 - y_top)
+
         rs = np.linspace(0, r_contact, 200)
         if add_leg: ax.plot([], [], color="k", ls=":", label="ref")
         h, = ax.plot(rs, -p_contact(rs, r_contact), ls=":")
 
+        if False:
+            r_contact_ana = np.sqrt(w_0 * R)
+
+            rs2 = np.linspace(0, r_contact_ana, 200)
+            if add_leg: ax.plot([], [], color="k", ls="--", label="ref2")
+            ax.plot(rs2, -p_contact(rs, r_contact_ana), ls="--", color=h.get_color())
+
         stress_yy = vtk_to_numpy(grid.GetPointData().GetArray("sigma"))[:,1]
         rs, avg_stress_yy = average_stress(xs, stress_yy)
-        if add_leg: ax.plot([], [], color="k", ls="--", label="ogs")
-        ax.plot(rs, avg_stress_yy, color=h.get_color(), ls="--")
+        if add_leg: ax.plot([], [], color="k", ls="-", label="ogs")
+        ax.plot(rs, avg_stress_yy, color=h.get_color(), ls="-",
+                label=r"$w_0 = {}$".format(w_0))
 
-        stress = vtk_to_numpy(grid.GetPointData().GetArray("stress_post"))
-        rs, avg_stress_yy = average_stress(xs, stress[:,1])
-        if add_leg: ax.plot([], [], color="k", label="post")
-        ax.plot(rs, avg_stress_yy, color=h.get_color(), label=r"$\Delta w = {}$".format(2*(1.0 - y_top)))
+        if False:
+            stress = vtk_to_numpy(grid.GetPointData().GetArray("stress_post"))
+            rs, avg_stress_yy = average_stress(xs, stress[:,1])
+            if add_leg: ax.plot([], [], color="k", label="post")
+            ax.plot(rs, avg_stress_yy, color=h.get_color(),
+                    label=r"$w_0 = {}$".format(2*(1.0 - y_top)))
 
         ax.scatter([r_contact], [0], color=h.get_color())
         ax.legend(loc="center left", bbox_to_anchor=(1.0, 0.5))
@@ -332,15 +344,17 @@ def stress_at_contact_area():
 
 
 fig, ax = plt.subplots()
-ax.scatter(ws, rs_contact)
+ax.scatter(ws, rs_contact, label="ogs")
 
 ws_ref = np.linspace(0, max(ws), 200)
 rs_ref = np.sqrt(ws_ref * R)
 
-ax.plot(ws_ref, rs_ref)
+ax.plot(ws_ref, rs_ref, label="ref")
+
+ax.legend()
 
-ax.set_xlabel("displacement of sphere centers / m")
-ax.set_ylabel("radius of contact areas / m")
+ax.set_xlabel("displacement of sphere centers $w_0$ / m")
+ax.set_ylabel("radius of contact area $a$ / m")
 
 fig.savefig("contact_radii.png")
 plt.close(fig)

web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md

@@ -0,0 +1,65 @@
++++
+project = "https://github.com/ufz/ogs-data/blob/master/Mechanics/Linear/Python/hertz-contact.prj"
+author = "Christoph Lehmann"
+date = "2018-08-06T11:41:00+02:00"
+title = "Hertz Contact using Python Boundary Conditions"
+weight = 156
+
+[menu]
+  [menu.benchmarks]
+    parent = "python-bc"
+
++++
+
+{{< data-link >}}
+
+## Problem description
+
+Two elastic spheres of same radius $R$ are brought into contact.
+The sphere centers are displaced towards each other by $w\_0$, with increasing
+values in every load step.
+Due to symmetry reasons a flat circular contact area of radius $a$ forms.
+
+{{< img src="../hertz-contact.png" >}}
+
+The contact between the two spheres is modelled as a Dirichlet BC
+on a varying boundary. The exact boundary and Dirichlet values for the
+$y$ displacements are determined in a Python script.
+Compared to the sketch above the prescribed $y$ displacements correspond
+to $w\_0/2$, because due to symmetry only half of the system (a section of the
+lower sphere) is simulated.
+
+
+## Analytical solution
+
+The radius of the contact area is given by
+$$
+\begin{equation}
+a = \sqrt{\frac{w\_0 R}{2}}
+\end{equation}
+$$
+
+The average pressure $\bar p$ over a the secant with distance $\xi$ to the
+center of the contact area (cf. vertical dashed line in the sketch above) is assumed to be
+$$
+\begin{equation}
+\bar p(\xi) = \kappa \sqrt{a^2 - \xi^2}
+\end{equation}
+$$
+with the prefactor $\kappa$ given by
+$$
+\begin{equation}
+\kappa = \frac{G}{R \cdot (1-\nu)}
+\end{equation}
+$$
+
+
+## Results
+
+Contact radii:
+
+{{< img src="../contact_radii.png" >}}
+
+Average pressure $\bar{p}$:
+
+{{< img src="../stress_at_contact.png" >}}