From f0ca87cc5a0705588a08e05cc5a091b64a8a7aa0 Mon Sep 17 00:00:00 2001 From: Christoph Lehmann Date: Wed, 8 Aug 2018 14:20:25 +0200 Subject: [PATCH] Documented Hertz Contact test --- .../Mechanics/Linear/Python/contact_radii.png | 3 - Tests/Data/Mechanics/Linear/Python/post.py | 44 +- .../Linear/Python/stress_at_contact.png | 3 - .../python-bc/hertz-contact/contact_radii.png | 3 + .../python-bc/hertz-contact/hertz-contact.md | 65 +++ .../python-bc/hertz-contact/hertz-contact.png | 3 + .../python-bc/hertz-contact/hertz-contact.svg | 382 ++++++++++++++++++ .../hertz-contact/stress_at_contact.png | 3 + 8 files changed, 485 insertions(+), 21 deletions(-) delete mode 100644 Tests/Data/Mechanics/Linear/Python/contact_radii.png delete mode 100644 Tests/Data/Mechanics/Linear/Python/stress_at_contact.png create mode 100644 web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png create mode 100644 web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md create mode 100644 web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png create mode 100644 web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg create mode 100644 web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png diff --git a/Tests/Data/Mechanics/Linear/Python/contact_radii.png b/Tests/Data/Mechanics/Linear/Python/contact_radii.png deleted file mode 100644 index d46b3899af4..00000000000 --- a/Tests/Data/Mechanics/Linear/Python/contact_radii.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:1866ea736d333d38fefe1cbba5c154b186b698b658e0b1a472211f87cba9f8b0 -size 20673 diff --git a/Tests/Data/Mechanics/Linear/Python/post.py b/Tests/Data/Mechanics/Linear/Python/post.py index 16b3630e53a..341ae9772bf 100755 --- a/Tests/Data/Mechanics/Linear/Python/post.py +++ b/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) diff --git a/Tests/Data/Mechanics/Linear/Python/stress_at_contact.png b/Tests/Data/Mechanics/Linear/Python/stress_at_contact.png deleted file mode 100644 index b5574cc9c83..00000000000 --- a/Tests/Data/Mechanics/Linear/Python/stress_at_contact.png +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3a29c1985003d0e0a62244d144c17dc7041f4e5124b58e26639b7f2f21b4b2fa -size 109883 diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png b/web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png new file mode 100644 index 00000000000..9fbace0274c --- /dev/null +++ b/web/content/docs/benchmarks/python-bc/hertz-contact/contact_radii.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:28b7641fdaea7c0b21428166d7e7468244cfc806f27d6d49221965f600323048 +size 22268 diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.md new file mode 100644 index 00000000000..6657f8bd5ba --- /dev/null +++ b/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" >}} diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png new file mode 100644 index 00000000000..1a80d39cbdf --- /dev/null +++ b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:02b46de2aa5e44e3fa45cde5c6a7fdf219ccfad8bafe3b78c2cf551303ea7c56 +size 29066 diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg new file mode 100644 index 00000000000..07dc91e4480 --- /dev/null +++ b/web/content/docs/benchmarks/python-bc/hertz-contact/hertz-contact.svg @@ -0,0 +1,382 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + w0 + + + a + contact area + + + simulation domain + + + R + R + + ΞΎ + + + diff --git a/web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png b/web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png new file mode 100644 index 00000000000..7ccc1781992 --- /dev/null +++ b/web/content/docs/benchmarks/python-bc/hertz-contact/stress_at_contact.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f90a55eb4d0f757cc46ea743fde8f29394bad228a802819e0e2d8ec0c14eb91e +size 86013