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

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.

• Your syntax is wrong in Clip (see documentation); try instead ep[u_, s_] := Clip[s/Ey[u], {0, 1}];. As a side note, you could easily find this from the error message Clip::rtwo and the context was not necessary (even if I find it very interesting in the present case!). – anderstood Apr 18 '18 at 18:51
• @anderstood Could I make a more stupid mistake? :) Thanks a lot – Foad Apr 18 '18 at 21:17
• At your wish: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_] := NIntegrate[ep[u, s], {u, l1, l2}, Method -> "LocalAdaptive"] + l1; ListLinePlot[Table[{s, xn[s, -100]}, {s, 0, 12, 1/10}]] – Mariusz Iwaniuk Apr 18 '18 at 22:04
• @MariuszIwaniuk Thanks a lot. @anderstood is there any way to tell Mathematica to calculate Integrate[ep[u, s], {u, l1, l2}] + l1 , iteratively varying l1 until it reaches a certain tolerance (e.g. 10^-5)? – Foad Apr 19 '18 at 7:12
• See edit with updated version. – anderstood Apr 19 '18 at 11:31

## 1 Answer

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...

• So basically I don't have to Limit l1 to minus infinity because at some finite negative point the epsilon is just zero because of the Mathematica computational limits. – Foad Apr 18 '18 at 21:31
• @Foad Yes (well, it's not really MMA computational limits but simply computation limits: computers don't have an infinite number of bits :P). But to make it clean you'd need a convergence test to find the right $l_1$ (such as: if I multiply $l_1$ by $2$, the relative change is smaller than $10^{-3}$). For example it might suffice to take $l_1=10^{-5}$ and that would then be unefficient to take $l_1=1$. Mariusz I... had posted a comment with another solution that provided the correct results; unfortunately he deleted it, I don't know why. – anderstood Apr 18 '18 at 21:36
• any chance you remember that answer? – Foad Apr 18 '18 at 21:41
• I think that it was something similar to mine except the + l1 was not in the integral. @MariuszIwaniuk – anderstood Apr 18 '18 at 21:46
• If it is bijective that should be doable, but you should ask a new question since Mathematica-wise, it is not related to the current question. @Foad – anderstood Apr 26 '18 at 12:11