FREE VIBRATIONS

Remove["Global`*"]

Highlights

Abstract

Introduction

Off[General::"spell", General::"spell1"]; Needs["Graphics`Arrow`"]; tera = (10^12); micro = (1/(10^6)); giga = (10^9); nano = (1/(10^9)); mega = (10^6); pico = (1/(10^12)); kilo = (10^3); femto = (1/(10^15)); milli = (1/(10.^3)); atto = (1/(10^18)); NSub
       = {c -> 2.99792458 (10^8), e -> (1.6021773/(10^19)), 
        h -> (6.6260755/(10.^34)), hbar -> (1.05457/(10.^34)), 
        K -> (1.380658/(10^23)), me -> (9.109389/(10^31)), 
        mu0 -> ((4. π)/(10^7)), eps0 -> (8.85419/(10^12))}; ESi = 0.179 tera; rhoSi = 2330; nuSi = 0.28; ESiO2 = 0.06 tera; rhoSiO2 = 2200; nuSiO2 = 0.3; ESi3N4 = 0.15 tera; rhoSi3N4 = 3100; nuSi3N4 = 0.24; EAu = 0.08 tera; rhoAu = 19300; nuAu = 0.42; ENi = 0.207 tera; rhoNi = 8900; nuNi = 0.41; EW = 0.345 tera; rhoW = 19300; nuW = 0.28; EPt = 0.17 tera; rhoPt = 21400; nuPt = 0.39; EIr = 0.52 tera; rhoIr = 22500; nuIr = 0.3; EHOPG = 0.2 tera; rhoHOPG = 2250; nuHOPG = 0.3; EPol = 0.001 tera; rhoPol = 1000; nuPol = 0.3; xx = 0.8; EPtIr = xx 0.8 EPt + (1 - xx) EIr; rhoPtIr = xx rhoPt + (1 - xx) rhoIr; nuPtIr = 0.3; opt
       = {PlotRange -> All, AxesOrigin -> {0, 0}, 
        Frame -> True, GridLines -> Automatic, 
        DefaultFont -> {"Times-Roman", 8}, 
        PlotStyle 
          -> {Dashing[{.0001}], Dashing[{.0200}], Dashing[{.0200}], 
            Dashing[{.0200}], Dashing[{.0200}], Dashing[{.0001}]}};

Operating Conditions

case0
       = {abm -> 1, k -> 1, f0 -> 10, Q -> 5, 
        f -> 10}; case1 = {w -> 2 π f, w0 -> 2 π f0} /. case0; 
    caseL = case0 ⋃ case1;

case0
       = {abm -> 1, k -> 1, f0 -> 10, Q -> 50, 
        f -> 10}; case1 = {w -> 2 π f, w0 -> 2 π f0} /. case0; 
    caseH = case0 ⋃ case1;

Equation of Motion

Numerical Solution

Equation of Motion

Clear[z, v, x1, x2, eq1, eq2]; eq1=x1'[t]==x2[t];eq2=x2'[t]==-((w0 x2[t])/Q) - 
              (w0^2) (x1[t] - abm Sin[2 π f t])/.caseL;

Numerical Solution

n = 2 Q /. caseL; cycle = (1/f) /. caseL; stop = n cycle;

Clear[z, v]; {z, v} = {x1, x2} /. 
        NDSolve[{eq1, eq2, x1[0] == 0, x2[0] == 0}, {x1, x2}, {t, 0, stop}, 
            MaxSteps -> 100000][[1
          ]];

Transient Regime

AGraphics = 
      Plot[z[t], {t, 0, stop}, PlotRange -> All, 
        DisplayFunction -> Identity];

BGraphics = 
      Plot[abm Q Sin[2 π f t] /. caseL, {t, 0, stop}, 
        PlotStyle -> Dashing[{.02}], DisplayFunction -> Identity];

CGraphics = 
      Plot[{abm Q (1 - Exp[-((w0 t)/(Q 2))]) /. caseL, 
          -abm Q (1 - Exp[-((w0 t)/(Q 2))]) /. 
            caseL}, {t, 0, stop}, DisplayFunction -> Identity];

Fig7x2 = 
      Show[{AGraphics, BGraphics, CGraphics}, 
        DefaultFont -> {"Times-Roman", 8}, 
        FrameLabel -> {"t in sec", "A in a.u"}, 
        Frame -> True, DisplayFunction -> $DisplayFunction];

Clear[z1]; z1[t_] := abm Q Sin[2 π f t - (π/2)];

DGraphics = 
      Plot[z1[t] /. caseL, {t, 0, stop}, PlotStyle -> Dashing[{.02}], 
        DisplayFunction -> Identity];

Fig7x3 = 
      Show[{AGraphics, CGraphics, DGraphics}, 
        DefaultFont -> {"Times-Roman", 8}, Frame -> True, 
        PlotLabel -> "Figure 3", 
        FrameLabel -> {"t in sec", "A in a.u"}, 
        DisplayFunction -> $DisplayFunction];

Bimorph-Cantilever Phase Diagram

Fig7x4 = 
      ParametricPlot[{abm Q Sin[2 π f t] /. caseL, z[t]}, {t, 0, 
          stop}, AspectRatio -> 1, 
        FrameLabel -> {"a", "A"}, Evaluate[opt]];

Displacement-Velocity Phase Diagram

Fig7x5 = 
      ParametricPlot[{z[t], v[t]}, {t, 0, stop}, AspectRatio -> 1, 
        FrameLabel -> {"z", "v"}, Evaluate[opt]];

Analytical Solution

Equation of Motion

Clear[z, A]; z[t_] := A Exp[I w t]\(\(eq10 = 
      ∂\_{t, 2}z[t] + \(w0\ ∂\_t z[t]\)\/Q + 
          w0\^2\ z[t] == abm\ w0\^2\ Exp[I\ w\ t]; \)\)

Steady-State Regime

sol1 = Solve[eq10, {A}]; A = A /. sol1[[1]]; Print["A(w) = ", A]

A0 = ComplexExpand[Abs[A]]; \(Print["A0[w] = ", A0 /. \@x_\^2 -> x]\)

A0 = 
      ((abm Q (w0^2))/(Sqrt[(w^2) (w0^2) + 
            (((-Q (w^2) + Q (w0^2)))^2)]));

Bimorph-Cantilever Phase

\(B = ComplexExpand[A]; 

    \(phi = ArcTan[
        \(I\ B[[2]]\)\/\(B
              [[1]] + B[\([3]\)]\)]; 
    \)\)Print["phi[w] = ", Simplify[phi]]

Clear[fun, phi]; fun[x_] := ArcTan[(1/(Q ((1/x) - x)))] /. caseLphi[x_] := -If[x < 1, fun[x], π + fun[x]]

Fig7x6 = 
      Plot[((180 phi[x])/π), {x, .01, 2}, 
        AxesOrigin -> {1, -90}, 
        FrameLabel -> {"w over w0", "phi in deg"}, 
        Evaluate[opt]];

Q-Dependent Resonance Frequency

AA = Simplify[A (A /. w -> -w)]; w0Q = 
      Solve[(∂_w) AA == 0, w][[3, 1, 2
        ]]; Print["w0Q[Q] = ", w0Q]

Fig7x7 = 
      Plot[(w0Q/w0) /. 
          caseL[[{1, 2, 3, 4, 6, 7}]], {
          Q, 5, 15}, FrameLabel -> {"Q", "w0Q over w0"}, 
        Evaluate[opt]];

Frequency-Dependent Amplitude

fig7x8 = 
      Plot[Evaluate[A0 /. w -> x w0 /. caseL], {x, 
          1 - (1/Q) /. caseL, 1 + (1/Q) /. caseL}, 
        AxesOrigin -> {(w0Q/w0) /. caseL, 0}, 
        PlotRange -> {0, 1.2 Q /. caseL}, 
        FrameLabel -> {"w over wo", "A in a.u"}, Evaluate[opt]];

Frequency at the Steepest Slope

\(A01 = 
      A0 /. caseL[[{1, 2, 3, 4, 7}]]; 

    sol1 = Re[
        N[\(Solve[∂\_{w, 2}A01 == 0, w]\)[[4, 1, 2
              ]]] /. caseL]; 

    sol2 = Re[
        N[\(Solve[∂\_{w, 2}A01 == 0, w]\)[[2, 1, 2
              ]]] /. caseL]; \)

listL1
       = {"Exact", N[(sol1/w0) /. caseL], 
        N[(A0/(abm Q)) /. w -> sol1 /. caseL], 
        N[(sol2/w0) /. caseL], 
        N[(A0/(abm Q)) /. w -> sol2 /. caseL]};

wm1 = w0 (1 - (0.35/Q)) /. caseL; wm2 = w0 (1 + (0.35/Q)) /. caseL;

listL2
       = {"Approx.", N[(wm1/w0) /. caseL], 
        N[(A0/(abm Q)) /. w -> wm1 /. caseL], N[(wm2/w0) /. caseL], 
        N[(A0/(abm Q)) /. w -> wm2 /. caseL]};

listL
       = {"listL", "fm1/f0", "A1/A0", "fm2/f0", "A2/A0"}; 
    TableForm[{listL, listL1, listL2}]

A01 = 
      A0 /. caseL[[{1, 2, 3, 4, 7}]]; 
    \(\(sol1 = 
      Re[N[\(Solve[∂\_{w, 2}A01 == 0, w]\)[[4, 1, 
              2]]] /. caseH]; \)\)\(\(sol2 = 
      Re[N[\(Solve[∂\_{w, 2}A01 == 0, w]\)[[2, 1, 
              2]]] /. caseH]; \)\)

listH1
       = {"Exact", N[(sol1/w0) /. caseH], 
        N[(A0/(abm Q)) /. w -> sol1 /. caseH], 
        N[(sol2/w0) /. caseH], 
        N[(A0/(abm Q)) /. w -> sol2 /. caseH]};

wm1 = w0 (1 - (0.35/Q)) /. caseH; wm2 = w0 (1 + (0.35/Q)) /. caseH;

listH2
       = {"Approx.", N[(wm1/w0) /. caseH], 
        N[(A0/(abm Q)) /. w -> wm1 /. caseH], N[(wm2/w0) /. caseH], 
        N[(A0/(abm Q)) /. w -> wm2 /. caseH]};

listH
       = {"listH", "fm1/f0", "A1/A0", "fm2/f0", "A2/A0"}; 
    TableForm[{listH, listH1, listH2}]