Analytical Model for MEMS Electret Energy Harvester with Long-stroke Tip-sliding Electrodes

We develop an equivalent circuit model for a MEMS vibrational energy harvester that uses electrets or permanent electrical charges to generate electrostatic induction currents from mechanical vibrations. An electrode pair of periodically arranged comblike fingers is electrically biased by the built-in potential of an electret, and the distribution of electrostatically induced charges is altered by the relative mechanical motion of the electrodes. The electrostatic force as well as the induction charges are described as a function of the boundary condition and implemented into a multiphysics equivalent-circuit model by using the nonlinear current sources of the simulation software LTspice. As a practical solution for avoiding computational error, we have eliminated the use of polygonal approximations for the conditional analytical model and newly introduced a geometrical modulation function based on sigmoidal functions, by which the analytical model has become mathematically smooth and twice-differentiable with respect to the displacement. The short-circuit waveforms of vibrational energy harvesting are reproduced by simulation and are in good agreement with the experimental results.


Introduction
In parallel with microsensors and actuators, energy harvesters are increasing their significance as the third application pillar in the field of MEMS owing to their potential market in the coming Internet-of-Things (IoT) era or the so-called trillion sensors, (1,2) where a tremendous number of wireless sensors are to be deployed in the environment to establish a portal for cyber-physical systems (CPS). (3) Most wireless sensors today use batteries as an energy source, and hence the users are bothered by frequent charging or replacement. Energy harvesters are therefore expected to be a game changer in the sense that IoT wireless sensors would be free from such costly maintenance routines. (4) Amongst various sources of energy, including sunlight, room light, electromagnetic waves, and heat, we opt to use environmental vibrations (5) as a redundant energy source that has long been left unused. Unlike the other types of energy sources, mechanical vibration energy passes through the thick walls of packages, and therefore, the harvester device can be placed within a hermetic seal for high reliability. Moreover, the hermetic seal can also provide the harvester device with a vacuum package for negligible air damping, thereby allowing it to generate a maximum output at a high efficiency.
The vibrational energy harvester is a transducer that converts kinetic energy into electrical energy. Regardless of power generation principles, including electrostatic, (6)(7)(8) electromagnetic, (9) and piezoelectric (10) types, the harvester output current is proportional to the velocity of the harvester mass. (5) Most vibrational energy harvesters are, therefore, designed to oscillate at a large mechanical stroke, which is enhanced by the mechanical resonance. (11) Constructing an analytical simulation model for such a long-stroke mechanism requires an approximation technique to render a nonlinear behavior compatible with numerical computation. For an electrostatic type in particular, the harvester characteristics, such as the electrostatic constraint force and induction current, are sensitive to the geometrical boundary condition of the working electrodes. At the equation level, one can use polygonal expressions of the formula to analytically continue the behavior function beyond the conditional boundaries; (12) for instance, two mathematical expressions are used for the electrostatic capacitance between comb-finger electrodes, which are used conditionally depending on the state of the comb either in the engage or disengage position. However, such acts frequently invite a halt in the computation owing to the mathematical indifferentiability at the analytic joints.
A solution to this problem is to use periodic functions such as trigonometric functions to artificially perturb the electrode geometry. (13,14) One can also use a numerical computation to prepare a look-up table for the electrostatic capacitance, for instance, as a function of the displacement, and use it to calculate the value at an arbitrary position by interpolation. (15) In either case, the modeling degree of freedom is reduced, particularly when the optimal design parameters are sought, by the rapid repetition of trial-and-error simulation.
In light of computational compatibility, we have developed a new analytical model for the electrostatic vibrational energy harvester that uses a built-in electret film. (16) The model is based on the electrostatic analysis of capacitance, force, and charge induction on combfinger electrodes. Owing to the use of sigmoidal functions, the electromechanical behavior can be seamlessly connected beyond the conditional boundaries, and therefore, the numerical computation becomes free of faults even when a large-stroke behavior is investigated.
In this work, we look into the fundamental analytical model for such electret comb-finger electrodes, and derive the basic equations for the electrostatic force and induction current. A mathematical method used to seamlessly continue the model by using a sigmoidal function is discussed. The analytical model is then implemented in an LTspice circuit simulator (17) to calculate the energy harvester current and voltage as a function of the applied acceleration. The simulation results were found to be in good agreement with the experimental measurements of an actual harvester device.

Principle of electrostatic induction
Amongst various materials, we have chosen silicon oxide as the electret material (16) because of its compatibility with silicon micromachining technology. Silicon oxide has been known to keep an electrical charge for a long time after being electrically activated at an elevated temperature by the biased-temperature (BT) process. (18) The temperature of such a process can be lowered to a moderate level of 600 ℃ when an impurity-rich silicon oxide film is used as a starting material. In this work, we included potassium ions in silicon oxide during thermal oxidation, as reported elsewhere. (19) Figure 1(a) illustrates a simplified power generation mechanism based on the electrostatic induction current. A pair of electrodes is placed with a small gap between them. The electrode on the right-hand side is fixed to the mechanical anchor, while the other is elastically tethered with a suspension so that it can be used as a movable electrode. The surfaces of the electrodes are coated with silicon oxide, which, on the fixed electrode, is BT-processed to form an electret skin.
The silicon oxide after the BT process is negatively charged at a high sheet density of 1 mC/m 2 or more, owing to the dangling bonds of SiO − created after the removal of K + ions. Most electrical fluxes associated with these negative charges are terminated by the holes in the bulk silicon immediately next to the oxide film, as illustrated in Fig. 1(a). When the movable electrode is inserted by the external force, some fluxes are redirected to other holes in the movable electrode behind the air gap, and the overall distribution of electrons or electrostatic capacitances is altered by electrons rearranged by the new geometrical boundary conditions of the electrodes. Through this process, the mechanical work given to the movable electrode is converted into electrical energy by the principle of electrostatic induction, and electrical current flows through the external load resistance R. Owing to the nature of electrostatic induction, the magnitude of the current is proportional to the change rate of the boundary condition, i.e., the velocity of the movable electrode, and therefore electrostatic energy harvesters are usually designed to enhance the oscillation amplitude of the electrode through the mechanical resonance of a mass-spring system. However, the stroke of the electrodes is limited by the device footprint or geometrical design, particularly when the electrodes are formed into a comblike structure, as shown in Fig. 1(a).
To avoid such a geometrical conflict, we have employed an alternative electrode, as illustrated in Fig. 1(b); the structure still appears similarly to a typical comb electrode but it does not move in the direction of intercomb engagement, rather, it oscillates horizontally to slide the comb tips. The electrical capacitance becomes maximal when the tips are aligned, and it takes its minimal value when a tip is brought to the valley between the fixed comb tips. Therefore, the electrical capacitance periodically changes more than once per mechanical stroke. The stroke is limited not by the comb finger length but by the elasticity of the supporting suspensions, and therefore, the oscillation amplitude can be amplified to a scale greater than those of typical comb electrodes.

Experimental results
On the basis of the above new electrode design, we developed a tip-sliding vibrational energy harvester, as shown in the photograph in Fig. 2. The harvester chip is made in a siliconon-insulator (SOI) wafer processed using photolithography steps and a deep reactive ion etching (DRIE) process. Its area and thickness are 2 × 3 cm 2 and 0.5 mm, respectively. Two blocks of tungsten (total mass m = 12 grams) are attached to the oscillating part on the top and bottom surfaces to increase the mass, whereby more kinetic energy can be effectively received from environmental vibrations.
The structures including the folded-beam suspensions are schematically illustrated in Fig.  3(a); the fishbone structures in the figure are the movable electrodes; the fixed electrodes are not shown for ease of viewing. The detailed view of the electrode pair is shown in Fig. 3(b). The fixed electrodes are arranged on both sides of the movable electrodes to increase the output. Such an electrode block is placed symmetrically about the oscillation directions. Electrode blocks deliver electrostatic induction currents that are out-of-phase to each other, and therefore, their output ports are electrically separated. They are connected to a full-wave rectifier to enhance the conversion efficiency.
As a practical implementation of the electrodes shown in Fig. 1(b), we used such dimensions as a comb height of H = 300 µm, a comb length of D = 30 µm, a comb tip width of T = 8 µm, a comb pitch of P = 40 µm, and a tip-to-tip air gap of g = 9 µm. A total of 6182 pairs of comb fingers were used in a single device. The suspensions were designed to have a net elasticity axis aligned to the shaker stroke. The two output ports were short-circuited with a small electrical resistor (R = 1 Ω). The short-circuit current was observed with an oscilloscope after current-voltage conversion. The actual resonant frequency was observed at 44.9 Hz, despite the designed value of 51.3 Hz, possibly owing to inevitable manufacturing errors such as overetching by DRIE on slender suspensions.
When the device was shaken at the resonant frequency, the oscillation of the suspended mass was amplified by the quality factor (Q-factor) to a stroke larger than ±120 µm at a very small excitation acceleration of 1.5 m/s 2 , as shown in Fig. 4(a). The shaker used in this work was not capable of controlling such small acceleration waveforms owing to the lack of a feedback control system, and therefore, we tuned the shaker excitation magnitude to its minimum level and observed the acceleration by using a reference accelerometer attached to the shaker table. The excitation frequency of 44.9 Hz was close to the power supply frequency of 50 Hz, and therefore, a small acceleration signal was embedded in the measurement noise. After washing out the signal by subtracting the effect of the powerline, the acceleration signal became distorted, as shown in Fig. 4(b); nevertheless, the effective acceleration should be sinusoidal.
Judging from the waveform of the short-circuit current, with more than six positive peaks in a half-cycle, we speculate that the movable electrode oscillated at a stroke of ±120 µm, which corresponds to a three-unit length of the comb tips. Under this condition, the output current reached a peak of 3.5 µA.
When the device was shaken at 50 Hz, which is off the resonance frequency, the amplitude decreased to a very low level of no greater than 1 µA, even though large vibrations of 40 mm/s 2 or more were used; the short-circuit current under this condition was less than the 1 µA peak, as shown in Fig. 4(b).

Electrostatic force and induction current
The electrostatic induction current from the electrets, as well as the electrostatic force acting on the comb electrodes, could be analytically described by using the Gauss's theorem applied to the model shown in Fig. 5. The detailed mathematical procedure has been reported elsewhere, (20) and therefore, we simply follow the same protocol to derive the results. In Fig. 5, we presume that each comb pair can be described by a conductive block model with the tip width T, height H (in the direction of the SOI layer thickness), and air gap length g. The movable tip on the upper side has already been displaced in the x-direction by the amount L, and therefore, the overlap of electrodes becomes L + x when the movable electrode is further displaced by x. The comb electrode on the lower side has a silicon oxide layer of thickness s and is electrically charged into an electret density of σ C/m 2 ) at the depth r measured from the oxide surface. The upper electrode also has a silicon oxide layer in the actual device, but we ignored it for the simplicity of modeling; the effect of the silicon oxide on the top is included in the effective air-gap length g. The parameter ε 0 is the dielectric constant of vacuum, and ε 1 and ε 2 are relative dielectric constants of silicon oxide and air, respectively. In this model, the field concentration at the edges of the electrodes is ignored for simplicity, but the effect is included in the mathematical model explained in Sect. 3.2. Different from the analytical model developed in Ref. 7, where the electret layer is set on the surface of the dielectric material, we introduce a variable r to specify the depth of the electret layer; the model becomes equivalent to that in Ref. 7 when r = 0.
Applying Gauss's law to the electrical fields E 1 in the lower side of the silicon oxide, E 2 in the upper side, and E 3 in the air gap, we write These electrical fields build up a voltage difference equal to the applied dc voltage V, and thus, ( ) By simultaneously solving these three equations, one will obtain the electric fields as The electrostatic potential energy stored in the lower part of the silicon oxide is expressed in terms of the field E 1 as ( )( ) is the volume that contains the field E 1. Likewise, the energies stored in the upper part of the silicon oxide and the air gap are respectively written as ( ) and ( ) We also presume another set of electrostatic fields that are confined within the lower sides of the silicon oxide, where the surface does not overlap with the counter electrode. Such a set of electrostatic fields is separately calculated using the Gauss's law as ( ) 0 1 / σ ε ε , and accordingly the stored energy is written as where T L x − − is the length of the comb tip that extends out of the overlap. In summary, the total electrostatic energy stored in one unit of the comb tips is When considering the virtual displacement ∆x, the total energy will increase by ∆U, as defined by At the same time, the upper electrode will gain additional charges ∆Q U as where the factor H∆x is the increment of the area due to the virtual displacement. The negative sign is due to the definition of the direction of the electrical field E 3 , which is taken to be positive in the direction against the voltage drop on applying V. By considering the conservation of energy, the increment of the electrostatic energy and the mechanical work done by the electrostatic force F e are equalized by the work done by the external voltage V that transferred the charge ∆Q U against the potential difference as We use the expressions for E 1 , E 2, and E 3 in Eq. (4) and derived the electrostatic force F e as In this equation, the electret charge density σ has been replaced with the electret potential V e, which can be directly measured at the silicon oxide surface with respect to the lower electrode. The electret voltage V e is defined as Note that the magnitude of force is thought to be independent of the displacement x as long as the comb tips have a finite overlap, i.e., in the range L x T L − ≤ < − , where L is the initial insertion length. We have omitted the details of the process of formula deformation between Eqs. (12) and (13) because the intermediate expressions are written in a polynomial form with more than twenty terms, and therefore, we used a mathematic tool to symbolically process the formulation to obtain the final simplified expression as Eq. (13).
The electrostatically induced charges on the top and bottom electrodes are respectively written as and ( ) ( ) An equivalent model for the electret device can be implemented as a two-port electrical device that responds to the applied voltages as a capacitance. Provided that the upper and lower electrodes are individually biased to V U and V L , respectively, we can replace the voltage V with the differential voltage as An equivalent circuit is created as a two-port device with nonlinear current sources as The effect of electrostatic induction is thus implemented by the mechanical motion of the electrode, / dx dt, and also by the time variation of the applied voltages. These two functions determine the mechanoelectric coupling for the vibrational energy harvester.

Mathematically differentiable force curve
The electrostatic force model formulated in Eq. (13) can yield the magnitude of the electrostatic force only when there is an overlap between the comb tips, as schematically illustrated in Fig. 5. When the movable tip is off the opposing tip area, we assume that no countercharge is induced on the top electrode and that the electrostatic force falls to null.
For a long stroke displacement of electrodes, the electrostatic force can be presented by the polygonal curve shown in Fig. 6. In the figure insets, we schematically show the electrostatic capacitance as a tip-to-tip geometry. When the movable top electrode is off the left-hand side edge of the bottom electrode (x L < − ), we assume no electrostatic force. The force takes a constant value modelled by Eq. (13) when the top electrode has a finite overlap area with the bottom electrode ( L x T L − ≤ < − ); the direction of the force in this range is to the right, where the top electrode would be zipped to align with the bottom electrode. At a point x T L = − , the force is thought to be virtually null, and then it flips its direction when the movable electrode slides over to the right ( 2 T L x T L − ≤ < − ), but it still retains the same force strength as that in Eq. (13). Finally, the force becomes null again when 2 x T L ≥ − . In summary, the polygonal force curve is written with conditional branches as ( ) The polygonal expression would be useful for a thought experiment, but it is impractical for computer simulation because polygonal functions with conditional branches are not twicedifferentiable when the displacement occurs across the bending corners. In such a case, the simulation computation is usually aborted.
For this reason, we introduce a modest modulation function to synthesize smooth transitions between the conditional ranges. The standard sigmoidal function ( ) a x ς is written as where a is the gain constant and determines the transition slope between 0 and 1, as shown in Fig. 7(a); the slope becomes steep by decreasing the positive value of a. Note that a sigmoidal function has flat ranges that stick to 0 for x << 0 and 1 for x >> 0, and that the function is continuous and twice-differentiable. By combining the three sigmoidal functions shown in Fig.  7(b), the staircase-like function shown is reproduced as The synthesized function S(x) takes values of +1 and −1 for the flat top and bottom terraces, respectively, and it asymptotically approaches null when |x| >> T, as shown in Fig. 7(c). The function also smoothly connects the terraces through the origin. Therefore, this function can work as a modulator for the electrostatic force in Eq. (20), which can be now rewritten as Owing to the combined sigmoidal functions, the direction of the force is also included in the formula as a function of the position x. For the movable electrode that has been inserted at L as the initial position, the position x in Eq. (24) must be replaced with . x L + Unlike the original polygonal function in Eq. (20), the approximated curve has smooth corners as indicated in Fig. 7(c). The slope of the curve around the origin is Here, we define the transition range 2 1 x x x ∆ = − by the distance at which 2 S ∆ = from ( ) 1 1 S x = − to ( ) 2 1 S x = . Then, we write Considering the electrical equipotential surface that is conformally distributed around the electrodes, we allow a finite space for this transition ∆x that is as equivalently large as the half pitch P/2. As an empirical guideline, we write by which the coefficient of the sigmoidal function is related to the design parameter of the electrodes. The effect of this parameter will be verified in Sect. 5.2.

Differentiable electrostatic induction current
The electrostatic induction currents given by Eqs. (18) and (19) are rewritten as where W(x) is the function that shows the overlap length of the electrodes in the x-direction.
Considering that the initial position of the movable electrode is L and that the overlap length takes its maximal value of T when x = T − L, as shown in Fig. 8, the function W(x) is written in a polygonal form as For the same reason as in numerical computation, the polygonal formation for the overlap length W(x) should be remodeled by using another smooth curve. Noting that the polygonal function in Eq. (31) can be simply modelled by integrating the outline of the sigmoidal curve S(x) shown in Fig. 7, we synthesize a modulated function for the overlap length as Substituting this into Eqs. (29) and (30), we obtain a set of electrostatic induction current models, which is compatible with numerical computation. For the movable electrode that has been inserted at L as the initial position in Fig. 8, the position x in Eq. (33) should be replaced with x + L.

Model comparison
Besides the sigmoidal function, the Gaussian function is another option to represent the smooth and differentiable transition curves. Koga et al. (21) used the Gaussian function to fit the polygonal representation for the capacitance model as  min min where C max and C min are the maximal and minimal capacitances, respectively, as shown in Fig.  9(a). The curve is favorable for rounding off the bending corners of the polygonal function for the capacitance curve. However, its derivative form exceeds the polygonal expressions, and thus it tends to overestimate the electrostatic force, as shown in Fig. 9(b). The sigmoidal function shown by the red curve, on the other hand, remains close to the original polygonal expressions in the flat ranges.

Multiphysics implementation
The mathematical model for the electrostatic force and electrostatic induction current is implemented as an equivalent circuit that runs on the LTspice platform. We use the nonlinear current source of LTspice as a multiphysics behavior model that can be programmed by using the unitless physical variables in algebraic equations. The fundamental programming technique has been reported elsewhere. (22) In this work, therefore, we added a new equivalent circuit model for the tip-sliding electret mechanism as an electrostatic actuator, as shown in Fig. 10. From the left to the right in this schematic are the anchor module, dashpot, and spring model, and the mass model, which deals with the Newtonian equation of motion (EoM). The EoM module is programmed to read in the force as a current input and return the calculated displacement and velocity as voltage outputs.
The electret actuator module in this work is shown in the dashed red box in the figure and programmed to return the electrostatic force to the EoM module as a function of the displacement and applied voltage. It also returns the electrostatically induced charges to the voltage ports. A chain of submodules is attached to the right-hand side of the electrode module; they are modules of comb tips that are used to calculate the direction of the electrostatic force by considering the relative position of the movable electrodes with respect to the fixed ones. According to the experimental results shown in Fig. 4, the movable tip is thought to oscillate within a stroke range that covers 3 tips on the left and another 3 tips on the right, and hence, we put a total of seven sliding comb-electrode modules, as shown in Fig. 10 for the rest of the simulation in this work. The number of this module can be increased to adapt to a larger displacement. Figure 11(a) shows the circuit symbol of the electret actuator module. Physical constants, such as the dielectric constant of vacuum, are defined by the constant voltage sources within the module. The design parameters, such as the number of combs, the tip width, and the comb length, are provided from outside as constant voltages that are regarded as unitless; for instance, the comb length of 30 µm is converted to a constant voltage of 3 × 10 −6 V. Such constant parameters are terminated with a 1 Ω resistor to be read out as a voltage for internal use. The actuator module is programmed to hand over these design parameters to the comb modules that are attached to the right-hand side.
The electrostatic actuator module is programmed to calculate the electrostatic force and electrostatic induction current, as implemented in Fig. 11(b), by using the analytical equations written in the LTspice language. A hint in reading these equations is to understand that the argument parameters are passed as either current or voltage. Therefore, every parameter or variable has the prefix of V or I; for instance, the number of comb tips, N, is declared by a voltage source, and therefore, it is referred to as a voltage, V(N), in an LTspice equation.
The electrostatic force in Eq. (13) is converted into the LTspice format as  The LTspice model shown in Eq. (35) does not include the direction of the force. The direction is calculated in the comb-tip modules attached to the right-hand side. The result is read in by the electret actuator module and then multiplied with the magnitude of the force to deliver the resultant current as ( ) where V(direction) is a parameter provided from the comb-tip modules. Figure 12(a) shows the block diagram of the cascaded comb tips used to calculate the net effects of the force direction and the overlap length of the electrode. We take one of the movable comb tips and calculate the relative position with respect to each fixed comb tip; the relative position is then referenced to simulate the force modulation function S(x) by Eq. (23) and the overlap length W(x) by Eq.
(33). Such information is collected through the cascaded bus line, as schematically shown in Fig. 12(b), to deliver the net force modulation to the electret actuator module. In the case shown in Fig. 12(a), only the i-th fixed comb tip faces the movable tip with a finite overlap area, and therefore, the contribution to the electrostatic force and charge comes only from this module. The overall contribution from the entire comb array is then calculated by multiplying them with the number of comb fingers, N.
where the absolute position of the ith movable tip, x P i T L − ⋅ − + , is used as the first argument, the tip width T as the second, and g/4 as the constant to determine the slope of the sigmoidal function. The contribution from the tip module on the right-hand side is recursively added by the second term V(Di). The analog signals S(x) and W(x) transferred through the cascaded tip modules are finally handed over to the electret actuator module, in which the net force is calculated by Eq. (36). The total capacitor area W(x) is used to calculate the electrostatic induction charges appearing on the drive electrodes Vd and Vb using Eqs. (29) and (30), respectively.

Multiphysics implementation
We used the design parameters listed in Table 1 in the developed equivalent circuit model and calculated the short-circuit current, as shown in Fig. 13. The parameters marked with asterisks are used as fitting parameters that reflect the effect of the microfabrication errors. In particular, we set V e to be smaller than that used in the BT process, by considering the time-dependent degradation. The resonant frequency (51.3 Hz) was slightly higher than the experimental result of 44.9 Hz, because of the spring constant that had been made slightly smaller than the design as a result of the side etching of the micromachining process. Nevertheless, we compare the results without adjusting the frequency. Figure 13(a) shows the simulated short-circuit current of the energy harvester that is excited at the resonant frequency. The excitation amplitude was adjusted to 4 mm/s 2 to make the simulation results fit the results of the experimental measurement shown in Fig. 4(a), such that six positive peaks are seen in half a cycle. When the excitation frequency was changed to 60.3 Hz, the output current decreased to a small level, as shown in Fig. 13(b). This matches the experimental  observation in Fig. 4(b). In both cases, the excitation amplitudes must be slightly larger than those used in the experiment, possibly because of the elastic constant of the suspension model, which was larger than that of the developed device. On the whole, the equivalent circuit model reproduced well the experimental results. In Fig. 14, we show the growth sequence of the short-circuit currents at a fixed frequency with increasing excitation acceleration. When the acceleration was as small as 0.5 mm/s 2 , we could only see a sign of oscillation in the short-circuit peak current of 0.01 µA or less; the current oscillated in phase with the acceleration, as shown in Fig. 14(a). The peak current grew to 0.08 µA as the excitation was increased to 1.5 mm/s 2 , as shown in Fig. 14(b), until only one comb tip was covered within the stroke of oscillation. When the excitation was further increased to 2.0 mm/s 2 , two comb tips seemed to have been included within one stroke, since two positive peaks were seen in a half-period, as shown in Fig. 14(c). With increasing acceleration, more comb tips joined in the oscillation, as shown in Fig. 14(d), after which the results became close to the resonant waveform shown in Fig. 13(a).
Note that the amplitude shown in Figs. 13 and 14 is not the acceleration of the suspended mass (x ) but that of the applied excitation force. At resonance, the velocity of the spring-mass system becomes in phase with the excitation acceleration, and hence, the movable mass travels at its peak velocity at the top and bottom peaks of the excitation acceleration. Therefore, the short-circuit current takes its maximum values in phase with the excitation, as shown in Figs. 14(c) and 14(d). The peak current is out of phase with the acceleration peak in Fig. 14(b), possibly because the resonant frequency at this small amplitude is different from that at a large amplitude owing to the nonlinearity of the electrostatic force. By increasing the excitation acceleration, the oscillation amplitude of the electrode is increased to increase the pulselike output.

Verification
In this work, the electrical field concentration at the electrode edge was not considered in level of the analytical model, as shown in Fig. 5, but the effect was included in the mathematical model based on the sigmoidal function. As an empirical guideline, we have chosen the parameter a of the sigmoidal function, Eq. (21), as a = P/8, where P is the pitch of the comb tips.
In order to quantitatively verify the effect, we used numerical simulation tool COMSOL to develop a two-dimensional model for the comb electrodes and calculated the electrical potential distribution ( ) , U x y in the air gap, as shown in Fig. 15. By taking the spatial second-order differential ( ( ) 2 , U x y ∇ ), we calculated the electrostatic induction charges on the electrodes, and estimated the electrostatic capacitance C as a function of the electrode displacement x, as shown in Fig. 16. We used two gap lengths of g = 4 µm and g = 10 µm for comparison. We used Eq. (33) to calculate the effective overlap of comb tips, and superposed the characteristic curves in the same figure. When the sigmoidal parameter was taken to be a = P/8, the curves were found to explain the COMSOL numerical simulation. From this result, the use of a mathematical model based on the sigmoidal function is thought to be an effective method of modeling. Figure 17 shows the waveforms of the short-circuit current at resonance obtained by experiment (44.9 Hz) and simulation (51.3 Hz); the simulation waveform had to be stretched to fit the experimental result to compensate for the frequency difference. The simulation curve explains well the outline of the experimental results. However, it differs from the measurement in detailed structures, including the small terrace followed by the steep slope at the rising edge of each peak; the difference was caused possibly because the model with the sigmoid functions overestimates the gradient of the electrostatic potential distributed around the comb edges.
For a quantitative comparison, we used the integrated power of the short-circuit current as an index. Provided that the short-circuit resistance was 1 Ω, the waveform obtained by experiment for the first 12 ms gave a total energy of 1.21 × 10 −12 J. The total energy calculated from the  simulation curve in the same time span, on the other hand, was 1.36 × 10 −12 J. The equivalent circuit model made a 12% overestimate, which is acceptable as the first-order approximation of the device design.

Conclusions
We reported a new approach to developing an equivalent circuit model of a tip-sliding vibrational energy harvester with a built-in permanent electrical charge layer (electrets). Starting from Gauss's law, we derived the analytical models for the electrostatic force as well as the electrostatic induction current. A program technique was developed to implement the analytical model in the LTspice platform. To avoid computational error associated with the analytical continuation in polygonal expressions, we have newly introduced a smooth and twicedifferentiable modulator using the sigmoidal functions. The electret model was developed as an electrostatic actuator model and combined with the solver for the Newtonian equation of motion. The simulation results were found to explain the experimentally observed behavior of an electret-type vibrational energy harvester. The simulation technique presented in this work should be useful for mixed-signal simulations for the vibrational energy harvester with peripheral circuits such as rectifiers and storage circuits.