modeling the elastic response of a rough surface using the error function

Question

I'm trying to model the elastic deformation of a part with rough surface using the error function Erf. A linear elastic part with theoretically smooth surface has a deformation of

$$\epsilon\left(u,\sigma\right)=\text{Clip}\left(\frac{\sigma}{E_y} ,0,1\right)\tag{1}$$

But as mentioned in this paper one can represent the effect of surface roughness as variable $E_y$:

We can use the error function to write the elasticity of a surface with normal distribution of asperity altitude:

$$E_y\left(u\right)=\frac{E_{y0}}{2}\left(1+\text{Erf}\left(\frac{u-m}{d\sqrt{2}}\right)\right) \tag{2}$$

If we assume the left hand side of the part with rough surface starts at $u=0$, then $m=0$ and for simplicity we assume $d=1$. We also assume the other side of the part is a smooth surface connected to the frame (no displacement) at point $u=l_2$. Now if we apply a stress $\sigma$ to the left side of the part in $-\infty$ the part will be squeezed to the point $x$ calculated as:

$$x\left( \sigma \right)=\lim_{l_1 \to -\infty}\left( \int_{l_1}^{l_2} \epsilon\left( u,\sigma\right)du +l_1\right)$$

To solve this in Mathematica first:

m = 0;
d = 1;
E0 = 10^9;
l2 = 10;
Ey[u_] := E0*(1 + Erf[(u - m)/(d*Sqrt[2])])/2;
ep[u_, s_] := Clip[s/Ey[u], 0, 1];

Now if I try to calculate the $x\left( s \right)$ analytically:

xa[s_] := Limit[Integrate[ep[u, s], {u, l1, l2}] + l1, l1 -> -Infinity]

Or numerically:

xn[s_] := Limit[With[{l11 = l1}, NIntegrate[ep[u, s], {u, l11, l2}] + l11],  l1 -> -Infinity]

On both cases when I try to calculate for example x[1] I receive the error:

Clip::rtwo: The argument 0 at position 2 is expected to be a list of a lower clip bound and an upper clip bound.

and

NIntegrate::nlim: u = l1 is not a valid limit of integration.

Which is very confusing for me. I would appreciate if you could let me know where is my mistake and I how I can solve it.

Answer 1

In the following I corrected the mistake in Clip and put the + l1 inside the integral:

m = 0;
d = 1;
E0 = 10^9;
l2 = 10;
Ey[u_] := E0*(1 + Erf[(u - m)/(d*Sqrt[2])])/2;
ep[u_, s_] := Clip[s/Ey[u], {0, 1}];
xn[s_, l1_?NumericQ] := 
 NIntegrate[ep[u, s] + l1/(l2 - l1), {u, l1, l2}]

Example:

xn[3.2, -10]
xn[3.2, -100]
xn[3.2, -10^5]
(* -5.63342 *)
(* -5.63342 *)
(* -5.63342 *)

Edit As requested in comment, a new version with a stop criterion:

xnCrit[s_, maxIter_: 100, initVal_: - 0.1, step_: 3, tol_: 10^-3] := 
 Module[{i, xp, xq},
  i = 1;
  xp = xn[s, initVal];
  xq = xn[s, initVal*step];
  While[Abs[xp - xq] > tol && i <= maxIter, xp = xq; 
   xq = xn[s, initVal*step^(i + 1)]; i += 1];
  Print[initVal*step^(i + 1)];
  If[i == 101, Print["not converged for s = " <> ToString@s]];
  xq]

Usage example: xnCrit[3.2]. You can play with the parameters...