I am trying to reproduce numerically the results found in this paper: https://ieeexplore.ieee.org/document/1707778

You don't necessarily need to read it. Basically I boils down to solving a two equations for two unknowns, deflection and voltage: ζ and V

\begin{equation*} (k_{0}+k_{1})\zeta+k_{2}\zeta^{3}=\mu f_{e}(\zeta) \tag{16} \end{equation*}

\begin{align*} f_{e}(\zeta)((k_{0}+k_{1})&+3k_{2}\zeta^{2})\\ &- \frac{df_{e}(\zeta)}{d\zeta}((k_{0}+k_{1})\zeta+k_{2}\zeta^{3})=0. \tag{18} \end{align*}

By solving (18) for ζ, and by substituting the value of ζ into (16), we determine the deflection and voltage.

Where: \begin{equation*} \mu=\frac{6\ell^{4}\epsilon_{0}\epsilon_{r}V^{2}}{Eh^{3}g_{0}^{3}}. \end{equation*}

\begin{equation*} k_{0}=\int_{0}^{1}(\bar w^{\prime \prime}(x))^{2}dx\\ k_{1}=N_{0} \int_{0}^{1}(\bar{w}^{\prime}(x))^{2}dx\\ k_{2}=N_{1}\left(\int_{0}^{1}(\bar w^{\prime}(x))^{2}dx\right)^{2}\\ f_{e}(\zeta)=\frac{1}{\mu}\int_{0}^{1}F_{e}(\zeta \bar w(x))\bar w(x)dx.\\ \end{equation*}

and

\begin{equation*} N_{1}=6\frac{g_{0}^{2}}{h^{2}}\\ N_{0}=\frac{N_{0}\ell^{2}}{EI}=12\frac{\tilde{\sigma}\ell^{2}}{Eh^{2}} \end{equation*}

\begin{equation*} F_{e}=\frac{\mu}{(1-\hat w)^{2}}\left(\mathcal{F}+\eta\frac{\partial \mathcal{F}}{\partial\eta}\right) \end{equation*}

\begin{equation*} \mathcal{F}(\beta,\ \eta)=1-0.36\frac{\beta}{\eta}+0.85(\frac{\beta}{\eta})^{0.76}+2.5\frac{\beta}{\eta^{0.76}}. \end{equation*}

where: \begin{equation*} \beta=\frac{h}{b},\ \eta=\frac{h}{g}. \tag{6} \end{equation*}

I am trying to implement it into Mathematica, but it is not working.

eq16[ζ_]:= -μf[ζ]+ζ^3 k2+ζ (k0+k1)
eq18[ζ_]:= f[ζ]((k0+k1)+3 k2 ζ^2)-D[f[ζ],ζ]((k0+k1)ζ+k2 ζ^3)

ℱ[η_,β_]:= 1 - 0.36 β/η + 0.85 (β/η)^0.76 + 2.5 β/η^0.76
F[w_]:= μ/(1-w)^2(ℱ[η,β]+η*D[ℱ[η,β],η])

f[ζ_] := 1/μ Integrate[F[ζ*w[x]]*w[x],{x,0,1}]

μ = (6*L^4*eps0*epsR*V^2)/(e*h^3*g^3)
k0 = Integrate[D[w[x],{x,2}]^2,{x,0,1}]
k1 = n0*Integrate[D[w[x],{x,1}]^2,{x,0,1}]
k2 = n1*(Integrate[D[w[x],{x,1}]^2,{x,0,1}])^2

n0 = 12 (sig*L^2)/(e h^2)
n1 = 6 g^2/h^2

(*Numeric Parameter*)
L = 100*10^-6;
b = 1*10^-6;
h = 2*10^-6;
g = 4*10^-6;
e = 169*10^9;
eps0 = 8.85*10^-12
epsR = 1
sig = 100*10^6*(1-0.066)
β=h/b
η=h/g


(*Trial Function*)
w[x] := 16 x^2 (1 - x)^2;

(*Solving the System*)
ζSolve = ζ/.Solve[eq18[ζ]==0,ζ]
VSolve = V/.Solve[eq16[ζSolve]==0,V]

It just keeps running and running and running... I guess that I need to solve it somehow numerically. In fact in the paper they write:

We emphasize that, once the trial function w(x) has been chosen, (18) reduces to a nonlinear algebraic equation for ζ, where the derivative df(ζ)/dζ can be computed for any ζ by numerical integration. In general, (16) cannot be solved analytically, and standard root finding techniques such as the bisection algorithm may be applied

...but how do you do that in Mathematica ?

Any help would be very much appreciated !!

P.S. The result for V should be around: 945

  • 1
    You have constants that differ by dozens of orders of magnitude. Numerical results will likely be very bad. Please, work in natural units, where all parameters are of order 1. – AccidentalFourierTransform Dec 3 at 19:58
  • @AccidentalFourierTransform This is true, but the author of the paper already thought about it, and hence all equations are dimensionless, so it should not be a problem. – james Dec 3 at 20:08
  • In your definition of eq16, do you want a minus sign in front of the $\mu f$ term? – LouisB Dec 3 at 21:00
  • @LouisB Very true. Thanks for the catch ! – james Dec 3 at 21:21
  • @james parameters not defined η,β. what are they needed for? – Alex Trounev Dec 3 at 23:13
up vote 3 down vote accepted

I read the article, made all the corrections.I managed to find the correct model that reproduces the data from table 3 of Romesh C. Batra, ASME, Fellow, Maurizio Porfiri, and Davide Spinello. Electromechanical Model of Electrically Actuated Narrow Microbeams. JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 5, OCTOBER 2006, p. 1175-1189

    eq16[\[Zeta]_, 
  V_] := -\[Mu][V]*f[\[Zeta]] + \[Zeta]^3 k2 + \[Zeta] (k0 + k1)
eq18[\[Zeta]_] := 
 f[\[Zeta]] ((k0 + k1) + 3 k2 \[Zeta]^2) - 
  D[f[\[Zeta]], \[Zeta]] ((k0 + k1) \[Zeta] + k2 \[Zeta]^3)
\[Mu][V_] := (6*L^4*eps0*epsR*V^2)/(e*h^3*g0^3)
\[ScriptCapitalF][\[Eta]_, \[Beta]_] := 
 1 - 0.36 \[Beta]/\[Eta] + 0.85 (\[Beta]/\[Eta])^0.76 + 
  2.5 \[Beta]/\[Eta]^0.76
F[x_] := 1/(1 - x)^2

(*Trial Function*)
w[x_] := 16 x^2 (1 - x)^2;

(*Numeric Parameter*)
L = 100*10^-6;
b = 1*10^-6;
h = 2*10^-6;
g0 = 4*10^-6;
e = 169*10^9;
eps0 = 8.854187817*10^-12;
epsR = 1;
nu = .066;
sig0 = 100*10^6;
sig = (1 - nu)*sig0;
p = {n0 = 12 (sig*L^2)/(e h^2),
   n1 = 6 g0^2/h^2, k0 = NIntegrate[D[w[x], {x, 2}]^2, {x, 0, 1}],
   k1 = n0*NIntegrate[D[w[x], {x, 1}]^2, {x, 0, 1}],
   k2 = n1*(NIntegrate[D[w[x], {x, 1}]^2, {x, 0, 1}])^2};
B = (\[ScriptCapitalF][\[Eta], \[Beta]] + \[Eta]*
       D[\[ScriptCapitalF][\[Eta], \[Beta]], \[Eta]]) /. {\[Eta] -> 
      h/(g0*(1 - \[Zeta]*w[x])), \[Beta] -> h/b} // Simplify;
G = Interpolation[
  Table[{\[Zeta], 
    NIntegrate[
     B*F[\[Zeta]*w[x]]*w[x], {x, 0, 1}]}, {\[Zeta], -0.5, .75, .01}]]
f[\[Zeta]_] := G[\[Zeta]]

(*Solving the System*)
 \[Zeta]0 = FindRoot[eq18[\[Zeta]] == 0, {\[Zeta], .62}]

(*Out[]= {\[Zeta] -> 0.620671}*)

 FindRoot[eq16[\[Zeta], V] /. \[Zeta]0, {V, 1}]

(*Out[]= {V -> 945.44}*)
  • Thanks a lot for your answer. Great idea with the interpolation ! Unfortunatly, when I run it, it does the interpolation, but somehow FindRoot does not work. My output is: "FindRoot[eq18[[Zeta]] == 0, {[Zeta], 0.4}]". It does not seem to evaluate it. – james Dec 4 at 7:58
  • I see that you put some constants (n0,n1,k1,k2) into brackets. Why ? – james Dec 4 at 8:13
  • @james I put the parameters in brackets for easier typing while debugging. What is your version, operating system and machine? – Alex Trounev Dec 4 at 11:38
  • I have Windows 10 and Mathematica 11.3 – james Dec 5 at 5:57
  • I also checked the article and I come to the same conclusion as you regarding the typo. I will try to check the math in the paper to find the error... but this could take a while – james Dec 5 at 5:59

Udated numerical data:

Clear["Global`*"]

(*Trial Function*)
w[x] := 16 x^2 (1 - x)^2

eq16[ζ_] := -μ f[ζ] + ζ^3 k2 + ζ (k0 + k1)
eq18[ζ_] := f[ζ] ((k0 + k1) + 3 k2 ζ^2) - D[f[ζ], ζ] ((k0 + k1) ζ + k2 ζ^3)

ℱ[η_, β_] := 1 - 0.36 β/η + 0.85 (β/η)^0.76 + 2.5 β/η^0.76
F[ww_] := μ/(1 - ww)^2 (ℱ[η, β] + η*D[ℱ[η, β], η])

I used ww since we already have a w[x]. Now what took forever for me was the following integration for x from 0 to 1. Indefinite integration works much faster and then manually apply the limits.

fz[x_] = 1/μ Integrate[F[ζ*w[x]]*w[x], x] // Simplify
(*Long answer in RootSum objects*)

f[ζ_] := fz[1] - fz[0] // Simplify

μ = (6*L^4*eps0*epsR*V^2)/(e*h^3*g^3)
k0 = Integrate[D[w[x], {x, 2}]^2, {x, 0, 1}]
k1 = n0*Integrate[D[w[x], {x, 1}]^2, {x, 0, 1}]
k2 = n1*(Integrate[D[w[x], {x, 1}]^2, {x, 0, 1}])^2

n0 = 12 (sig*L^2)/(e h^2)
n1 = 6 g^2/h^2

(*Numeric Parameter*)
L = 100*10^-6;
b = 1*10^-6;
h = 2*10^-6;
g = 4*10^-6;
e = 169*10^9;
eps0 = 8.85*10^-12
epsR = 1
sig = 100*10^6*(1 - 0.066)
β = h/b
η = h/g

Find ζ

Plot[Evaluate[eq18[ζ] // N // Chop], {ζ, -.5, .7}, PlotRange -> All]

enter image description here

FindRoot[Re[eq18[ζ]] == 0, {ζ, .55}] // Chop
(*{ζ -> 0.531648}*)

ζ = ζ /. %

Solve[(eq16[ζ] // Chop) == 0, V]
(*{{V -> -795.934}, {V -> 795.934}}*)
  • Very nice answer ! Thanks a lot ! I updated the numeric constants since I did a mistake when copying from the paper. However, the result is still not the same... – james Dec 4 at 8:14
  • @james What should be the result? – Alex Trounev Dec 4 at 11:46
  • @AlexTrounev It should be 945. I also spotted a little copy mistake in ℱ[η_,β_] but it is now fixed in the question. – james Dec 4 at 13:22
  • @james then we have to check all the data more carefully. – Alex Trounev Dec 4 at 14:07
  • Hi ! So I have contacted the author of the paper. He told me that we have to adapt η because the gap decreases with increasing voltage. Basically on needs to replace η = h/g with η = h/(g-w[x]). I tried it with your code, but I don't get any result. – james Dec 6 at 10:59

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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