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
eq16, do you want a minus sign in front of the $\mu f$ term? $\endgroup$