440

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Magnetic Induction Micromachine—Part III: Eddy Currents and Nonlinear Effects Hur Koser, Member, IEEE, and Jeffrey H. Lang, Fellow, IEEE

Abstract—The magnetic induction micromachine fabricated in Part II (see this issue or Ref. [2]) was not laminated, as designed in Part I (see this issue or Ref. [1]). Consequently, eddy currents in the stator core, and the associated nonlinear saturation, significantly decreased its performance from that predicted in Part I. To investigate and explain these phenomena and their consequences, this paper models the behavior of the solid-stator-core machine fabricated in Part II using a finite-difference time-domain numerical analysis. The inherent stiffness in the time-domain integration of Maxwell’s equations is mitigated via reducing the speed of light artificially by five orders of magnitude, while taking special care that assumptions of magneto-quasi-static behavior are still met. The results from this model are in very good agreement with experimental data from the tethered magnetic induction micro motor. [1437] Index Terms—Eddy currents, finite-difference time-domain (FDTD), magnetic induction, micromotor, saturation. Fig. 1. A rendered visualization of the tethered induction micro machine structure (from [2]).

I. INTRODUCTION

T

HE tethered magnetic induction micromotor [2] is mainly a test bed for the microelectromechanical systems (MEMS) fabrication challenges. Torque measurements from this device were obtained using Kapton tethers and a separation layer that results in a rotor-stator air gap of about 70 m. The results, as reported in [2] and summarized here in Fig. 1, are about an order of magnitude less than the corresponding predictions of the linear design program described in [1] (see Fig. 2). The assumptions of that model include the linearity of the magnetic materials and the absence of eddy currents within the stator; eddy currents within the rotor are captured by the boundary layer formulation. Given these assumptions, the input current levels are chosen to keep the magnetic flux density within NiFe below its saturation value. This approach is valid as long as the eddy currents inside the stator are negligible, as is the case for a fully laminated stator structure. However, initial fabrication limitations did not permit a laminated stator; hence, the performance of the machine suffered. The resulting drop in torque, for example, can be seen by comparing Figs. 1 and 2. Again, this deterioration in performance originates from the presence of significant stator eddy currents, which Manuscript received September 28, 2004; revised November 28, 2005. This work was supported by the Army Research Office under Research Grant DAAG55-98-1-0292 and by the Defense Advanced Projects Research Agency under Research Grant DABT63-98-C-0004. Subject Editor O. Tabata. H. Koser is with the Electrical Engineering Department, Yale University, New Haven, CT 06520 USA (e-mail: [email protected]). J. H. Lang is with the Electrical Engineering and Computer Science Department, and the Laboratory for Electrical and Electronic Systems, Massachusetts Institute of Technology, Cambridge, MA 02139 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/JMEMS.2006.872240

require further modeling. The purpose of this paper is to model and explain the torque measured using the nonlaminated test machine. It should be stressed that the design process of a magnetic machine at the macroscale would not normally involve a detailed study of eddy currents; as such machines are finely laminated. However, the unusually high operating frequencies and the solid core stator of the micromachine necessitate such a study, which yields valuable insight that would help in the design of laminated structures in all high-frequency magnetic MEMS applications. In this paper, we develop a fully nonlinear magnetic model of the micro machine based on a finite-difference time-domain method. This modeling approach is quite convenient for the study of many electromagnetic phenomena and is especially suited for high-frequency operation regimes and nonlinear material applications. Based on this nonlinear model, we explain the data in Fig. 1. Finally, we conclude with a brief discussion of the next generation stators, which will be laminated to eliminate detrimental eddy current effects.

II. MAGNETIC DIFFUSION REVISITED In Part I [1], we present a coupled system of rotor and stator models for the magnetic induction micromachine. These models are then used to create a design tool to optimize the performance of the micromachine within prespecified fabrication limits. The modeling and the subsequent design optimization are performed

1057-7157/$20.00 © 2006 IEEE

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

441

Fig. 2. Measured torque from the tethered magnetic induction micromotor for several current levels. The device consists of a six pole stator and a rotor suspended above it via Kapton tethers. The error bars correspond to the maximum fit uncertainty.

under the assumption of a nonconductive NiFe stator. This assumption is a valid starting point, given that the eventual micromachines will be densely laminated [12]. In reality, though, NiFe has a relatively high conductance (about 5 10 S/m), and the frequencies of operation for the tethered micromotor yield skin depths that are smaller than the stator dimensions. This means that eddy current effects within a solid stator micromotor will be very significant and mostly responsible for the disparity between the predicted design torque and the lower measured torque. This hypothesis is supported by the results of a linear eddy current model that we have developed to study the effects of eddy currents inside the NiFe. A. Linear Eddy Current Model Maxwell’s equations for a magneto-quasi-static and stationary system reduce to magnetic diffusion equations [9]. In Part I [1], we have already employed the magnetic diffusion equations to obtain boundary layer relations for the rotor/conductor/air-gap section of the micro motor. Here, the same equations, in vector potential form, have been employed to the whole machine. Due to the complex structure of the stator geometry, these equations can be most easily solved using a finite difference or a finite element approach. Once again, we exploit the axial symmetry of the micromotor, and consider radial cross-sections of it as mapped onto a two-dimensional Cartesian plane. In steady-state, a sinusoidally varying current input to the wire slots within the stator results in a periodic vector potential whose fundamental harmonic has the same

frequency as the input excitation. As the wavelength of the traveling excitation spans four wire slots, the computation space may be reduced to include only a single wire slot, with appropriate boundary conditions (a complex phase difference of 2) on the left and the right, capturing a quarter of a wavelength difference. B. Flux Crowding The finite-difference modeling outlined above has been implemented for the particular geometry of the test magnetic induction micromachine. The magnetic permeability of electroplated NiFe is taken constant around 3000 , a medium-range value that is inferred from experimentally measured B-H curves of NiFe rings (see Fig. 8). The electrical conductivity of NiFe has been taken to be 5.0 10 S/m, and that of copper around 5.8 10 S/m. Fig. 3 illustrates the magnetic flux density within the micromotor just around a wire slot, computed with the linear eddy current model briefly outlined above. Here, the input current level is such that magnetic flux density inside the NiFe everywhere would be well below the saturation value (0.86 T) in the absence of eddy currents. When eddy currents are present (as in Fig. 3), however, the flux density inside the stator just around the wire slot increases with frequency well beyond the point where material nonlinearities would begin to take effect. This phenomenon (which we will henceforth refer to as “flux crowding”) appears to be a direct consequence of the traveling nature of the input excitation; it does not occur if the excitory waveform is simply a spatially uniform current varying sinusoidally in time.

442

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Fig. 3. Predictions of torque versus frequency from Model I for the same tethered magnetic micromotor tested. Comparison with Fig. 1 makes it clear that the torque values predicted by Model I are about an order of magnitude higher than the measured results; the peak frequencies are also higher than the measurements. The source of the discrepancy is the effect of eddy currents within the stator, which are not captured by Model I.

C. Need for a Nonlinear Model The simple linear model presented above suggests an explanation for the low values of measured torque. The high-frequency nature of the traveling excitation input leads to eddy currents within the highly conductive magnetic material, where the skin depth ends up smaller than the geometrical dimensions of the motor. This, in turn, leads to flux crowding within the rotor and especially the stator, squeezing the field lines toward the exterior NiFe surfaces. As the field lines are squeezed, the magnetic flux density within the NiFe surface around the wire slots rises dramatically above the saturation level. Of course, in reality, the inner surfaces will saturate, and the “saturation wave” will penetrate deeper, until the magnetic flux driven by the windings is cancelled by the combination of induced eddy currents and flux carried within the NiFe. Inside the tethered motor, the saturation wave travels with the excitation around the stator. As a result of eddy currents and saturation, a dramatic reduction in output torque is observed. It is clear that any model that has a hope of explaining the measured torque data of Fig. 1 must capture not only the presence of eddy currents but also the nonlinearity in the magnetization of NiFe. Such a model should also be able to include the effects of multiply excited modes, since the level of nonlinearity is likely to be intense enough to cause significant mixing of modes. Most commercial magnetics simulation packages, such as ANSOFT, solve a nonlinear eddy current problem by assuming that only the fundamental component of the vector potential

contributes to the solution. This is essentially equivalent to solving the magnetic diffusion equation with repeated Newton’s iterations over a user-defined B-H curve. One problem with this approach is that it does not capture the effects of higher harmonics, which are significant in the nonlinear micro machine. Also, mathematically speaking, using the linear magnetic diffusion equations in a nonlinear case is unwarranted. Therefore, it is not surprising that ANSOFT’s torque prediction for the same motor geometry and inputs as the tested tethered motor is far off in nonlinear runs, e.g., more than 300% higher than the measured value for the peak at 6 A (see Fig. 4). Our search for a more comprehensive, nonlinear, commercially available magnetics simulation package yielded no better candidates. It was decided that creating our own magnetics simulation software was the best solution. III. A FINITE-DIFFERENCE TIME-DOMAIN APPROACH A fully nonlinear magnetics solution is a demanding task computationally. A common way to tackle a nonlinear steadystate problem is to use the so-called shooting method [10]. The problem with this matrix approach is that all field variables and the material properties that change at every spatial point at every phase within a period must somehow be stored in computer memory. Even if memory is abundant, the method quickly gets overwhelmingly expensive in terms of the CPU time needed to solve the resulting system of equations. Improved algorithms, such as the generalized minimum residuals approach based on

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

443

Maxwell’s equations in the absence of free electrical charge and polarization (1) (2) (3) (4)

Fig. 4. Magnitude of B field (first harmonic) within the magnetic induction micro machine, as predicted by the linear eddy-current model. The scale units on the right are in tesla.

matrix-free Krylov-subspace methods, have recently been employed [5]; unfortunately, such simulations still take on the order of days to complete. In order to be practical, a nonlinear simulation tool for magnetic diffusion must be both accurate and very fast—on the order of less than an hour per simulation. One useful way to solve for magnetic diffusion in the microscale involves a direct, explicit time integration of the relevant partial differential equations over a finite-difference mesh. In this manner, matrix inversion and storage issues are all but eliminated. Given a clever setup of computation parameters, a finite-difference time-domain (FDTD) approach has the potential to provide the practicality that we seek in a nonlinear simulation method. The general FDTD method for electromagnetics, first outlined by Yee [6], solves the full set of Maxwell’s equations; hence, the approach could be used to study and explain virtually all electromagnetic phenomena [7], from photonic band gap structures to microwave circuits. However, applying FDTD algorithms to magnetic diffusion is not as straightforward as might first appear. In what follows, we first introduce our implementation of a FDTD magnetic diffusion algorithm. Later, we discuss some of the issues associated with FDTD in general and our implementation in particular. Finally, we present micromotor simulation results and discuss how they compare to experiments. A. Discretization Setup Once again, we will exploit the axial symmetry of the micro motor; the model will solve for electromagnetic fields in radial cross-sections as mapped onto a two-dimensional Cartesian plane. However, there are two major differences between our approach in this section and the linear eddy current model discussed above and in [1]. First, we need to use the full set of

This means that Ampere’s law includes the electrical induction term. Here, we choose a convention in which the total current density is represented as the sum of externally applied current density and the induced current; in those locations where a source sets the total current, the term is not present, and in those locations where induction sets the current, is not present. Second, the field variables in this time-domain approach take on real values. Therefore, complex phase cannot be used to implement symmetry. As a result, the smallest necessary computation space for the micro motor encompasses two wire slots, with odd symmetric boundary conditions enforced on either side, as illustrated in Fig. 5. 1) Electric Field Update: Yee’s algorithm calls for a discretization scheme in which magnetic and electric field points are separated both in space and time by half a discretization, as shown in Fig. 6. We shall evaluate electric field values at integer multiples of the time step, and find magnetic field values half a time step later, in a leapfrogging manner. Interestingly, it turns out that this method ensures that first-order finite differences for derivates are second-order accurate. Moreover, the finite difference grid, with each electric field point encircled by magnetic field lines, is divergence-free (in both and ) by construction [7]. The expression for the electric field update is obtained by integrating (1) over the surface of a given unit cell on the - plane at location and converting the right-hand side (RHS) into a line integral along around that cell, with field variables evaluated at time step (1/2). This yields

(5) Here,

we

shall

represent as and use a semi-implicit approximation to write and as and , respectively. Then, applying (5) to the two-dimensional discretization in Fig. 6, and solving for , we find

(6)

444

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Fig. 5. Comparison of ANSOFT’s simulated torque predictions at 6 A with measured torque values from the tethered magnetic induction motor. ANSOFT overpredicts the torque by more than 300%.

Fig. 6. The smallest computation space needed for the FDTD scheme, and the boundary conditions applied. For clarity, a diffuse mesh is shown.

where the function (7) has been defined for notational convenience. Assuming that magnetic field values at time step (1/2) have already been computed and that the input is specified for all time, (6) is explicit; it can be solved by direct substitution. 2) Magnetic Field Update: Another equation is required to find the magnetic field values at each spatial location; for that, we shall employ Faraday’s law (4) and follow the same approach as above. This time, however, in order to convert the

Fig. 7. The general meshing scheme used in an FDTD algorithm. Over each unit cell, corresponding field values and material properties are taken constant. Hence, material boundaries fall on the edges of the unit cells shown here.

curl of the electric field into a line integral, the contours must be around a surface pointing along the direction, as shown in Fig. 7. We begin by integrating both sides of (4) over the surface of Fig. 7 at time step ; converting the RHS into a line integral, we obtain (8)

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

445

Assuming that the surface of integration extends in the dimension, the RHS evalutes to . One must pay special care in evaluating the left-hand side of (8), since both the material properties and the flux density field values are not necessarily uniform over the surface of integration. Evaluation then yields

Fig. 8.

Surface for magnetic field update.

The magnetic flux density along the direction at unit cell can be found in exactly the same manner as above. The update equation for is

(14) (9) by using the continuity of We can solve for tangential material interfaces, which implies

across

(10) Notice that we have allowed the magnetic permeability of each unit cell to vary with time. Equation (8) finally reduces to

(11) The problem with implementing (11) is that we do not yet know at time step (1/2), as we would first need to find it. Hence, (11) is not fully explicit as desired. Short of solving a matrix inversion problem, this implicitness may be mitigated by using as the initial guess for and iterating (11) between the flux density and the permeability until eventually converges. Depending on the size of the time step, our experience is that a few such iterations generally suffice. is often so small that approximating In fact, the time step for the purposes of the above equation works quite well. Adopting this approximation, (11) simplifies to

3) Magnetic Permeability Update: Now that is known ev, the magnetic permeerywhere at time ability of the NiFe across the computational grid can be updated using the B-H curve of the magnetic material. For this purpose, a lookup table was created by fitting a function to the measured magnetization curve shown in Fig. 8. The fitting function used in Fig. 8 is of the form (15) where the parameters through are varied to obtain the best fit to the data. The variable is very close to 2.0. Notice that on the far right of Fig. 8 where full saturation is evident, (15) settles to a relative magnetic permeability of unity. The slope at the origin is also , as desired. Given and hence is found by direct interpolation. Care is taken to ensure that the magnetic field and the magnetic flux density vectors are collinear, i.e., (16)

4) Implementation Issues: A finite-difference scheme that is based on an explicit discretization of both space and time is subject to what is known as the Courant stability condition [7]. In the case of linear materials, a stability analysis of the FDTD equations reveals that (17)

(12) where (13)

is defined for notational simplicity.

is a necessary (but not sufficient) condition to guarantee stability of our numerical scheme. What is more, must be much shorter than the shortest wavelength or characteristic distance present in the problem. In the case of the micromotor, has to be chosen smaller than the shortest skin depth. This presents a problem. Assuming a skin depth of about 20 m in NiFe around 50 kHz, must be on the order of at most 10 m, resulting, according to (17), in a time step on the order of a few tens of femtoseconds. Since the minimum time span one needs to integrate

446

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Fig. 9. The measured magnetization curve of the NiFe around 7.5 kHz (stars) and the corresponding fit to data (thin solid curve). The slope of the curve at the origin and at high H values is  . Given a particular j j, the corresponding magnetic field magnitude is found using the fit (as illustrated by the thick solid line).

B

before a steady-state solution can be obtained is on the order of the excitation period (20 s), the simulation takes on the order of 10 time steps. With a unit cell count on the order of 10 and each unit cell requiring (10) operations for field computations, the number of total double precision operations needed is larger than 10 . It would hence take an average PC several years to conclude a single simulation, which is completely impractical. The tight constraint on stability arises from the inherent stiffness in the partial differential equations (PDEs): electromagnetic waves travel essentially at the speed of light (1 ) in air, yet the magnetic diffusion waves travel much slower inside NiFe. For instance, even at full saturation that drives , the diffusion wave speed is over 200 000 times slower than , the speed of light. Moreover, the input excitation has a phase speed of , which is also much slower than . Indeed, given the stator winding pattern, the input excitation completes a period over four wire slots, which, at a midradius of 1.5 mm at 50 kHz, corresponds to a phase speed of just under 1000 m/s, 300 000 times slower than light. If this inherent disparity between the different characteristic times could be mitigated, the Courant constraint could be relaxed to give a larger integration time step. In this modeling work, we adopt an implementation in which the speed of light is artificially reduced many orders of magnitude, in order to reduce the stiffness of the PDEs [11]. Since the emphasis is on magnetic diffusion phenomena, the reduction of involves increasing the electric permittivity of all materials in our computation space, without altering their magnetic permeabilities. Care is taken to make sure that is still many times

faster than both the input traveling wave and the magnetic diffusion waves within all materials. This approach has enabled a stable time step as large as 2 ns—an improvement of over 60 000 times. With this method, a simulation takes less than an hour to conclude on an average PC. To give a perspective, the same simulation would have taken about seven years without the modification. Another numerical technique that helps in achieving stability is to avoid sharp gradients by turning on the input excitation slowly and letting the system settle into its steady-state behavior gradually. In our particular implementation, we modulate the input current by a factor of (1 ), with chosen typically around a few microseconds . This allows enough time for initial magnetic fields to diffuse in and start to saturate the material, before the input reaches its maximum amplitude, which enables a larger stable time step while minimizing the overshoot and the consequent low-damped oscillations that would otherwise dominate the simulated output behavior of the micromotor. Depending on the input frequency, this numerical technique speeds up the time integration further by an order of magnitude. The FDTD simulation approach described here is also capable of computing eddy current and resistive losses within materials (see Section V-B). It is interesting to note in this context, however, that eddy current losses are lower in the presence of material nonlinearities, thanks again to reduced flux levels within the NiFe. Hysteresis losses, on the other hand, do not come into play until many tens of megahertz, and are therefore neglected in the loss computations.

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

447

Fig. 10. Some key frames of the magnitude of magnetic flux density within the tethered magnetic micromotor. Excitation period is 20 microseconds. Saturation in NiFe around the wire slots is clearly evident in these simulation results.

IV. MODELING THE MICROMOTOR The nonlinear FDTD equations have been solved within the computation space of Fig. 5. The stator and rotor geometry and

material properties have been mostly selected to match those of the actual tethered magnetic micromotor tested; see Fig. 11 and Section V below. The snapshots in Fig. 9 illustrate the evo-

448

Fig. 11.

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Time evolution of torque for a sample current input (2 A, 30 kHz). The dashed line corresponds to the steady-state torque value extracted.

lution of magnetic flux density within the micro motor, over a time span of about 1.5 periods at 50 kHz. Though not separately labeled, the wire slots and the stator-rotor air gap are clearly distinguishable due to high field boundaries within the surrounding NiFe surfaces. Recall that there exist two winding phases within each wire slot, carrying currents in quadrature. Since the skin depth within copper at the operating frequencies of the micromotor is larger than the smallest dimension of the wires, an effective uniform current density within the wire slot is assumed. In Fig. 9, the excitation amplitude is 6 A within each winding, corresponding to one of the experimental conditions. The heavy saturation regions around the wire gaps, and especially the bottleneck for flux at the stator teeth, indicate that most of the stator pole and backiron volumes do not contribute to flux linkage at all. The same can be said for most of the rotor backiron, as well. Already, one can see that intense saturation is the main reason for low values of torque measured with the tethered micromotor. The second and the third columns of pictures in Fig. 9 are organized such that they are almost half a period apart in time within each row. By the time one full period has completed (third column), we can observe an interesting phenomenon: a new saturation pattern around a wire slot begins to ripple out before an earlier saturation region has time to completely dissipate.

Fig. 12. A radial cross-section schematic of the magnetic induction micromotor and the relevant dimensions used to simulate it.

Several simulations at different radii are performed, and the total torque is approximated through a trapezoidal sum, which yields the instantaneous total torque acting on the rotor. The time-evolution of torque is then plotted, as in Fig. 10, and the steady-state value is extracted from the mean torque once the simulation settles.

A. Computing Torque Torque output from the micromotor is calculated in the same manner as outlined in [1]. The Maxwell stress tensor is integrated around the rotor section (from on the left to on the right) in the computation space of Fig. 5 to find the shear force acting on the rotor. The computed shear force is multiplied by the radius and number of poles to obtain the instantaneous torque at that radius (18)

V. RESULTS AND DISCUSSION A. Certain Uncertainties Recreating the exact experimental conditions in numerical simulations is not possible. To begin, our model is a two-dimensional approximation of the three-dimensional micromachine. Certain three-dimensional effects such as eddy currents turning corners at the inner and outer radii are not captured in our model. Another issue is that local nonuniformities in each material’s properties—such as voids within the electroplated NiFe—are not considered in these simulations. Moreover, parameters such

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

449

Fig. 13. Torque predictions from simulations for various tooth gaps. Measured values are overlaid for reference. The rotor conductor thickness is taken to be 8 m, and the temperature within the micromotor is constant at 250 C above room temperature.

as dimensions and operating temperature vary across the micromotor, and using their average values is yet another approximation that introduces potential errors. In this section, we discuss those parameters whose variation or uncertainty has the most prominent effect on simulated results. 1) Geometric Dimensions: Fig. 11 depicts the relevant dimensions of the tethered micro motor that was tested. Notice that some of the dimensions are either approximate or they span a range of values. The mask set was designed for a 10 m tooth gap, for instance, but homemade photo-reduction chromium masks end up having a substantial edge roughness, and aspect ratio limitations of SU-8 end up increasing the gap. Consequently, the tooth gap, on the average, is much closer to 15 m than 10 m. The rotor conductor, on the other hand, was fabricated to be around 8 m, and the real average thickness is probably very close to that value, though the layer exhibits thickness variations due to edge effects. Because the most prominent saturation effects take place at the maximum current, the measured torque at 6 A provides the ideal set of data with which to study the effects of geometric variations. Figs. 12 and 13 illustrate the simulated effects of variation in the tooth gap and rotor conductor thickness, respectively, while all other parameters are kept constant. A tooth gap of 15 m, in combination with a rotor conductor thickness of 8 m, does indeed yield the best agreement with the measured torque. It is interesting to note here that, in the absence of eddycurrents and saturation, the uncertainty in the tooth gap does not change the solution. However, in the nonlinear case, output torque is quite sensitive to stator tooth gap. Simulating the mi-

cromotor with an average tooth gap is an optimistic approximation, as it ignores localized changes in saturation and flux density and flux rerouting, all of which eventually effect the output torque. Hence, simulation results should be evaluated with this inevitable uncertainty in mind. Finally, consider the rotor-stator air gap, which is depicted to be around 70 m in Fig. 11. In fact, the Kapton die that houses the rotor core sits on top of a spacer layer that is 65 m thick, but the rotor core itself is slightly recessed within its housing. This is probably because the rotor was manually dropped into its housing over a glass slide and glued from top; as the glue dried, it shrank and pulled the rotor up with it. The recess is nonuniform around the circumference of the rotor, with a mean value of about 5 m. The resulting average rotor-stator air gap is about 70 m, which is the dimension used in simulations. Once again, using this average value inevitably introduces systematic errors in the nonlinear simulations. 2) Magnetization Curve Characteristics: Besides some of the geometric dimensions, another variable that introduces uncertainty is the magnetization curve. Specifically, the B-H characteristic of electroplated NiFe has been measured using thin rings of this material wound in a transformer configuration, which permits an ac measurement. Due to voltage resolution limitations of the experimental apparatus, the lowest measured frequency is around a few kilohertz. Ideally, however, a B-H curve at dc is used to determine saturation effects; this way, eddy current effects are excluded from the measurement. The magnetization curve in Fig. 8, however, was measured at 7.5 kHz, a frequency at which the skin depth within the NiFe rings is beginning to get comparable to their thickness. How, then, are

450

Fig. 14.

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Torque predictions from simulations for various rotor conductor thicknesses. Measured values are overlaid for reference.

we justified in using that curve as our input B-H relationship for the FDTD simulations? The answer is that the micromotor operates so deep inside the saturation region that it does not matter too much whether we use the magnetization curve at a midfrequency or at dc. In order to characterize the sensitivity of the numerical setup to uncertainties in the magnetization characteristics of NiFe, we simulated the micromotor with the magnetization curve of Fig. 8 and the same curve shifted lower by a factor of two. In the worst case (see Fig. 14), the in difference in simulated steady-state torque has been found to be 7%, with the typical discrepancy being around 4%. In fact, simulated steady-state torque results are quite insensitive to shifts in the B-H curve of electroplated NiFe, as long is known accurately and is much smaller than the as typical field magnitudes within NiFe. This insensitivity could be applied as a numerical technique to simulate the NiFe wafer, which saturates at least two orders of magnitude earlier in than the electroplated NiFe. Recall that one of our assumptions in arriving from the field update [see (11)] to the final, explicit (12) and (14) is that the magnetic permeability does not change by much over a time step. With that assumption, we have (19) for any and . However, the approximation in (19) becomes invalid for the B-H characteristics of the NiFe wafer, because it saturates too fast, and the magnetic permeability actually changes appreciably during a unit time step. The result is numerical instability in the FDTD algorithm. This problem could be addressed by reducing the time step substantially to make (19) valid. Alternatively, shifting the B-H curve instead to relieve the stiffness in the problem accommodates a much

larger time step; see Fig. 15. The output torque result based on the shifted B-H curve yields identical results as that based on the original B-H curve, without the numerical instability. Once again, we can get away with such a lateral shift in the B-H curve because the input current levels for the micromotor are of the large enough to guarantee operation well over the resulting curves. 3) Temperature: Since the rotor and the stator surfaces face each other across a very small air gap, incorporating thermocouples on those surfaces during torque measurements is not possible. Even adding a thermocouple on the rotor backiron is not practical, since we must know the moving rotor mass accurately to determine torque [2]. Using infrared cameras to deduce operating temperatures is also not possible because the rotor core blocks the region of interest from view. Hence, temperature remains a variational parameter in our simulation studies. The relative insensitivity of the results to B-H curves shifted along the H axis comes in handy in the consideration of uncertain temperatures. To first order, the effect of temperature on the [8]. As long NiFe magnetization characteristics is to shift as the operating temperature is not very near the Curie temperaremains mostly constant. ture (450 C for 80–20% NiFe), As will be discussed below, the maximum temperature at the stator and rotor surfaces is estimated to be below 300 C, and the resulting shift in the B-H curve [8] does not significantly alter the simulation results. The effect of temperature on the B-H curves is therefore neglected. A change in temperature, however, changes the electrical conductivity of both NiFe and copper, and this is where temperature effects become critical. An increase in temperature reduces the NiFe conductivity and increases the skin depth within

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

451

Fig. 15. The time evolution of simulated torque output from the magnetic micro motor, with two different B-H curves for the electroplated NiFe. The slightly lower steady-state value of torque corresponds to the B-H curve, which is shifted toward lower H by a factor of two from the original curve.

Fig. 16. The evolution of torque with two different B-H curves for the NiFe wafer. The two transients are essentially identical, except that the original B-H curve eventually results in numerical instabilities.

the material, which, in turn, reduces saturation and tends to increase the torque. Such a temperature increase also reduces the rotor copper conductivity, which has the effect of shifting the peak of the torque toward higher frequencies. This is because an increase in the rotor conductor resistance reduces the time constant of the rotor. Also, the peak torque values are lowered as the magnitude of the induced currents inside the rotor copper drops off slightly with increased resistivity. Fig. 16 illustrates the overall simulated effect of temperature on the micro motor performance. Here, the rotor and the stator are assumed

isothermal; given the small dimensions, the intimate contact of the two dies, and the long operating times at each current level, this is a reasonable assumption. An overall temperature of 250 C appears to provide the best match to the measured data; the corresponding simulation data also predicts the peak frequency correctly. An exhaustive study of temperature around this “nominal” value has indicated that a temperature variation of 50 C about 250 C results in slightly worse, yet plausible, fits to measurements; see Fig. 17. It is reasonable for this temperature not to be significantly higher than

452

Fig. 17.

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Simulated torque values at different overall temperatures within the micromotor.

the glass transition temperature of SU-8 (around 200 C), as the stator of the micromotor that yielded the data in Fig. 1 lost functionality beyond 6 A of input current amplitude.

wafer) and the rotor conductor by summing up resistive loss densities in these volumes. This yields

B. Temperature at Lower Current Levels

is the depth of the computation space (1 mm), and where the extent of the summation is chosen to include the region of interest. Incorporating (20) into the FDTD simulation program, we have computed estimated power losses within the test device. Fig. 18 illustrates how the power dissipation due to stator winding and contact resistance compares to the rest of the dissipation processes. We have chosen to depict the comparison given the worst case eddy current loss, which occurs at the highest applied frequency given a current amplitude (winding losses are independent of frequency). Notice that even in the worst case scenario, the winding losses are still an order of magnitude larger than the rest of the contributions to dissipated power.1 In a way, this is welcome news, as it implies that in an eventually integrated magnetic micromachine where contact resistance will not be a major issue, effective heat sinking techniques could achieve a very small temperature rise within the device. Together with our assumption of an isothermal micro machine, Fig. 18 suggests that the operating heat loss of the test device depends mainly on the square of the input current amplitude. Assuming thermal coefficients independent of temperature, we find that the micromachine temperature rise is proportional to input current amplitude squared. In other words, in

We could make similar temperature fits to the measurements at lower current levels. However, consistency once again requires that we have a temperature model that applies to all current inputs. The factors that have an impact on heating within the losses within the stator windmicro machine are resistive ings, resistive losses within the rotor copper, and eddy current losses within the NiFe core; hysteresis losses are inconsequential in this frequency range. The eddy current losses depend both on frequency and input current amplitude. Resistive losses in the stator include the resistance of the copper windings at either phase, as well as the contact resistance of the current probes. The total resistance of the stator current path has been measured to be around 0.5 ; given the conductivity of copper (5.8 10 S/m) and the cross-sectional area of the windings (200 m 65 m inside the wire slots and inner end-turns, wider at the outside), (a good estimate is about 0.05 ) of that is less than 0.1 winding resistance. Hence, stator resistive losses are dominated by the contact resistance. Interestingly, it turns out that the overall dissipation within the tested tethered micromachine is also dominated by contact resistance. Since FDTD simulations compute, among other variables, the electric field everywhere inside the computation space, it is straightforward to extract iinstantaneous eddy current losses within the NiFe core (both electroplated and the stator

(20)

1Interestingly, Fig. 18 also shows that, unlike winding losses, eddy current losses do not exactly scale as I in the presence of saturation.

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

Fig. 18.

453

Simulated torque values at temperatures between 200 and 300 C above room temperature make good fits to the measured data at 6 A.

order to estimate the operating temperatures for current levels lower than 6 A, we need to scale the temperature differences from ambient in Fig. 17 with the ratio of the current amplitudes squared. Fig. 19 illustrates the results. The worst case discrepancy between the measurements and simulation results is about 30%, in the case of 3 A current amplitude. This is a significant improvement over commercial finite-element analysis (FEA) program results (Fig. 4). We believe part of the discrepancy at 3 A may be due to poorly balanced stator phases during the experiment at this current amplitude, which would reduce the measured torque. Also, the simulations are two-dimensional and cannot model such effects as eddy-currents turning at the inner and outer radii of the micromachine. Such three-dimensional effects act to reduce the measured torque, and they are more prominent at lower current levels where saturation is less severe. These effects could account for part of the discrepancy, especially at lower current amplitudes. VI. CONCLUSION Compared with the results from commerial FEA software packages such as ANSOFT, the FDTD results and their agreement with the measured data in a consistent manner provide evidence for the value and the validity of the general modeling approach presented here. Not only does the modified FDTD method improve simulation accuracies by an order of magnitude, it is also quite fast. The modified FDTD approach is a promising candidate for a microscale nonlinear electromagnetic simulation package; as implemented, the method can be used to study any magnetic actuator and microsystem. It can also be

modified easily for nonlinear microsystems in the electrical domain, such as electrostatically actuated micromirrors for optical applications. The extension of the FDTD algorithm into three dimensions is also straightforward, though it would be significantly more demanding in terms of memory and CPU time. The power of the FDTD approach is not limitless, however. In linear cases, established matrix methods can usually provide results much faster. The FDTD method is better suited to study transient effects and nonlinear dynamics associated with stiff partial differential equations that may otherwise be impractical to solve using other methods. The results of this paper indicate that stator eddy currents must be eliminated in order to achieve the high performance numbers predicted in [1]. In the case of the magnetic induction micromachine, the NiFe material properties, together with the operating conditions of the device, determine the lamination thickness necessary to achieve close to full power. Estimating the necessary lamination thickness is straightforward; it involves solving a linear magnetic diffusion equation of the form (21) within a thin NiFe slab (see Fig. 20), subject to the boundary condition of essentially equal field magnitudes at the edges. The general field solution within the magnetic slab is then given by a normalized function along the thickness of the slab. If the micro gas turbine rotates at 2.4 Mrpm, the corresponding synchronous electrical frequency is 120 kHz for a six-pole stator. Assuming a 10% slip in the generator mode,

454

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

Fig. 19. Comparison of simulated power loss mechanisms within the test device. Power dissipation due to stator winding and contact resistance dominates the overall heat generation process.

the operating frequency at the peak torque will be 108 kHz, resulting in a necessary lamination thickness of approximately 17 m to recover 90% of the torque in the absence of eddy-currents. This estimation does not take the packing factor of laminations into account; when factored in, a practically achievable packing factor of 50% will reduce the output torque accordingly. Early versions of the laminated devices are intended to test fabrication challenges and limitations; for simplicity, they will not incorporate laminations at either the stator or rotor backiron cores. Depending on the lamination thickness that can be achieved, we believe that most of the design torque may be recovered. In order to achieve full torque, however, the entire machine must be laminated.

It is our conviction that research for better magnetic materials must be conducted in parallel with lamination efforts. The eddy currents within the NiFe can be brought down by reducing the electrical conductivity of the electroplated NiFe in the first place. In macroscale manufacturing, a small percentage of silicon compounds is introduced into the steel to reduce its conductivity. In the microscale fabrication, deliberate introduction of low-density nonmagnetic and nonmetal ionic contaminants that dissolve in water, such as calcium ions, into the electroplating bath may help reduce the conductivity of the NiFe substantially. Such impurities will no doubt degrade the magnetic properties (such as decrease the magnetic permeability) of the NiFe slightly, but the accompanying reduction in electrical

KOSER AND LANG: MAGNETIC INDUCTION MICROMACHINE—PART III

Fig. 20.

455

Simulated torque results for 2–5 A. Temperature values are scaled down from the 6 A fits by the square of the current amplitude ratios.

conductivity is well worth it. Besides, a lower magnetic permeability actually helps increase the skin depth within the magis not netic material. As long as the saturation flux density substantially lowered, the resulting NiFe will significantly improve the performance of the magnetic induction micromachine, making mildly laminated stators successful. ACKNOWLEDGMENT The authors would like to thank J. White and S. D. Senturia for insightful discussions and F. Donovan and all other friends who donated valuable computing time that made this paper possible. REFERENCES Fig. 21. Potential lamination schemes for the stator (Figure courtesy in part of Florent Cros).

[1] H. Koser and J. H. Lang, “Magnetic induction micro machine—Part I: Design and analysis,” J. Microelectromech. Syst., vol. 15, no. 2, pp. 415–426, Apr. 2006.

456

JOURNAL OF MICROELECTROMECHANICAL SYSTEMS, VOL. 15, NO. 2, APRIL 2006

[2] F. Cros, H. Koser, M. G. Allen, and J. H. Lang, “Magnetic induction micro machine—Part II: Fabrication and testing,” J. Microelectromech. Syst., vol. 15, no. 2, pp. 427–439, Apr. 2006. [3] W. Hemmert, M. S. Mermelstein, and D. M. Freeman, “Nanometer resolution of three-dimensional motions using video interference microscopy,” in Proc. IEEE Int. Conf. MEMS, Orlando, FL, Jan. 1999, pp. 302–308. [4] S. D. Senturia, Microsystem Design. Norwell, MA: Kluwer Academic, 2001. [5] S. Li and H. Hofmann, “Numerically efficient steady-state finite element analysis of magnetically saturated electromechanical devices using a shooting-Newton/GMRES approach,” in Proc. IEEE Int. Electric Machines Drives Conf., Piscataway, NJ, 2001, pp. 275–279. [6] K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag., vol. AP-14, pp. 302–307, 1966. [7] A. Taflove, Computational Electrodynamics. Norwood, MA: Artech House, 1995. [8] R. M. Bozorth, Ferromagnetism. Piscataway, NJ: IEEE Press, 1993. [9] H. H. Woodson and J. R. Melchner, Electromechanical Dynamics. New York: Wiley, 1968. [10] A. Granas, R. B. Guenther, and J. W. Lee, “The shooting method for the numerical solution of a class of nonlinear boundary value problems,” SIAM J. Numer. Anal., vol. 16, no. 5, pp. 828–836, Oct. 1979. [11] R. Holland, “FDTD analysis of nonlinear magnetic diffusion by reduced c,” IEEE Trans. Antennas Propog., vol. 43, no. 7, pp. 653–659, 1995. [12] J. Park, F. Cros, and M. G. Allen, “A sacrificial layer approach to highly laminatedmagnetic cores,” in Proc. IEEE Int. Conf. MEMS, Las Vegas, NV, Jan. 2002, pp. 380–383.

Hur Koser (M’03) received the B.S. degree in physics degree in 1999, and the B.S., M.Eng., and Ph.D. degrees electrical engineering and computer science from the Massachusetts Institute of Technology, Cambridge, in 1998, 1999, and 2002, respectively. Since 2003, he has been an Assistant Professor of Electrical Engineering at Yale University, New Haven, CT. His other research interests are in the design and development of magnetic micromachines and other micropower devices, microsensors and actuators, and microfluidic devices with an emphasis on biomedical applications. Dr. Koser has received the NSF Career Award to work on the hydrodynamics of ferrofluids in 2005.

Jeffrey H. Lang (S’78–M’79–SM’95–F’98) received the S.B., S.M., and Ph.D. degrees in electrical engineering from the Massachusetts Institute of Technology (MIT), Cambridge, in 1975, 1977, and 1980, respectively. Currently, he is a Professor of Electrical Engineering at MIT. He has been an MIT faculty member since receiving his Ph.D. degree and his research and teaching interests focus on the analysis, design and control of electromechanical systems with an emphasis on rotating machinery, microsensors and actuators, and flexible structures. He has written over 160 papers and holds 10 patents in the areas of electromechanics, power electronics, and applied control. Dr. Lang has been awarded four Best Paper prizes from various IEEE societies. He is a former Hertz Foundation Fellow and a former Associate Editor of Sensors and Actuators.