Verilog-A large-signal model for a fixed-fixed beam RF MEMS resonator

/* Verilog-A large-signal model for a fixed-fixed beam RF MEMS resonator
Content of RF_MEMS_RESONATOR.va. Open-source code covered by BSD license.

The Verilog-A large signal model can be used for operating point, small signal (AC, SP)
and large signal (HB, PSS, QPSS) simulations, using ADS, SpectreRF or Qucs. It includes
large-signal electromechanical effects, as well as small-signal Brownian noise. It does
not include electro-thermal effects yet.

1st of October, 2011
Check for update: www-personal.umich.edu/~vcaeken

Reference: "Modeling RF MEMS Devices", by K. Van Caekenberghe,
IEEE Microwave Magazine, vol. 13, no. 1, pp. 83-110, Jan.-Feb. 2012 */

`include "disciplines.vams"
`include "constants.vams"
`define EPS 1e-10
// `define POW(x,y) pow( (x) , y )
`define POW(x,y) pow( (x+`EPS) , y )

module RF_MEMS_RESONATOR (in,out,beam);
    inout in, out, beam;
    electrical in, out, beam, dVindt, dVoutdt;
    kinematic z, velocity;

    // Model parameters
    parameter real m=2.275e-13 from (0:inf]; // beam mass [kg]
    parameter real b=2.721e-9 from (0:inf]; // damping coefficient [N.s/m]
    parameter real k_1=1.505e3 from (0:inf]; // k_i: stiffness parameters [N/m^i]
    parameter real k_2=0 from [-inf:inf];
    parameter real k_3=-6.2e15 from [-inf:inf];
    parameter real k_4=0 from [-inf:inf];
    parameter real k_5=1.2e30 from [-inf:inf];
    parameter real g_0=0.33e-6 from (0:inf]; // gap [m]
    parameter real C_0=0.185e-15 from (0:inf]; // capacitance [F]

    real C_m, E_max, eta, f_m, F_e, F_s, i_out_max, kappa, k_1_eff, k_3_eff, L_m, Q_m, R_m, Vin, Vout, v_in_max, z_c;

    analog begin

    @ ( initial_step ) begin

        $strobe("\n%M: b (Damping coefficient [N.s/m]) = %E",b);
        $strobe("%M: m (Effective beam mass [kg]) = %E",m);
        $strobe("%M: k_1 (Stiffness parameter [N/m^1]) = %E",k_1);
        $strobe("%M: k_2 (Stiffness parameter [N/m^2]) = %E",k_2);
        $strobe("%M: k_3 (Stiffness parameter [N/m^3]) = %E",k_3);
        $strobe("%M: k_4 (Stiffness parameter [N/m^4]) = %E",k_4);
        $strobe("%M: k_5 (Stiffness parameter [N/m^5]) = %E",k_5);
        $strobe("%M: C_0 (Transducer capacitance [F]) = %E",C_0);
        $strobe("%M: g_0 (Transducer gap [m]) = %E",g_0);

        k_1_eff=k_1-2.0*C_0*`POW(V(beam),2.0)/`POW(g_0,2.0);
        $strobe("%M: k_1_eff (Effective stiffness parameter [N/m^1]) = %E",k_1_eff);
        k_3_eff=k_3-4.0*C_0*`POW(V(beam),2.0)/`POW(g_0,4.0);
        $strobe("%M: k_3_eff (Effective stiffness parameter [N/m^3]) = %E",k_3_eff);
        f_m=sqrt(k_1_eff/m)/(2.0*`M_PI);
        $strobe("%M: f_m_e (Mechanical resonant frequency [Hz]) = %E",f_m);
        Q_m=sqrt(k_1_eff*m)/b;
        $strobe("%M: Q (Mechanical Q factor []) = %E",Q_m);
        eta=V(beam)*C_0/g_0;
        if (eta == 0) eta=`EPS;
        $strobe("%M: eta (Electromechanical coupling coefficient (DC-biased beam) [A.s/m]) = %E",eta);
        R_m=sqrt(k_1_eff*m)/(Q_m*`POW(eta,2.0));
        $strobe("%M: R_m (Motional resistance (DC-biased beam) [Ohm]) = %E",R_m);
        C_m=`POW(eta,2.0)/k_1_eff;
        $strobe("%M: C_m (Motional capacitance (DC-biased beam) [F]) = %E",C_m);
        L_m=m/`POW(eta,2.0);
        $strobe("%M: L_m (Motional inductance (DC-biased beam) [H]) = %E",L_m);
        kappa=3.0*k_3_eff/(8.0*k_1_eff)*sqrt(k_1_eff/m);            
        $strobe("%M: kappa (Amplitude-frequency coefficient (DC-biased beam) [1/(m^2.s)]) = %E",kappa);
        z_c=sqrt(4.0*sqrt(k_1_eff/m))/(sqrt(3.0*sqrt(3.0)*Q_m*(-kappa)));
        $strobe("%M: z_c (Critical displacement (DC-biased beam) [m]) = %E",z_c);
        i_out_max=eta*sqrt(k_1_eff/m)*z_c;
        $strobe("%M: i_out_max (Maximum output current (DC-biased beam) [A]) = %E",i_out_max);
        v_in_max=i_out_max*R_m;
        $strobe("%M: v_in_max (Maximum input voltage (DC-biased beam) [V]) = %E",v_in_max);
        E_max=0.5*k_1_eff*`POW(z_c,2.0);
        $strobe("%M: E_max (Maximum reactive energy (DC-biased beam) [J]) = %E",E_max);

    end
    
    Vin=V(in,beam);
    Vout=V(out,beam);
    $strobe("\n%M: V_in [V] = %E", Vin);
    $strobe("%M: V_out [V]= %E", Vout);

        // F_e: voltage-controlled electrostatic force [N]
    F_e=0.5*C_0*g_0/`POW(g_0-Pos(z),2.0)*`POW(V(in,beam),2.0)-0.5*C_0*g_0/`POW(g_0+Pos(z),2.0)*`POW(V(out,beam),2.0);
    $strobe("%M: F_e (Electrostatic force [N]) = %E", F_e);

        // F_s: spring force [N]
        F_s=k_1*Pos(z)+k_2*`POW(Pos(z),2.0)+k_3*`POW(Pos(z),3.0)+k_4*`POW(Pos(z),4.0)+k_5*`POW(Pos(z),5.0);
    $strobe("%M: F_s (Spring force [N]) = %E", F_s);

    // Nonlinear state-space description
    Pos(velocity):ddt(Pos(z))==Pos(velocity);
    Pos(z):ddt(Pos(velocity))==1/m*(-b*Pos(velocity)-F_s+F_e+white_noise(4.0*`P_K*$temperature*b, "BROWNIAN_NOISE"));

    `define SLEW   
    `ifdef SLEW   
        V(dVindt)<+ddt(slew(V(in,beam),1p,-1p));
        V(dVoutdt)<+ddt(slew(V(out,beam),1p,-1p));
    `else
        V(dVindt) <+ddt(V(in, beam));
        V(dVoutdt)<+ddt(V(out,beam));
    `endif

    I(out,beam)<+(C_0*g_0/(g_0+Pos(z))*V(dVoutdt)-C_0*g_0/`POW(g_0+Pos(z),2.0)*Pos(velocity)*V(out,beam));
    I(in,beam)<+(C_0*g_0/(g_0-Pos(z))*V(dVindt)+C_0*g_0/`POW(g_0-Pos(z),2.0)*Pos(velocity)*V(in,beam));
    I(beam)<+I(out,beam)+I(in,beam);

    end

endmodule