Scripting examples

Approximating the Hertzian model with a spherical indenter

There is no closed form for the Hertzian model with a spherical indenter. The force \(F\) does not directly depend on the indentation depth \(\delta\), but has an indirect dependency via the radius of the circular contact area between indenter and sample \(a\) [Sne65]:

\[\begin{split}F &= \frac{E}{1-\nu^2} \left( \frac{R^2+a^2}{2} \ln \! \left( \frac{R+a}{R-a}\right) -aR \right) \label{eq:1}\\ \delta &= \frac{a}{2} \ln \! \left(\frac{R+a}{R-a}\right) \label{eq:2}\end{split}\]

Here, \(E\) is the Young’s modulus, \(R\) is the radius of the indenter, and \(\nu\) is the Poisson’s ratio of the probed material.

Because of this indirect dependency, fitting this model to experimental data can be time-consuming. Therefore, it is beneficial to approximate this model with a polynomial function around small values of \(\delta/R\) using the Hertz model for a parabolic indenter as a starting point [Dob18]:

\[F = \frac{4}{3} \frac{E}{1-\nu^2} \sqrt{R} \delta^{3/2} \left(1 - \frac{1}{10} \frac{\delta}{R} - \frac{1}{840} \left(\frac{\delta}{R}\right)^2 + \frac{11}{15120} \left(\frac{\delta}{R}\right)^3 + \frac{1357}{6652800} \left(\frac{\delta}{R}\right)^4 \right)\]

This example illustrates the error made with this approach. In nanite, the model for a spherical indenter has the identifier “sneddon_spher” and the approximate model has the identifier “sneddon_spher_approx”.

The plot shows the error for the parabolic indenter model “hertz_para” and for the approximation to the spherical indenter model. The maximum indentation depth is set to \(R\). The error made by the approximation of the spherical indenter is more than four magnitudes lower than the maximum force during indentation.

_images/model_spherical_indenter.png

model_spherical_indenter.py

 1import matplotlib.pylab as plt
 2from mpl_toolkits.axes_grid1 import make_axes_locatable
 3from matplotlib.lines import Line2D
 4import matplotlib as mpl
 5import numpy as np
 6
 7from nanite.model import models_available
 8
 9# models
10exact = models_available["sneddon_spher"]
11approx = models_available["sneddon_spher_approx"]
12para = models_available["hertz_para"]
13# parameters
14params = exact.get_parameter_defaults()
15params["E"].value = 1000
16
17# radii
18radii = np.linspace(2e-6, 100e-6, 20)
19
20# plot results
21plt.figure(figsize=(8, 5))
22
23# overview plot
24ax = plt.subplot()
25for ii, rad in enumerate(radii):
26    params["R"].value = rad
27    # indentation range
28    x = np.linspace(0, -rad, 300)
29    yex = exact.model(params, x)
30    yap = approx.model(params, x)
31    ypa = para.model(params, x)
32    ax.plot(x*1e6, np.abs(yex - yap)/yex.max(),
33            color=mpl.cm.get_cmap("viridis")(ii/radii.size),
34            zorder=2)
35    ax.plot(x*1e6, np.abs(yex - ypa)/yex.max(), ls="--",
36            color=mpl.cm.get_cmap("viridis")(ii/radii.size),
37            zorder=1)
38
39ax.set_xlabel(r"indentation depth $\delta$ [µm]")
40ax.set_ylabel("error in force relative to maximum $F/F_{max}$")
41ax.set_yscale("log")
42ax.grid()
43
44# legend
45custom_lines = [Line2D([0], [0], color="k", ls="--"),
46                Line2D([0], [0], color="k", ls="-"),
47                ]
48ax.legend(custom_lines, ['parabolic indenter',
49                         'approximation of spherical indenter'])
50
51divider = make_axes_locatable(ax)
52cax = divider.append_axes("right", size="3%", pad=0.05)
53
54norm = mpl.colors.Normalize(vmin=radii[0]*1e6, vmax=radii[-1]*1e6)
55mpl.colorbar.ColorbarBase(ax=cax,
56                          cmap=mpl.cm.viridis,
57                          norm=norm,
58                          orientation='vertical',
59                          label="indenter radius [µm]"
60                          )
61
62plt.tight_layout()
63plt.show()

Fitting and rating

This example uses a force-distance curve of a zebrafish spinal cord section to illustrate basic data fitting and rating with nanite. The dataset is part of a study on spinal cord stiffness in zebrafish [MKH+19].

_images/fit_and_rate.png

fit_and_rate.py

 1import matplotlib.gridspec as gridspec
 2import matplotlib.pylab as plt
 3
 4import nanite
 5
 6# load the data
 7group = nanite.load_group("data/zebrafish-head-section-gray-matter.jpk-force")
 8idnt = group[0]  # this is an instance of `nanite.Indentation`
 9# apply preprocessing
10idnt.apply_preprocessing(["compute_tip_position",
11                          "correct_force_offset",
12                          "correct_tip_offset"])
13# set the fit model ("sneddon_spher_approx" is faster than "sneddon_spher"
14# and sufficiently accurate)
15idnt.fit_properties["model_key"] = "sneddon_spher_approx"
16# get the initial fit parameters
17params = idnt.get_initial_fit_parameters()
18# set the correct indenter radius
19params["R"].value = 18.64e-06
20# perform the fit with the edited parameters
21idnt.fit_model(params_initial=params)
22# obtain the Young's modulus
23emod = idnt.fit_properties["params_fitted"]["E"].value
24# obtain a rating for the dataset
25# (using default regressor and training set)
26rate = idnt.rate_quality()
27
28# overview plot
29plt.figure(figsize=(8, 5))
30gs = gridspec.GridSpec(2, 1, height_ratios=[5, 1])
31
32ax1 = plt.subplot(gs[0])
33ax2 = plt.subplot(gs[1])
34
35# only plot the approach part (`1` would be retract)
36where_approach = idnt["segment"] == 0
37
38# plot force-distance data (nanite uses SI units)
39ax1.plot(idnt["tip position"][where_approach] * 1e6,
40         idnt["force"][where_approach] * 1e9,
41         label="data")
42ax1.plot(idnt["tip position"][where_approach] * 1e6,
43         idnt["fit"][where_approach] * 1e9,
44         label="fit (spherical indenter)")
45ax1.text(.2, 2.05,
46         "apparent Young's modulus: {:.0f} Pa\n".format(emod)
47         + "rating: {:.1f}".format(rate),
48         ha="center")
49ax1.legend()
50# plot resiudals
51ax2.plot(idnt["tip position"][where_approach] * 1e6,
52         (idnt["force"] - idnt["fit"])[where_approach] * 1e9)
53
54# update plot parameters
55ax1.set_xlim(-4.5, 3)
56ax1.set_ylabel("force [pN]")
57ax1.grid()
58ax2.set_xlim(-4.5, 3)
59ax2.set_ylim(-.2, .2)
60ax2.set_ylabel("residuals [pN]")
61ax2.set_xlabel("tip position [µm]")
62ax2.grid()
63
64plt.tight_layout()
65plt.show()