3
$\begingroup$

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$:

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ 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!). $\endgroup$ – anderstood Apr 18 '18 at 18:51
  • $\begingroup$ @anderstood Could I make a more stupid mistake? :) Thanks a lot $\endgroup$ – Foad Apr 18 '18 at 21:17
  • 1
    $\begingroup$ 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}]] $\endgroup$ – Mariusz Iwaniuk Apr 18 '18 at 22:04
  • $\begingroup$ @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)? $\endgroup$ – Foad Apr 19 '18 at 7:12
  • 1
    $\begingroup$ See edit with updated version. $\endgroup$ – anderstood Apr 19 '18 at 11:31
5
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Foad Apr 18 '18 at 21:31
  • 1
    $\begingroup$ @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. $\endgroup$ – anderstood Apr 18 '18 at 21:36
  • $\begingroup$ any chance you remember that answer? $\endgroup$ – Foad Apr 18 '18 at 21:41
  • 1
    $\begingroup$ I think that it was something similar to mine except the + l1 was not in the integral. @MariuszIwaniuk $\endgroup$ – anderstood Apr 18 '18 at 21:46
  • 1
    $\begingroup$ 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 $\endgroup$ – anderstood Apr 26 '18 at 12:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.