The nanodrum 2: force-distance curves in small gradient approximation

Necessary functions to determine the force

hpre2[ρ_, h0_, Rp_, Rt_, λ_, 
        c_] := (2 (c^2) (λ^3) BesselI[0, (c/λ)] + 
            2 (c^2) (λ^3) BesselI[0, (Rp/λ)] - 
            4 (c^2) (λ^3) BesselI[0, (ρ/λ)] - 
            (c^3) (λ^2) BesselI[1, (c/λ)] + 
            2 c h0 Rt (λ^2) BesselI[1, (c/λ)] + 
            (c^2) Rp (λ^2) BesselI[1, (Rp/λ)] - 
            2 h0 Rp Rt (λ^2) BesselI[1, (Rp/λ)] + 
            c λ BesselI[1, 
                (c/λ)] (-2 (c^2) λ BesselI[0, 
                      (c/λ)] + 
                  2 (c^2) λ BesselI[
                      0, (ρ/λ)] + ((c^2) - 
                        2 h0 Rt) (c BesselI[1, (c/λ)] - 
                        Rp BesselI[1, (Rp/λ)])) BesselK[0, 
                (Rp/λ)] + 
            2 (c^3) (λ^2) BesselI[1, (c/λ)] Log[
                (c/λ)] - 
            2 (c^2) Rp (λ^2) BesselI[1, (Rp/λ)] Log[
                (c/λ)] - 
            2 (c^2) Rp (λ^2) BesselI[1, (Rp/λ)] Log[
                (Rp/λ)] + λ (c BesselI[1, (c/λ)] - 
                  Rp BesselI[1, (Rp/λ)]) BesselK[
                0, (ρ/λ)] (2 (c^2) λ BesselI[0, 
                      (c/λ)] - 
                  2 (c^2) λ BesselI[0, 
                      (Rp/λ)] - ((c^3) - 2 c h0 Rt) BesselI[1, 
                      (c/λ)] + 
                  Rp BesselI[1, 
                      (Rp/λ)] ((c^2) - 2 h0 Rt + 
                        2 (c^2) (-Log[(c/λ)] + 
                              Log[(Rp/λ)]))) + 
            Rp BesselK[1, 
                (Rp/λ)] (2 (c^2) (λ^2) BesselI[0, 
                      (Rp/λ)] BesselI[
                      0, (ρ/λ)] + (c BesselI[1, 
                            (c/λ)] - 
                        Rp BesselI[1, 
                            (Rp/λ)]) (λ BesselI[
                            0, (ρ/λ)] ((c^2) - 2 h0 Rt + 
                              2 (c^2) (-Log[(c/λ)] + 
                                    Log[(Rp/λ)])) + 
                        c ((c^2) - 2 h0 Rt) BesselI[1, 
                            (c/λ)] (Log[(Rp/λ)] - 
                              Log[(ρ/λ)])) + 
                  2 (c^2) λ BesselI[0, 
                      (c/λ)] (-λ BesselI[0, 
                            (Rp/λ)] + (-c BesselI[1, 
                                  (c/λ)] + 
                              Rp BesselI[1, (Rp/λ)]) (Log[
                                (Rp/λ)] - 
                              Log[(ρ/λ)]))) + 
            2 (c^2) Rp λ BesselI[1, (Rp/λ)] BesselK[0, 
                (c/λ)] (-λ BesselI[0, 
                      (Rp/λ)] + λ BesselI[
                      0, (ρ/λ)] + 
                  Rp BesselI[1, 
                      (Rp/λ)] (Log[(Rp/λ)] - 
                        Log[(ρ/λ)])) - 
            2 (c^3) (λ^2) BesselI[1, 
                (c/λ)] Log[(ρ/λ)] + 
            4 (c^2) Rp (λ^2) BesselI[1, 
                (Rp/λ)] Log[(ρ/λ)] + 
            c (-2 (c^2) λ BesselI[0, 
                      (c/λ)] + ((c^2) - 
                        2 h0 Rt) (c BesselI[1, (c/λ)] - 
                        Rp BesselI[1, (Rp/λ)])) BesselK[1, 
                (c/λ)] (λ BesselI[0, 
                      (Rp/λ)] - λ BesselI[
                      0, (ρ/λ)] + 
                  Rp BesselI[1, 
                      (Rp/λ)] (-Log[(Rp/λ)] + 
                        Log[(ρ/λ)])))/(2 Rt (c BesselI[1, (c/λ)] - 
                Rp BesselI[1, 
                    (Rp/λ)]) (λ (-2 λ + 
                      c BesselI[0, (Rp/λ)] BesselK[1, 
                          (c/λ)] + 
                      Rp BesselI[0, (c/λ)] BesselK[1, 
                          (Rp/λ)]) + 
                Rp BesselI[1, 
                    (Rp/λ)] (λ BesselK[0, (c/λ)] + 
                      c BesselK[1, 
                          (c/λ)] (Log[(c/λ)] - 
                            Log[(Rp/λ)])) + 
                c BesselI[1, 
                    (c/λ)] (λ BesselK[0, (Rp/λ)] + 
                      Rp BesselK[1, 
                          (Rp/λ)] (-Log[(c/λ)] + 
                            Log[(Rp/λ)]))))

FullSimplify[D[hpre2[ρ, h0, Rp, Rt, λ, c], ρ], 
      Assumptions -> {ρ > 0, Rp > 0, Rt > 0, 
          Rp > Rt, λ > 0, c > 0}]

hpre2der[ρ_, h0_, Rp_, Rt_, λ_, 
        c_] := (-((2 (c^3) (λ^2) BesselI[1, 
                      (c/λ)])/ρ) + ((4 (c^2) Rp (λ^2) BesselI[1, (Rp/λ)])/ρ) - 
            4 (c^2) (λ^2) BesselI[1, (ρ/λ)] + 
            2 (c^2) Rp λ BesselI[1, 
                (Rp/λ)] (-((Rp BesselI[1, 
                            (Rp/λ)])/ρ) + 
                  BesselI[1, (ρ/λ)]) BesselK[0, 
                (c/λ)] + 
            2 (c^3) λ BesselI[1, (c/λ)] BesselI[
                1, (ρ/λ)] BesselK[0, (Rp/λ)] + 
            c (-2 (c^2) λ BesselI[0, 
                      (c/λ)] + ((c^2) - 
                        2 h0 Rt) (c BesselI[1, (c/λ)] - 
                        Rp BesselI[1, 
                            (Rp/λ)])) (((Rp BesselI[1, (Rp/λ)])/ρ) - BesselI[1, (ρ/λ)]) BesselK[1, 
                (c/λ)] + 
            Rp BesselK[1, 
                (Rp/λ)] (((2 (c^2) λ BesselI[0, (c/λ)] (c BesselI[1, (c/λ)] - Rp BesselI[1, (Rp/λ)]))/ρ) + 2 (c^2) λ BesselI[0, (Rp/λ)] BesselI[
                      1, (ρ/λ)] - (((c BesselI[1, (c/λ)] - Rp BesselI[1, (Rp/λ)]) (((c^3) - 2 c h0 Rt) BesselI[1, (c/λ)] + ρ BesselI[1, (ρ/λ)] (-(c^2) + 2 h0 Rt + 2 (c^2) Log[(c/Rp)])))/ρ)) - (c BesselI[1, (c/λ)] - Rp BesselI[1, (Rp/λ)]) BesselK[
                1, (ρ/λ)] (2 (c^2) λ BesselI[0, 
                      (c/λ)] - 
                  2 (c^2) λ BesselI[0, (Rp/λ)] - 
                  c ((c^2) - 2 h0 Rt) BesselI[1, (c/λ)] + 
                  Rp BesselI[1, 
                      (Rp/λ)] ((c^2) - 2 h0 Rt + 
                        2 (c^2) Log[(Rp/c)])))/(2 Rt (c BesselI[
                    1, (c/λ)] - 
                Rp BesselI[1, 
                    (Rp/λ)]) (λ (-2 λ + 
                      c BesselI[0, (Rp/λ)] BesselK[1, 
                          (c/λ)] + 
                      Rp BesselI[0, (c/λ)] BesselK[1, 
                          (Rp/λ)]) + 
                Rp BesselI[1, 
                    (Rp/λ)] (λ BesselK[0, (c/λ)] + 
                      c BesselK[1, (c/λ)] Log[(c/Rp)]) + 
                c BesselI[1, 
                    (c/λ)] (λ BesselK[0, (Rp/λ)] + 
                      Rp BesselK[1, (Rp/λ)] Log[
                          (Rp/c)])))

FullSimplify[
      D[D[hpre2[ρ, h0, Rp, Rt, λ, c], ρ], ρ], 
      Assumptions -> {ρ > 0, Rp > 0, Rt > 0, 
          Rp > Rt, λ > 0, c > 0}]

hpre2derder[ρ_, h0_, Rp_, Rt_, λ_, 
        c_] := (((2 (c^3) (λ^2) BesselI[1, (c/λ)])/(ρ^2)) - ((4 (c^2) Rp (λ^2) BesselI[1, (Rp/λ)])/(ρ^2)) - 
            2 (c^2) λ (BesselI[0, (ρ/λ)] + 
                  BesselI[
                    2, (ρ/λ)]) + (((c^2) Rp BesselI[1, (Rp/λ)] (2 Rp λ BesselI[1, (Rp/λ)] + (ρ^2) (BesselI[0, (ρ/λ)] + BesselI[2, (ρ/λ)])) BesselK[0, (c/λ)])/(ρ^2)) + 
            (c^3) BesselI[1, 
                (c/λ)] (BesselI[0, (ρ/λ)] + 
                  BesselI[2, (ρ/λ)]) BesselK[0, 
                (Rp/λ)] + 
            c (-2 (c^2) λ BesselI[0, 
                      (c/λ)] + ((c^2) - 
                        2 h0 Rt) (c BesselI[1, (c/λ)] - 
                        Rp BesselI[1, 
                            (Rp/λ)])) (-((Rp BesselI[1, 
                            (Rp/λ)])/(ρ^2)) - ((BesselI[0, (ρ/λ)] + BesselI[2, (ρ/λ)])/(2 λ))) BesselK[1, (c/λ)] + 
            Rp BesselK[1, 
                (Rp/λ)] ((c^2) BesselI[0, 
                      (Rp/λ)] (BesselI[0, (ρ/λ)] + 
                        BesselI[
                          2, (ρ/λ)]) + (1/(λ (ρ^2))) ((c BesselI[1, (c/λ)] - 
                            Rp BesselI[1, 
                                (Rp/λ)]) (-2 (c^2) (λ^2) BesselI[0, (c/λ)] + ((c^2) - 
                                  2 h0 Rt) ((ρ^2) BesselI[
                                      0, (ρ/λ)] + 
                                  c λ BesselI[1, 
                                      (c/λ)] - λ ρ BesselI[1, (ρ/λ)]) + 
                            2 (c^2) ρ (-ρ BesselI[
                                      0, (ρ/λ)] + λ BesselI[1, (ρ/λ)]) Log[
                                (c/Rp)]))) + (1/(2 λ)) ((c BesselI[1, (c/λ)] - 
                      Rp BesselI[1, (Rp/λ)]) (BesselK[
                        0, (ρ/λ)] + 
                      BesselK[
                        2, (ρ/λ)]) (2 (c^2) λ BesselI[0, (c/λ)] - 2 (c^2) λ BesselI[0, (Rp/λ)] - 
                      c ((c^2) - 2 h0 Rt) BesselI[1, (c/λ)] + 
                      Rp BesselI[1, 
                          (Rp/λ)] ((c^2) - 2 h0 Rt + 
                            2 (c^2) Log[
                                (Rp/c)]))))/(2 Rt (c BesselI[1, 
                    (c/λ)] - 
                Rp BesselI[1, 
                    (Rp/λ)]) (λ (-2 λ + 
                      c BesselI[0, (Rp/λ)] BesselK[1, 
                          (c/λ)] + 
                      Rp BesselI[0, (c/λ)] BesselK[1, 
                          (Rp/λ)]) + 
                Rp BesselI[1, 
                    (Rp/λ)] (λ BesselK[0, (c/λ)] + 
                      c BesselK[1, (c/λ)] Log[(c/Rp)]) + 
                c BesselI[1, 
                    (c/λ)] (λ BesselK[0, (Rp/λ)] + 
                      Rp BesselK[1, (Rp/λ)] Log[
                          (Rp/c)])))

FullSimplify[
      D[D[D[hpre2[ρ, h0, Rp, Rt, λ, 
              c], ρ], ρ], ρ], 
      Assumptions -> {ρ > 0, Rp > 0, Rt > 0, 
          Rp > Rt, λ > 0, c > 0}]

hpre2derderder[ρ_, h0_, Rp_, Rt_, λ_, 
        c_] := (-((16 (c^3) (λ^2) BesselI[1, 
                      (c/λ)])/(ρ^3)) + ((32 (c^2) Rp (λ^2) BesselI[1, (Rp/λ)])/(ρ^3)) - 
            4 (c^2) (3 BesselI[1, (ρ/λ)] + 
                  BesselI[
                    3, (ρ/λ)]) + ((2 (c^2) Rp BesselI[1, (Rp/λ)] (-8 Rp (λ^2) BesselI[1, (Rp/λ)] + (ρ^3) (3 BesselI[1, (ρ/λ)] + BesselI[3, (ρ/λ)])) BesselK[0, (c/λ)])/(λ (ρ^3))) + ((2 (c^3) BesselI[1, (c/λ)] (3 BesselI[1, (ρ/λ)] + BesselI[3, (ρ/λ)]) BesselK[0, (Rp/λ)])/λ) + ((c (-2 (c^2) λ BesselI[0, (c/λ)] + ((c^2) - 2 h0 Rt) (c BesselI[1, (c/λ)] - Rp BesselI[1, (Rp/λ)])) (8 Rp (λ^2) BesselI[1, (Rp/λ)] - (ρ^3) (3 BesselI[1, (ρ/λ)] + BesselI[3, (ρ/λ)])) BesselK[1, (c/λ)])/((λ^2) (ρ^3))) + (1/(λ^2))((c BesselI[1, (c/λ)] - 
                    Rp BesselI[1, (Rp/λ)]) (3 BesselK[
                        1, (ρ/λ)] + 
                    BesselK[
                      3, (ρ/λ)]) (-2 (c^2) λ BesselI[0, (c/λ)] + 2 (c^2) λ BesselI[0, (Rp/λ)] + 
                    c ((c^2) - 2 h0 Rt) BesselI[1, (c/λ)] + 
                    Rp BesselI[1, 
                        (Rp/λ)] (-(c^2) + 2 h0 Rt + 
                          2 (c^2) Log[(c/Rp)]))) + 
            4 Rp BesselK[1, 
                (Rp/λ)] (((4 (c^2) λ BesselI[0, (c/λ)] (c BesselI[1, (c/λ)] - Rp BesselI[1, (Rp/λ)]))/(ρ^3)) + (((c^2) BesselI[0, (Rp/λ)] (3 BesselI[1, (ρ/λ)] + BesselI[3, (ρ/λ)]))/(2 λ)) + (((c BesselI[1, (c/λ)] - Rp BesselI[1, (Rp/λ)]) (-2 c ((c^2) - 2 h0 Rt) (λ^2) BesselI[1, (c/λ)] + (ρ^2) (ρ BesselI[1, (ρ/λ)] - λ BesselI[2, (ρ/λ)]) ((c^2) - 2 h0 Rt - 2 (c^2) Log[(c/Rp)])))/((λ^2) (ρ^3)))))/(8 Rt (c BesselI[1, (c/λ)] - 
                Rp BesselI[1, 
                    (Rp/λ)]) (λ (-2 λ + 
                      c BesselI[0, (Rp/λ)] BesselK[1, 
                          (c/λ)] + 
                      Rp BesselI[0, (c/λ)] BesselK[1, 
                          (Rp/λ)]) + 
                Rp BesselI[1, 
                    (Rp/λ)] (λ BesselK[0, (c/λ)] + 
                      c BesselK[1, (c/λ)] Log[(c/Rp)]) + 
                c BesselI[1, 
                    (c/λ)] (λ BesselK[0, (Rp/λ)] + 
                      Rp BesselK[1, (Rp/λ)] Log[
                          (Rp/c)])))

findcontact[h0_, Rp_?NumberQ, Rt_, λ_, κ_, w_] := 
      Part[Part[
          FindRoot[(hpre2derder[c, h0, Rp, Rt, λ, 
                  c] - ((1/Rt) - (Sqrt[((2*w)/κ)]))), {c, 
              0.000001, 0, Rp}], 1], 2]

findcontact[3, 33.5, 12, 4.69, 1.1, 0.0]

Force via the surface stress tensor

force[h0f_, Rpf_, Rtf_, λf_, κf_, wf_, Rintf_] := 
      Block[{hder, hderder, hderderder, g, 
          cont = findcontact[h0f, Rpf, Rtf, λf, κf, wf]}, 
        hder = hpre2der[Rintf, h0f, Rpf, Rtf, λf, cont]; 
        hderder = hpre2derder[Rintf, h0f, Rpf, Rtf, λf, cont]; 
        hderderder = hpre2derderder[Rintf, h0f, Rpf, Rtf, λf, cont]; 
        g = 1 + (hder^2); ((-2*π*
              Rintf*(((1/2)*κf*((((hderder^2)/(g^3))) - (((hder^2)/((Rintf^2)*g)))) - κ(f/(((λf))^2)))*
                    (hder/(Sqrt[g])) + κf*(-(hder/((Rintf^2)*(Sqrt[g]))) - (((hder^2)*hderder)/(Rintf*(g^(3/2)))) + 
                        (hderder/(Rintf*(Sqrt[g]))) - ((3*hder*(hderder^2))/(g^(5/2))) + 
                        (hderderder/(g^(3/2))))*(1/g))))]

Check how the force depends on r (ideally not at all--errors are due to the small gradient approximation)

Plot[force[3, 33.5, 12, 4.69, 1.1, 0.0, x], {x, 12, 33.5}, 
        AxesLabel ->  {"Rint / nm" ,  "F / nN"}, 
        Background -> GrayLevel[1.0], 
        ImageSize -> 400];

Plot[force[10, 33.5, 12, 4.69, 1.1, 0.0, x], {x, 12, 33.5}, 
        AxesLabel ->  {"Rint / nm" ,  "F / nN"}, 
        Background -> GrayLevel[1.0], 
        ImageSize -> 400];

Force as a function of λ, Rp, and Rt, respectively

Force as a function of λ for h0 = 3 nm

Plot[force[3, 33.5, 12, L, 1.1, 0.0, 33], {L, 0.5, 6}, 
        AxesLabel ->  {"lambda / nm" ,  "force / nN"}, 
        Background -> GrayLevel[1.0], 
        ImageSize -> 400];

Force as a function of the tip radius Rt for h0 = 3 nm

Plot[force[3, 33.5, x, 4.69, 1.1, 0.0, 33], {x, 6, 20}, 
        AxesLabel ->  {"tip radius / nm" ,  "force / nN"}, 
        PlotRange  ->  {0, 0.8}, Background -> GrayLevel[1.0], 
        ImageSize -> 400];

Plot[force[3, 33.5, x, 4.69, 1.1, 0.0, 33], {x, 6, 20}, 
        AxesLabel ->  {"tip radius / nm" ,  "force / nN"}, 
        Background -> GrayLevel[1.0], 
        ImageSize -> 400];

Force as a function of the pore radius Rp for h0 = 3 nm

Plot[force[3, x, 12, 4.69, 1.1, 0.0, 19], {x, 20, 100}, 
        AxesLabel ->  {"pore radius / nm" ,  "force / nN"}, 
        PlotRange  ->  {0, 1}, Background -> GrayLevel[1.0], 
        ImageSize -> 400];

Force-indentation curves -- fitting an experiment

Plot[{0.0142825 - 
            Check[(force[10, 225, 20, λ, 1.1, 0.0, 220]) - (force[5, 225, 20, λ, 1.1, 0.0, 220]), -1]/5}, {λ, 
          13, 14}, PlotRange -> {-0.005, 0.005}, 
        Background -> GrayLevel[1.0]];

Plot[{0.0036 - 
            Check[(force[10, 225, 20, λ, 0.025, 0.0, 220]) - (force[5, 225, 20, λ, 0.025, 0.0, 220]), -1]/5}, {λ, 
          3, 4}, PlotRange -> {-0.0005, 0.0005}, 
        Background -> GrayLevel[1.0]];

PP1 = 
        Plot[(0.0142825*xxx), {xxx, 1, 20}, PlotStyle -> {Hue[0.6]}, 
          DisplayFunction -> Identity];
PP2 = 
        Plot[Check[(force[h1, 225, 20, 13.5, 1.1, 0.0, 
                220]), -1], {h1, 5, 25}, 
          DisplayFunction -> Identity];
PP3 = 
        Plot[(0.00363496*xxx), {xxx, 1, 20}, 
          PlotStyle -> {Hue[0.6]}, 
          DisplayFunction -> Identity];
PP4 = 
        Plot[Check[(force[h2, 225, 20, 3.3, 0.025, 0.0, 
                220]), -1], {h2, 5, 25}, 
          DisplayFunction -> Identity];

Show[PP1, PP2, DisplayFunction -> $DisplayFunction, 
        ImageSize -> 600,  PlotRange  ->  {0, 0.5}, 
        AxesLabel ->  {"indentation / nm" ,  "force / nN"}, 
        Background -> GrayLevel[1.0]];

Show[PP3, PP4, DisplayFunction -> $DisplayFunction, 
        ImageSize -> 600,  PlotRange  ->  {0, 0.1}, 
        AxesLabel ->  {"indentation / nm" ,  "force / nN"}, 
        Background -> GrayLevel[1.0]];