

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


MODELING AND CONTROL OF A PERMANENTLY MAGNETIZED SPHERE’S MOTION BY FIELD MANIPULATION By MATTHEW GLENN DUVALL Bachelor of Science Mechanical Engineering Oklahoma State University Stillwater, Oklahoma, USA 1994 Master of Science Mechanical Engineering Oklahoma State University Stillwater, Oklahoma, USA 1997 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy December, 2006 COPYRIGHT c By MATTHEW GLENN DUVALL December, 2006 MODELING AND CONTROL OF A PERMANENTLY MAGNETIZED SPHERE’S MOTION BY FIELD MANIPULATION Dissertation Approved: Eduardo A. Misawa Dissertation Advisor Prabhakar R. Pagilla Frank W. Chambers Martin Hagan Dean of the Graduate College iii ACKNOWLEDGMENTS The author would like to sincerely thank the following individuals: Firstly, my academic advisor Dr. Eduardo A. Misawa whose support and advice greatly facilitated this project. Secondly, my academic committee, specifically Dr. Prabhakar R. Pagilla. His guidance, mentorship, and not to mention the many hours of time sacrificed to examining draft manuscripts proved to be invaluable to my educational experience. I would also like to thank Dr. Gary Young for his support during my teaching assistantship which provided a crucial aspect of the doctoral candidacy experience. In addition I would like to thank Dr. Kenneth L. Pottebaum of Seagate Technology, LLC for his professional guidance and counseling regarding the balance of school and academia. Moreover, I would like to thank Dr. Ryan T. Ratliff of the Boeing Company for his assistance in, and clarification of, various technical difficulties. Furthermore, I would like to thank Leon E. Boomershine andWenjie Li. Without assistance from these two individuals it is doubtful that the experimental hardware setup for this project would have been completed in a timely fashion. Finally, I must offer an extreme amount of gratitude to my wife, Heather D. Duvall, for her unflinching support in my academic endeavors and the many sacrifices for which she has had to endure. Without the support of these individuals, it is inconceivable that this body of work would be completed and this particular scholastic goal accomplished. iv TABLE OF CONTENTS Chapter Page 1 Introductory Work 1 1.1 Impetus for the Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Overview of Relevant Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Contributions of this Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Overview of this Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Dynamic System Model 8 2.1 Dynamic Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Inverse Current Solution for Constrained 2D Motion with a Reduced Coil Set . . . 14 2.3 Current Trajectory Smoothing Technique . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 State Space Model and Closed Loop Control Considerations . . . . . . . . . . . . . . 26 3 Electromagnetic Field Description 32 3.1 The Single Filamentary Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Multiturn Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Quad Coil Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Coil Resistance and Inductance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 Design and Optimization of Experiment 48 4.1 Functional for Minimization and Algorithm . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Minimization Results and Experimental Setup . . . . . . . . . . . . . . . . . . . . . 55 4.3 Experimental Platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Experimental Validation of Platform Design . . . . . . . . . . . . . . . . . . . . . . . 59 5 Closed Loop Control of Sphere Motion 65 5.1 Coil Switching Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2 Closed Loop Control Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 Closed Loop Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 v 6 Conclusions and Future Work 89 6.1 Conclusions from this Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.2 Recommended Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 BIBLIOGRAPHY 95 A Partial Spatial Derivatives of g Field Components Resulting From a Single Rectangular Filamentary Loop 99 B Power Circuit Diagram 103 C Closed Loop Matlab/Simulink Simulation Source Code 104 C.1 Model Parameters Script Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 C.1.1 Coil Resistance and Weight Calculation Script . . . . . . . . . . . . . . . . . 108 C.1.2 Coil Self Inductance Calculation MEX Script . . . . . . . . . . . . . . . . . . 109 C.1.3 Coil Mutual Inductance MEX Script . . . . . . . . . . . . . . . . . . . . . . . 116 C.1.4 G Field Calculation MEX Script . . . . . . . . . . . . . . . . . . . . . . . . . 121 C.2 Simulink Top level Model and Subsystems . . . . . . . . . . . . . . . . . . . . . . . . 134 C.3 Embedded Matlab Function Block Scripts . . . . . . . . . . . . . . . . . . . . . . . . 147 C.3.1 Square path embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . 147 C.3.2 Point path embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . . 148 C.3.3 Inverse current calculation embedded function. . . . . . . . . . . . . . . . . . 149 C.3.4 Current selector embedded function. . . . . . . . . . . . . . . . . . . . . . . . 209 C.3.5 K matrices embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . . 212 C.3.6 Wall forces embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . . 217 vi LIST OF FIGURES Figure Page 2.1 Typical M − H curve for a ferromagnetic material. . . . . . . . . . . . . . . . . . . . 10 2.2 Schematic of a four coil arrangement to provide forces for planar motion of a sphere. The individual coil coordinates indicate direction of positive current polarity (right hand rule). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Graphical representation of the applied planar magnetic force. . . . . . . . . . . . . . 15 2.4 Solution loci for the 1&4 coil combination in the J space. . . . . . . . . . . . . . . . 20 2.5 Initial current solution for coil combination 1 and 3 during circular trajectory. . . . . 22 2.6 Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces. . . . . . . . . . . . . . . . . . . . . . 22 2.7 Correct current solution for coil combination 1 and 3 during circular trajectory. . . . 23 2.8 Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces are now absent. . . . . . . . . . . . . 24 2.9 A plot of sgn(D1,3), where D = det[eQ ]. A value of +1 indicates that transformation is a rotation, while a value of −1 is indicative of a reflection. . . . . . . . . . . . . . . 24 2.10 Smoothed current solution for coil combination 1 and 3 during circular trajectory. . 26 2.11 Initial and smoothed current solution for coil combination 2 and 3 during circular trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.12 Initial and smoothed current solution for coil combination 2 and 4 during circular trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.13 Initial and smoothed current solution for coil combination 1 and 4 during circular trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1 Schematic of a single filamentary conducting loop. The local coordinate system is fixed to the center of the loop with the central axis aligned with the z coordinate direction. Positive current flow obeys the right hand rule. The numbers indicate sections of conductor while the letters denote the dimensions of the rectangle. . . . . 33 vii 3.2 Schematic of a rectangular coil with a “perfect wind” packing pattern. The diameter of the wire is denoted as “d”, while the number of turns/layer and number of layers are denoted t and l, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.1 Calculated current derivative from inverse solution of coil combination 1&3. Displayed results correspond to system with R1,3 = 0.547 , L1,1 = 6.79mH, and L1,3 = 0.58mH (calculated). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Instantaneous power required for inverse current solution of coil combination 1&3. . 51 4.3 Plot of k 1,3k for solution sets a and c. The norm of solution sets b and d are identical to sets a and c, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 A flowchart illustrating the algorithm used to find the number of turns tf and number of layers lf which results in minimal power used for a given ain and bin. . . . . . . . 53 4.5 A flowchart illustrating the algorithm used to find the inner x coil dimension ain,f and inner y coil dimension bin,f which results in minimal power used for a given t, l, and . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.6 A flowchart illustrating the complete algorithm used to find the parameters (tf , lf , ain,f , bin,f ) resulting in minimal power. The nested algorithms are indicated by “broken line” boxes. 55 4.7 Experimental apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.8 Schematic of experimental setup and signal flow. . . . . . . . . . . . . . . . . . . . . 57 4.9 Schematic of tray cross section. The sphere resides below the centerline of the fluid gap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.10 Experimental apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.11 Experimental setup for measurement of the flux density produced by a single coil in the y direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.12 Experimental setup for measurement of the flux density produced by a single coil in the z direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.13 A comparison of the predictive model for the field component Gy to measured data points. The brackets around the points indicate ±3 distributions of the measurement technique at that point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.14 A comparison of the predictive model for the field component Gz to measured data points. The brackets around the points indicate ±3 distributions of the measurement technique at that point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.1 Calculated current trajectories for different coil combinations. . . . . . . . . . . . . . 67 viii 5.2 Block diagram of closed loop control simulation with a classical “proportional” controller is implemented. The upper diagram is the closed loop system while the lower diagram is a detailed view of the controller. . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 Block diagram of closed loop control simulation with a feedforward block. A classical “proportional” controller is implemented in concert with a feedforward block which estimates the desired forces on the spherical particle. The upper diagram is the closed loop system while the lower diagram is a detailed view of the controller. . . . . . . . 71 5.4 Simulation of closed loop PD controller resultant trajectory following a spiral reference. 72 5.5 Simulation of closed loop PD controller resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.6 Simulation of closed loop feedforward controller resultant trajectory following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.7 Simulation of closed loop feedforward controller resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.8 Experimental closed loop controller resultant trajectory following a spiral reference. . 75 5.9 Experimental closed loop controller resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.10 Experimental closed loop controller with feedforward block resultant trajectory following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.11 Experimental closed loop controller with feed forward block resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . 77 5.12 Experimental closed loop controller resultant trajectory following a square reference. 77 5.13 Experimental closed loop controller resultant current and position error following a square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.14 Experimental closed loop controller with feedforward resultant trajectory following a square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.15 Experimental closed loop controller with feedforward resultant current and position error following a square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.16 Experimental closed loop controller resultant trajectory attempting a point regulation. 80 5.17 Experimental closed loop controller resultant current and position error attempting a point regulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.18 Experimental closed loop controller resultant trajectory attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 ix 5.19 Experimental closed loop controller resultant current and position error attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.20 Experimental closed loop feedforward controller resultant trajectory attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.21 Experimental closed loop feedforward controller resultant current and position error attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . 83 5.22 Experimental closed loop feedforward controller resultant trajectory attempting to track the fast spiral trajectory. (kP = 0.5) . . . . . . . . . . . . . . . . . . . . . . . . 84 5.23 Experimental closed loop feedforward controller resultant current and position error attempting to track the fast spiral trajectory. (kP = 0.5) . . . . . . . . . . . . . . . . 84 5.24 Experimental closed loop controller resultant trajectory following the fast square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.25 Experimental closed loop controller resultant current and position error following the fast square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.26 Experimental closed loop feedforward controller resultant trajectory following the fast square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.27 Experimental closed loop feedforward controller resultant current and position error following the square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.28 Experimental open loop controller resultant trajectory attempting to follow a spiral. 87 5.29 Experimental open loop controller resultant current and position error attempting to follow a spiral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 B.1 Circuit diagram for coil power op amps and signal conditional op amps. Note that these are inverting amplifiers and that each coil has a power circuit. . . . . . . . . . 103 C.1 Top level Simulink Model (“closed loop”). . . . . . . . . . . . . . . . . . . . . . . . . 134 C.2 Spiral trajectory subsystem (closed loop!spiral trajectory). . . . . . . . . . . . . . . 135 C.3 Spiral trajectory path definition (closed loop!spiral trajectory!spiral). . . . . . . . 136 C.4 Square trajectory subsystem (closed loop!square trajectory). . . . . . . . . . . . . . 137 C.5 Point regulator subsystem (closed loop!point regulator). . . . . . . . . . . . . . . . 138 C.6 Controller subsystem (closed loop!Controller). . . . . . . . . . . . . . . . . . . . . . 139 C.7 Inverse current calculation (closed loop!Controller!inverse current calculation). . . 140 C.8 Inverse voltage conversion (closed loop!Controller!inverse voltage conversion). . . 141 C.9 Required force subsystem (closed loop!Controller!required force). . . . . . . . . . 142 x C.10 Plant model subsystem (closed loop!plant model). . . . . . . . . . . . . . . . . . . 143 C.11 Electrical dynamics subsystem (closed loop!plant model!electrical dynamics). . . 144 C.12 Magnetic force subsystem (closed loop!plant model!magnetic force). . . . . . . . . 145 C.13 Sphere dynamics subsystem (closed loop!plant model!sphere dynamics). . . . . . 146 xi LIST OF TABLES Table Page 2.1 Correct Line Loci Groupings for e J Space Solution Coordinates . . . . . . . . . . . . 19 2.2 Modified Sign Groupings for Smooth e J1,3 = ˜j 1 1 1 ,˜j 3 3 3 Solution Coordinates . . 25 2.3 Modified Sign Groupings for Smooth e J2,3 = ˜j 2 2 2 ,˜j 3 3 3 Solution Coordinates . . 28 2.4 Modified Sign Groupings for Smooth e J2,4 = ˜j 2 2 2 ,˜j 4 4 4 Solution Coordinates . . 29 2.5 Modified Sign Groupings for Smooth e J1,4 = ˜j 1 1 1 ,˜j 4 4 4 Solution Coordinates . . 29 4.1 System Parameters to be Used for Optimization . . . . . . . . . . . . . . . . . . . . 49 4.2 Optimized System Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3 Targeted wound coil parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 xii CHAPTER 1 Introductory Work 1.1 Impetus for the Research In recent decades many new technologies have been proposed and investigated in the biomedical sciences. Several studies have suggested the use of magnetic nanoparticles in vivo to the human body. The benefits of such methods are believed to be invaluable because of the minimally invasive capabilities of anticipated procedures. Envisioned applications of nanoparticles include, but are not limited to, cell labeling, magnetic separation, targeted drug delivery, hyperthermic cell treatment, and MRI contrast enhancement [1]. Significant study has been performed regarding these concepts and is available in references [2, 3, 4, 5]. A second novel application of a magnetically reactive device being introduced to the internal workings of a living body and being manipulated by an applied ex vivo magnetic field is magnetic implant and catheter guidance. The practice of this procedure on the human brain is referred to as stereotactic neurosurgery and a significant effort has been put forth on this technique [6, 7]. Similarly, a process utilizing nano sized particles and applied external magnetic fields for intraocular retinal repair is being explored as well [8]. Finally, the fascinating and emerging science of nanoscaled machines has offered the potential employment of microdevices in medical procedures. To reduce device complexity, an external magnetic propulsion technique is envisioned [9, 10]. The common attribute present in all of the aforementioned systems is the use of an applied magnetic field. This field is utilized to exert a desired force and/or torque on an object susceptible to such fields. Upon that realization, one may attempt to first recognize the analogy and similarities to the classic magnetic levitation problem. Indeed a parallel may be observed that manipulation of a particle or microdevice along a unidirectional trajectory is comparable to the wellknown technique of suspending or manipulating objects by either the application of alternating fields or utilizing a direct field to maintain the static position of the target object [11, 12]. However, in the aforementioned biomedical processes, the ability to either sustain a static position in a moving fluid medium possessing varying components of velocity in multiple coordinate directions, or propel 1 the target object(s) along a predetermined spatial trajectory elevates this problem to one beyond classical levitation. Achieving this type of dynamics requires application of forces with components in multiple coordinate directions. The resulting fields, field gradients, and the means by which to generate them are far more complex than that of the levitation problem, which often requires force generation only along a single direction. A system such as this presents an interesting control problem for which a suitable descriptive model is necessary before a sufficient treatment and analysis may occur. The technique of moving an implant or microdevice by fields from static electromagnets is not a new concept [6, 7, 13]. Indeed a substantial amount of work has been performed using a cubic arrangement of coils to facilitate three dimensional motion. A recent technique as described in [13] utilizes six active coils to impart a predetermined force and torque on a target object. The studied approach powers all coils simultaneously with the three opposite the primary attractive coils being operated with a lower current of reverse polarity. This would result in an increase in the field gradients at the location of the target and in effect raise the magnetic force on the object. While this technique has been explored, calculation of the required currents is quite tedious. The device operates in an open loop control scheme and employs minimization techniques to numerically solve for the six coil currents. The difficulty in the current solution for this system rests with the fact that it is underdetermined (i.e., there are more equations than unknowns). For the scenario of moving a spherical object where the orientation is unconstrained, a solution of six unknowns (the currents) from three equations is needed. Hence, many solutions exist, requiring the use of optimization methods to single out a minimal current set. This work explores the technique of reducing the number of coils to only those which are necessary to generate the force. This results in a determined system, allowing for an exact calculation of the current set. The analysis will be limited to the less complicated problem of moving a sphere two dimensionally in a fluid medium using static electromagnets. The system consists of four electromagnets, of which only two are used at any time to provide planar motion to the particle. The pair of active coils which will be utilized is determined by the force requirement on the sphere. Before exploring the dynamics of this system in detail, a synopsis of the work available in the literature will be presented. 1.2 Overview of Relevant Research The use of magnetic particles in a biomedical application is not new and the available work is quite extensive. Early concepts consisted of using magnetic fluids to perform magnetic tracing and 2 mapping of blood flow systems [14]. This technique is “passive” in nature relying on the high permeability characteristic of the introduced fluid for detection by an electromagnetic pickup array. A second concept explored by [14] involves actively guiding and controlling an unmixed and undiluted material in the bloodstream. The envisioned technique would be helpful in sealing off defects in the circulatory system allowing for the initiation of clotting and scarring. Work on this has been limited with very little success due to the nature of ferrofluidic material not being easily controlled thru magnetic propulsion methods. More recent proposals involve manipulation of particles rather than a fluid based material. The envisioned techniques by which an active guidance system utilizing magnetic propulsion may be employed consist of targeted drug delivery and guidance of micro devices in blood vessels. For the technique of targeted drug delivery, the vast majority of present methods utilize uniform magnetic fields in concert with fluid flow dynamics. In this method, the magnetic particles are moved to the desired location by way of fluid system and applied uniform magnetic fields are then applied to hold the particles in the desired location [1, 2, 5]. While these techniques are promising, a more robust solution would be the employment of an active system to guide the particle to the desired location, independent of the flow dynamics. While this would be the ideal system, present technology is somewhat limited due to the fact that extremely small particles require large field gradients to produce a significant degree of body force on the object. While this does not preclude the use of a magnetic propulsion system, the large power requirements necessary for implementation can be a significant engineering challenge. Guidance of larger scale particles and microdevices is a concept on which magnetic propulsion and guidance techniques are readily applicable. Previous investigations of magnetic guidance of microdevices have been undertaken. Once such study has been given in [9, 10]. In this work, an initial investigation to evaluate the feasibility of utilizing a conventional magnetic resonance imaging machine to impart large enough forces on a ferromagnetic sphere was performed. It was the intent of this study to assess the potential for employment of a microdevice in a human circulatory system, therefore the size of the sphere was chosen accordingly. The investigators utilized a 1/8” 1010 carbon steel sphere immersed in water and subjected to uniaxial flow in a 1/4” diameter tube. The MRI was utilized to keep the global position of the sphere static. It was determined that with a field gradient of 0.018T/m, the sphere could withstatnd a maximum flow of 0.370 ± 0.0064l/min. This value is comparable to that of the human circulatory system a good distance away from the heart, but still significantly larger than the flow rates in regions furthest from the heart. As previously stated, the fluid flow only was present in a single coordinate direction therefore evaluation of device guidance in multiple directions was 3 not achieved in this work. In addition, the static position control was achieved in an open loop fashion with the current in the electromagnet coils adjusted accordingly to arrive at equilibrium of the sphere in the flow. Feedback position control of a microdevice was not demonstrated in this work. A second novel concept of a microdevice guided by applied magnetic fields investigated in recent years is available in [15, 16]. This work studies motion of a “screw” type device propelled through both fluid and organic tissue by application of a rotating magnetic field. The device is essentially comprised of a cylindrical magnet with a spiral blade residing on its surface, analogous to a mechanical screw or auger. The cylindrical magnet is magnetized in a direction perpendicular to its central axis (a radial direction). The magnetization in this coordinate direction allows for a torque about the central axis of the cylinder upon application of a magnetic field perpendicular to the central magnet axis. The torque about the cylinder will be present until the magnetic moment of the cylinder aligns itself with the applied field. Subsequent rotation of the applied field results in a rotation of the cylinder about its central axis. The forward (or backward) motion of the device then occurs in the direction of the central cylindrical axis and is a product of the mechanical interaction of the spiral blade with the medium in which the device is immersed. The mechanism may be viewed much in the same way as that of an auger boring its way through a material body. Directional changes of the device motion are facilitated by a change in direction of the axis of rotation of the rotating magnetic field. Experimental investigations of this concept utilized a 15mm long device with a diameter of 1.2mm and a 0.15mm spiral blade with a pitch of 45 . Results of this study were quite interesting with maximum velocities achieved valued at approximately 20mm/s and 0.022mm/s in fluids with kinematic viscosities of 500mm2/s and 5 × 105mm2/s, respectively. Additionally, maneuvering of the auger device was accomplished in a mazed patterned waterway with 15mm wide passage and filled with a silicone oil having a viscosity of 1000mm2/s. Finally, additional experimental work has been performed evaluating the ability of a device such as this to “bore” its way through gel or a soft solid. Results in this area were quite satisfactory for the initial study into this technique. While this is a novel technique, it differs subtly from the method proposed in this work. The method proposed in [15, 16] utilizes an applied torque and a mechanical “propeller” effect to achieve device motion. This is different than the technique to be investigated by this work in which the direction of motion is determined by the applied field and the device is effectively pulled through the medium in which it resides. The body of work most closely related to the study presented in this investigation may be found in [17, 6, 7, 13, 18]. The cited research investigates a system known as the magnetic stereotaxis 4 system. This system is capable in theory of moving a small intracranial implant along a path of minimum neurological damage in the brain by utilizing applied field gradients. The insert is then positioned inside a deepseated tumor and inductively heated, hence heating and eliminating only the tumor. The envisioned implant would have a dimensional scale < 3mm. The initial investigations of this system incorporated a single electromagnet situated on the end of a robotic manipulator. The imaging method employed by this technique consist of a biplanar fluoroscope which was chosen due to its insensitivity to the applied magnetic fields. Further work into this apparatus emphasized the shortcomings of the moveable electromagnet technique and a system employing an array of stationary electromagnets was the selected path of continued study. More specifically, the stationary electromagnet assembly is a cubic arrangement on which the central axis of each electromagnet is collinear with the global coordinate system. This system is similar to the electromagnetic actuator technique explored in this work. Indeed the mechanism by which the magnetized seed is moved is the method studied. The target object is essentially pulled through the medium in which it resides. For this system however, all six electromagnets are energized. This operation technique results in an underdetermined system in current in which more variables are present than equations available to solve. As a result, optimization techniques must be employed to evaluate and search for a current solution satisfying the desired force on the target object. One advantage to this approach of utilizing all electromagnets is the ability to increase the field gradient at the target location, however it must be recognized that the overall system power is greatly increased due to utilization of all coils. In contrast, this work will examine a technique which will employ only the minimum number of coils required to generate the desired force on the target object. Another aspect of the magnetic stereotaxis system is that it is only capable of operation in an “open loop” mode. The target object is propelled in an impulsive fashion with the location being reached by just “nudging” the seed into position. This is presumably due to the difficulties associated with the optimal calculation of the current vector for the movement instance. By utilizing the minimal coil set method, it can be shown that a direct formulation giving the required coil currents can be made. This will be demonstrated with more detail in a later passage. 1.3 Contributions of this Study There are several contributions produced by this body of work. First, a detailed first principle model of the particle dynamics and electromagnetics will be presented. The electromagnetic model will be one formulated from a simplified coil representation consisting of finite current loops. Use 5 of this method allows for incorporation of the model into a controller in which closed loop control may be achieved. Secondly, a new technique based on the minimal coil set will be presented and explored. Advantages to this technique lie in the fact that a direct formulation may be made to calculate the required coil currents for producing the necessary or desired force on the target object. In addition, the coil switching algorithm will be rooted in a method by which the knowledge of the actual coil currents is unnecessary. Only the geometric field relations of the coils and the target object position is required. Thirdly, an experimental apparatus for investigation of this concept is designed and presented. The test configuration is of “bench top” scale and is optimally designed for minimal use of power. Finally, extensive closed loop experiments in which a small magnetized sphere is guided along 2D trajectories are performed and the results presented. As this is the initial work for achieving closed loop control in 2 dimensions utilizing electromagnetic actuators, a classical control technique is employed along with the minimal coil set current solution method. 1.4 Overview of this Study As has been previously stated, this study will limit itself to the simplified system of a magnetized sphere residing in a fluid medium. Additionally the sphere will be constrained to have only planar motion. Therefore, this study will examine the details associated with a system of four coil which are engineered to provide controlled motion of the sphere in the plane. This body of work will be presented in six separate chapters. Chapter 1 provides the motivation for the work and any existing relevant research. These details are available in the previous passages. The theoretical derivation of the system model will be presented in Chapter 2. In this chapter, the dynamic model will be derived from first principles. It describes the particle dynamics as well as the electromagnetics with the geometric field vector functions presented in a general form. Additionally, a state space representation will be presented. Also, the inverse current solution will be formulated. This will correspond to the system being studied which consists of the particle being constrained to a 2D plane. Subsequently, the technique by which the coil switching will be determined and the resulting calculated current trajectories smoothed will be discussed. A detailed electromagnetic field description will be presented in Chapter 3. This will be a more rigorous treatment of the general geometric vector field functions associated with the electromagnets which were introduced in Chapter 2. The coil field model will be based on a discrete filamentary loop description. Initially, the field resultant from a single filamentary loop will be presented, followed by a multiturn expansion. Subsequently, a formulation describing the net field and field gradient 6 generated by the quad coil assembly will be presented and discussed. Following that, relationships for calculating the electromagnet resistance, self inductance, and mutual inductance will be derived. Finally, experimental verification of the predicted field will be presented and discussed. The design and optimization of the experimental electromagnet assembly will be presented in Chapter 4. Calculation of the electromagnetic characteristics used in the optimization will employ the results of Chapters 2 and 3. A functional for minimization based on minimal power required for moving a magnetized sphere in a circular orbit is derived and presented. Additionally, a minimization technique based on steepest descent is discussed and employed. Finally, the details of the optimally designed electromagnetic coil set are presented. The results of closed loop control of the sphere’s 2D motion are presented in Chapter 5. A closed loop control law utilizing a classical PID controller with and without feedforward in concert with the inverse current solution is presented in Chapter 2. A numerical simulation is formulated and the results utilized to discuss further considerations about the coil switching technique. Additionally, three different path trajectories are examined in simulation. These consist of a spiral beginning at the origin and approaching a circular orbit, a linear path following the outline of a square, and the point regulation problem in which a single coordinate is provided and the path to approach it is not specified. Simulation results for both controllers with all three trajectories are presented and discussed. Finally, the control method was implemented in experiment and the results are presented and discussed. 7 CHAPTER 2 Dynamic System Model 2.1 Dynamic Model Development The equation of motion describing the movement of the target object (sphere) may be devised through a Newtonian formulation given here as m ¨~S (t) = ~FBias + ~FMag + ~FDrag, (2.1) where ~FBias denotes any bias force, ~FMag denotes the applied magnetic force, ~FDrag denotes a fluid drag force from the motion of the object relative to the fluid medium it resides in, m denotes the mass of the object, and ¨~S (t) denotes the vector acceleration of the object in an inertial reference frame. The corresponding trajectory ~S (t) describing the object’s motion may be ascertained by solving the differential equation (2.1). For the system consisting of a spherical target residing in a fluid medium, the first term on the right hand side of (2.1) may be defined to represent the resulting force of gravity on the sphere and the buoyant force of the fluid medium on the sphere ~FBias = ( sph − fl)Vsph~g, (2.2) where sph denotes the density of the sphere, fl denotes the density of the fluid medium in which the sphere resides, Vsph denotes the volume of the sphere, and ~g denotes the acceleration of the sphere due to gravitational attraction. The second term on the right hand side of (2.1) may be expressed by approximating the sphere as a magnetic dipole. This approximation for small objects is widely accepted and has seen some experimental correlation [2, 8]. When a magnetic field is applied, a resulting force and torque act on the dipole. In this case, assuming that the fluid medium has a permeability close to that of free space, the force and torque on the sphere are ~FMag(H) = 4 μ0r3 3 M· r H (2.3a) ~ Mag(H) = 4 μ0r3 3 (M×H), (2.3b) where M is the magnetization of the body, μ0 is the permeability of free space (valued at 4 × 10−7T · m/A), r is the radius of the sphere, and H is the applied magnetic field strength at the 8 location of the sphere [19, 20]. Note that boldface indicates a vector field. In the case of the particle residing in a medium where the relative permeability μfl 1, then μ0 ! μ0μfl. Recognizing that rotational resistance of the spherical object will be small in most fluid media, it is apparent from (2.3b) that the body will have a tendency to align with the applied field H. It may be noted that in some applications, particularly those involving an elongated magnetic object as suggested in [13], the alignment of the object relative to the applied force ~FMag may be considered as an additional constraint on the target. In this study, however, this requirement is disregarded due to the symmetry of the object. At this point in the model formulation, the type of material utilized for the spherical target must be addressed, as this affects the magnetization M achieved when the field H is applied. Primarily two types of materials are proposed for use, paramagnetic or ferromagnetic. The paramagnetic class of material is characterized by a linear relationship between the resulting magnetization and the internal magnetic field of the object. This linear constant is referred to as the magnetic susceptibility and is typically denoted as m. The resulting magnetic force on the paramagnetic sphere is [19] ~FMag,PM(H) = 2 μ0r3( m) ( m + 3) r H·H . The second class of material, ferromagnetic, is characterized by a nonlinear relationship between the applied magnetic field and the magnetization of the object. More specifically, the M −H curve is sigmoidal in shape indicative of these materials reaching a magnetic saturation. Additionally, the sigmoidal M − H curve may also possess a region of hysteresis about the origin depending on the classification of the ferromagnetic material as “hard” or “soft”. An example of a typical ferromagnetic M − H curve is available in Figure 2.1. Hard ferromagnetic materials are typified as those in which the M − H hysteresis loop is large and upon removal of the applied H, the sample is permanently magnetized (e.g. a permanent magnet). This is often referred to as magnetic remanence. In contrast, soft ferromagnetic materials have a small hysteresis loop and very little residual magnetization upon removal of any applied magnetic field (e.g. transformer iron) [21]. This work will address only the ferromagnetic material type. For the ferromagnetic sphere residing in a medium with permeability close to that of free space, the resultant magnetic force for a given amount of magnetization is ~FMag,FM(H) = 4 μ0r3 3 M H H· r H, (2.4) where M is the magnitude of magnetization and H is the magnitude of the applied field H at the sphere’s location [19]. For the permanently magnetized sphere, M may be treated as a constant. In 9 M H indicates initial magnetization Figure 2.1: Typical M − H curve for a ferromagnetic material. the Cartesian coordinate system, taking H = HXˆi + HYˆj + HZˆk whereˆi , ˆj , and ˆk are unit vectors in the respective coordinate directions, the magnitude of the applied field is given by H = H = q H2X+ H2 Y + H2Z . For the final force on the right hand side of (2.1), it has been readily accepted in the literature [1, 2, 3, 4, 5, 8, 10] that the flow about a spherical particle or object in these applications corresponds to a low Reynolds’ number condition (i.e. Re 1) and that the drag force may be formulated by applying the Stokes’ flow unbounded free stream drag force relation ~FDrag = 6 rμ(˙~S fl − ˙~S ), (2.5) where μ represents the viscosity of the fluid medium, ˙~S fl represents the free stream velocity and the quantity (˙~S fl − ˙~S ) denotes the velocity of the moving object relative to the free stream velocity. Indeed there is experimental evidence that, at least for the cases where particles are manipulated in a stagnant or slowly moving medium, this approximation correlates with the observed phenomena [2, 8]. It must be acknowledged, however, that flow scenarios may exist in which the Stokes approximation (2.5) may lose some accuracy. For instance, the Reynolds’ number for flow over a sphere 10 is Re = 2 fl ˙S flr μ , where fl is the density of the fluid. The extreme smallness of the particle may still allow for a calculation of Re 1 even with a significant magnitude of flow velocity around the sphere. If the object is of a larger size (e.g. a microdevice or catheter guide), the result of Re 1 may no longer hold, pushing the drag coefficient into a nonlinear function of Re. Additionally, the relative size of the sphere to that of the channel in which it is moving and the geometry of that channel will have an impact on the calculation of the drag coefficient. As an example, consider the system of a sphere moving along the central axis of a cylindrical duct or a square duct. For the cylindrical duct, the drag force on the sphere is given as FDrag,cyl = 6 μUr 1 + k r R0 , (2.6) where U is the velocity of the sphere, R0 is the radius of the cylindrical duct, and k is a correction coefficient valued at k = 2.10444 [22]. For the square duct, the drag force on the sphere is given by FDrag,sqr = 6 μUr 1 + k r l , (2.7) where l is the half width of the square duct and k = 1.903266 [22]. Likewise, the wall effects on an object in close proximity to the channel surface have been neglected. These will influence the drag force and often impart lift and moment components on the sphere. Depending on the flow boundary conditions, often the most convenient way to evaluate these effects is through numerical simulation [23, 24]. Additionally, it is well known that fluids in living bodies do not behave in a “Newtonian” fashion. Additional forces may be present, i.e. collisions with other small bodies. Depending on the size of the target object, these dynamics may need to be accounted for. Finally, the condition of turbulent flow in the medium is neglected. Turbulent flow of the fluid medium will lead to significant deviation from the Stokes’ approximation. The wall and turbulent flow effects will be ignored for the purposes of this study and left for analysis at a later time. For now, equations (2.2), (2.4), and (2.5) may be used in (2.1) to solve for the motion of the sphere. A component which remains to be addressed is the method to generate the required field H in (2.4). The technique utilized in this study for generation of the magnetic field on the sphere is that of an air core electromagnet. It will be assumed that the required motion of the target and hence the required field changes will be slow enough that the magnetic quasistatic approximation may be used to ignore any transient effects (e.g. eddy currents, etc.). In this system, stationary electromagnets will be employed. Because of this, multiple coils must be used to facilitate parallel and antiparallel 11 forces along multiple coordinate directions. This is due to the fact that when the sphere is unrestrained in the rotational sense and allowed to freely align with the field, a single electromagnet will act in an attractive capacity alone. For three dimensional motion an obvious configuration is to make use of six coils arranged in a cubic fashion similar to that of the magnetic implant guidance system by Stereotaxis, Inc. [6, 7, 13]. The resulting field from a closed path current density conforms to a solution of Maxwell’s field equations which are readily available in [19, 20]. Application of these equations and subsequent reduction of this set to the law of BiotSavart yields a relation expressing the field H of the Pth individual coil as a linear function of the current applied to the coil HP = jP [GXˆi + GYˆj + GZˆk] = jPGP (X, Y,Z), (2.8) where jP is the applied current,ˆi , ˆj , and ˆk are unit vectors in the respective coordinate directions, and Gp is a nonlinear vector field multiplier, which is a function of the location of the sphere relative to the coil and the geometry of the coil itself. The combined field from the employment of Q electromagnets may be ascertained through superposition of the individual coil fields H = [J]T 2 66664 G1 ... GQ 3 77775 = [J]T [G], (2.9) where [J] is a Q × 1 column vector of the applied currents ([J] 2 RQ) and GP is the geometry and position dependent vector field from the Pth coil evaluated in an inertial coordinate system. Using (2.9) the following quantity in (2.4) is derived M H H· r H =M H [J]T [kFM,X]ˆi + [kFM,Y ]ˆj + [kFM,Z]ˆk [J] =M H [J]T [~kFM][J], (2.10) where kFM,X = [GX] @GX @X T + [GY ] @GX @Y T + [GZ] @GX @Z T , kFM,Y = [GX] @GY @X T + [GY ] @GY @Y T + [GZ] @GY @Z T , kFM,Z = [GX] @GZ @X T + [GY ] @GZ @Y T + [GZ] @GZ @Z T , H = r [J]T h [GX][GX]T + [GY ][GY ]T + [GZ][GZ]T i [J], 12 [GS] = 2 66664 GS,1 ... GS,Q 3 77775 , (2.11a) @GS @S T = @GS,1 @S , · · · , @GS,Q @S . (2.11b) Note that in equations (2.11a) and (2.11b), S = X, Y,Z. All field components and partial derivatives of the field components are evaluated in a global inertial coordinate system. Substituting (2.2), (2.4), (2.5), and (2.10) into (2.1) yields the equation of motion m ¨~S = ( s − fl)Vs~g + 4 μ0r3 3 M H [J]T [~kFM][J] + 6 rμ(˙~S fl − ˙~S ). (2.12) Observe that (2.12) relates a given column vector of input current to the resultant motion of the sphere. Because this system utilizes electromagnet coils, the inductive effects of the electromagnet array must be taken into account. Hence, the equation of motion must be augmented by a suitable description of the electrical dynamics for a complete description of the system. Approximating each coil as an RL circuit (a resistor in series with an inductor) results in the following relationship between the input voltage and current flow of the electromagnets vm = Rmjm + XQ n=1 Lm,n djn dt , (2.13) where Q is the number of coils, m = 1, 2, . . . ,Q, n = 1, 2, . . . ,Q, vm denotes the input voltage of coil m, Rm denotes the resistance of coil m, jm denotes the current in coil m, and djn/dt denotes the time rate of change of the current in coil n. Additionally, Lm,n signifies the mutual magnetic inductance of coil n on coil m or the self inductance of coil m when n = m [19]. It should be noted that while the self inductance is always positively valued, the sign of the mutual inductance is dependent on the orientation of coil m to n. Observing that (2.13) is a linear system, it may be rewritten in state space form [J˙] = −[L−1][R][J] + [L−1][V ] (2.14) where [V ] and [J] are Q × 1 column vectors of the voltages and currents in Q coils, [R] is a Q × Q diagonal matrix of the coil resistances, and [L] is a Q×Q symmetric matrix of inductance with the diagonal elements being the self inductance terms and the off diagonal elements being the mutual inductance terms. Note that [L] is symmetric because Ln,m = Lm,n and also note that [L] is invertible [19]. Equation (2.12) is used in concert with (2.14) to provide a full description of the sphere subject to Stokes’ flow and a magnetic field generated by Q discrete electromagnets with input voltages [V ]. 13 Z X Y x1 x3 x2 x4 y1 z1 y2 z2 z y 4 4 y3 z3 plane of motion (YZ) j1 j2 j3 j4 1 4 3 2 c0 sphere Figure 2.2: Schematic of a four coil arrangement to provide forces for planar motion of a sphere. The individual coil coordinates indicate direction of positive current polarity (right hand rule). 2.2 Inverse Current Solution for Constrained 2D Motion with a Reduced Coil Set In an effort to further explore a system using stationary coils to manipulate the position of a magnetized sphere in a fluid medium, this study will examine the two dimensional planar motion case with the sphere constrained to the Y −Z plane (see Fig. 2.2). As previously mentioned, a cubic arrangement of coils can facilitate three dimensional motion of the target sphere [6, 7, 13]. Similarly, a set of four coils arranged on the edges of a plane with the central axes of the coils lying in the plane will facilitate forced motion in that plane. A schematic of this type of system is shown in Fig. 2.2. The inverse current solution entails determining the current trajectories [J] which would cause the sphere to move along a desired motion trajectory ~S . Indeed, the inverse current solution may be considered as a more thorough examination of (2.4), acknowledging in this case that ~FMag is the force required to “offset” the inertial, bias, and drag forces incurred by the sphere moving along ~S . It is obvious from (2.4) that ~FMag,FM is a function linear in [J] but is nonlinear in the components of H, so a direct algebraic solution is not feasible. With this noted, it is evident that in the configuration where the sphere is constrained to a central plane and when the coil configuration to provide both parallel and antiparallel motion is employed, the system is underdetermined. In other words, a 14 θ FZ FY fMag,FM,Y fMag,FM,Z FMag,FM Figure 2.3: Graphical representation of the applied planar magnetic force. solution in four variables (the coil currents) is to be ascertained from just two equations (the vector components of (2.4) in the Y and Z directions). Due to the fact that electromagnets will act in an attractive capacity alone, one may reduce the number of utilized coils to only those which are necessary for the required force. In the scenario of planar motion utilizing the system in Fig. 2.2, this entails reducing the active coils to that of “corner” combinations (i.e. coils 1&4, 1&3, 2&4, and 2&3), commensurate with utilizing only adjacent coils. This will reduce the number of solution variables from four to two and make the system determined. For example, suppose that there is a known required magnetic force ~FMag,FM in the Y − Z plane. This force may be represented in a cylindrical coordinate system as ~FMag,FM = fMag,FM[cos( )ˆj + sin( )ˆk] where fMag,FM is the force magnitude and gives the angular direction of the force measured from the (global) Y axis with the counter clockwise direction taken as positive (see Fig. 2.3). Recognizing that the coils will only attract the sphere, for a required force in the direction of = 45 , the necessary coils to achieve this will be numbers 1 and 4 (see Fig. 2.2) with the actual current amplitudes and polarity dependent on the location of the sphere. In other words, the coil combination for a force pointing into a given quadrant will be the coils which bound that quadrant. One may postulate that the corner combinations of 1&4, 1&3, 2&3, and 2&4 correspond to force directions of 0 90 , 90 180 , 180 270 , and 270 360 , respectively. However, this quadrant 15 assignment to pairs of coils is approximate due to the nonlinear nature of the attractive magnetic fields. The components of the required force ~FMag,FM may be expressed as a function of one another fMag,Z = tan( )fMag,Y , (2.15) where once again is the angular direction of the required force measured from the Y axis with the counter clockwise direction taken as positive. As can be seen from (2.4) and (2.12), the Y and Z components of magnetic force for the ferromagnetic sphere are fMag,Y = 4 μ0r3 3 M H [J]T [kFM,Y ][J], (2.16a) fMag,Z = 4 μ0r3 3 M H [J]T [kFM,Z][J]. (2.16b) Substituting (2.16a) and (2.16b) into (2.15), rearranging, and then simplifying yields [J]T [ ][J] = 0, (2.17) where [ ] = 8>>< >>: [kFM,Y ], for = 2 ± n , n = 0, 1, 2. . . . ; [kFM,Z] − tan( )[kFM,Y ], otherwise. The quantity [J] is a 2×1 vector of the currents in the reduced coil set and the vector fields [kFM,Y ] and [kFM,Z] are 2×2 matrices containing the G and @G/@S terms of the reduced set in accordance with (2.11a) and (2.11b). The matrices [kFM,Y ] and [kFM,Z] are given by the quantities in (2.10), however GS and @GS/@S are now 2 × 1 matrices with the first element corresponding to coil 1 or 2 and the second element corresponding to coil 3 or 4, depending on the coil combination being examined. Once again note that all field components and partial derivatives of the field components are evaluated in the global inertial coordinate system. Equation (2.17) corresponds to a quadratic surface in the J (current) space, centered about the origin. For a current solution other than J = [0 0]T to exist, [ ] must be indefinite or semidefinite [25]. In other words, for a nontrivial solution to exist [ ] cannot be positive or negative definite. Several techniques are available to ascertain the definiteness of [ ]. One method involves examining the signs of the principle minors of [ ]. If the principle minors are not exclusively < 0 or > 0, then [ ] is semidefinite or indefinite [25]. Alternatively, for the 2×2 matrix [ ] one may examine the product of eigenvalues. If the eigenvalue product is 0, then [ ] is semidefinite or indefinite. To facilitate the inverse current solution, it is necessary to require the eigenvalues to be real. This may be achieved by forcing symmetry through the relation [J]T [ ][J] = [J]T [( + T )/2][J] = [J]T [ ][J]. 16 Recognizing that the 2×2 matrix [ ] of the reduced coil set will yield two eigenvalues, the eigenvalues of [ ] will be denoted as 1 and 2. For a sphere located in the Y −Z plane, a force may be effected on the sphere by any coil combination with an eigenvalue product of 1 2 0. As an example, the sign of the eigenvalue products for coil combination 1&4 was examined over a range of −45 135 for three different locations in the Y −Z plane. The initial system analyzed and demonstrated here consists of four rectangular filamentary current loops with a major dimension of 1 unit, a minor dimension of 0.5 units, and an offset “c” from the global origin (see Fig. 2.2) of 1.5 units. The effective force angle range evaluated at the Y − Z locations of (0.5, 0.5), (0.0), and (−0.5,−0.5) are −30.51 120.51 , −2.43 82.43 , and 10.89 79.11 , respectively. These angular ranges indicate the directions of force that coil combination 1&4 is capable of achieving at the given sphere locations. As alluded to earlier, the range of force directions a particular coil combination may impart on the sphere is dependent on the nonlinear field expressions which are functions of the sphere’s position as well as coil geometry. The effective range of force direction decreases as the sphere moves away from the electromagnet pair. To solve for the required currents to achieve the desired (planar) force ~FMag,FM, the orthonormal basis [eQ ] composed of the normalized eigenvectors of [ ] may be used to form the orthogonal transformation [J] = [eQ ][ e J]. Application of this to [J]T [ ][J] = 0 results now in the diagonal quadratic [ e J]T [eD ][ e J] = 0 where [eD ] is a diagonal matrix of the eigenvalues of [ ]. In the e J space, two line loci may now be formed from the diagonal quadratic ˜j 2 i 1 +˜j 2 ii 2 = 0 8>>< >>: ˜j ii = ±˜j i r − 1 2 ; 2 6= 0, (2.18a) ˜j i = ±˜j ii r − 2 1 ; 1 6= 0, (2.18b) where the subscripts i = 1, 2 and ii = 3, 4 denote which coil is used in the calculation. As a reminder, recall that [ ] is also a function of which coils are included in the reduced coil set. Note that (2.18) actually provides two separate formulations for the line loci with a caveat on the existence of a 17 nonzero eigenvalue. Utilizing (2.18), the [ e J] space current vector may now be rewritten as [ e J] = 8>>>>>>>>< >>>>>>>>: ˜j i 2 64 1 ± q − 1 2 3 75 =˜j i[ ± 1 ]; 2 6= 0, (2.19a) ˜j ii 2 64 ± q − 2 1 1 3 75 =˜j ii[ ± 2 ]; 1 6= 0. (2.19b) Applying the orthogonal transformation to (2.16a) yields fMag,Y = 4 μ0r3 3 M eH [ e J]T [eQ ]T [kFM,Y ][eQ ][ e J], (2.20) where eH = q [ e J]T [eQ ]T [G][eQ ][ e J], [G] = GXGT X + GY GTY + GZGT Z. Substitution of (2.19a) into (2.20), simplifying, and utilizing (2.18a) results in the [ e J] solution of ˜j i = ± 3fMag,Y 4 μ0r3M q [ ± 1 ]T [eQ ]T [G][eQ ][ ± 1 ] [ ± 1 ]T [eQ ]T [kFM,Y ][eQ ][ ± 1 ] , (2.21a) ˜j ii = ±˜j i r − 1 2 . (2.21b) Note that (2.21) must still conform to the caveat of 2 6= 0. For the condition that 2 = 0, a similar process may be employed utilizing (2.18b) and (2.19b) to provide the solution ˜j i = ±˜j ii r − 2 1 , (2.22a) ˜j ii = ± 3fMag,Y 4 μ0r3M q [ ± 2 ]T [eQ ]T [G][eQ ][ ± 2 ] [ ± 2 ]T [eQ ]T [kFM,Y ][eQ ][ ± 2 ] , (2.22b) where 1 6= 0. Additionally, solutions may be formulated using (2.16b). This will result in forms similar to (2.21) and (2.22) with fMag,Y ! fMag,Z and [kFM,Y ] ! [kFM,Z]. To select the proper signs for the solution coordinates, the correct groupings on loci intersections must be observed. These are available in Table 2.1, where in the solution value˜j ( , = ±, = i, ii), = ± denotes the sign just to the right of the “=” symbol in (2.21) or (2.22) and = ± denotes the sign used for [ ], as defined in (2.19). It can be seen from (2.21), (2.22), and Table 2.1 that for the planar motion reduced set, there are four different solution coordinates. The coordinate system e J is a transformed system and therefore has no physical meaning. To complete the calculation of the current solution, the e J coordinate solutions must be transformed back to the J space via [J] = [eQ ][ e J]. 18 Table 2.1: Correct Line Loci Groupings for e J Space Solution Coordinates solution set sign e Ja (˜j ++ i ,˜j ++ ii ) e Jb (˜j −+ i ,˜j −+ ii ) e Jc (˜j +− i ,˜j −− ii ) e Jd (˜j −− i ,˜j +− ii ) A plot of the intersecting loci indicating the solution coordinates in the J space for the previously mentioned filamentary loop configuration, a sphere with a unit parameter multiplier (4 μ0r3M/3 = 1) located at the global origin, and a required magnetic force of ~FMag = 1ˆj + 1ˆk is shown in Fig. 2.4. As can be seen, the four solution coordinates may be ascertained by the intersection of the component loci of equation(s) (2.16) or the intersection of the line loci for equation (2.17) with just one of the component loci of (2.16). Indeed it is the latter technique which yields the results of (2.21) and (2.22). The calculated solution coordinates using (2.21) for solutions Ja, Jb, Jc, and Jd are (11.208, 11.208), (−11.208,−11.208), (−25.519, 25.519), and (25.519,−25.519), respectively. As can be seen, they correlate with the loci intersections. The inverse current solution presented in this work was derived for the ferromagnetic material type. A similar derivation may be used for the paramagnetic material and will give the same number of solutions, four. However, the solution is of a simpler nature because the shapes of the component loci for a given force are analytically known. In this case the current solution coordinates correspond to intersecting conic sections. Extending this line of thought would result in the inverse current solution for the paramagnetic material being manipulated by a cubic arrangement of coils being the intersection of quadric surfaces. While the inverse current solution given in relations (2.21) and (2.22) do indeed result in coil currents providing the desired magnetic force on the sphere, as will be shown the sign convention given in Table 2.1 is not sufficient to provide a smooth current trajectory as a function of time. It is advantageous however for the current profile to be smooth to minimize inductive effects incurred by discontinuities in the desired current profile. Further analysis into this formulation is necessary to achieve this and will be provided in the next section. 19 −50 −40 −30 −20 −10 0 10 20 30 40 50 −50 −40 −30 −20 −10 0 10 20 30 40 50 j 1 j 4 line loci f Mag,Z locus f Mag,Y locus Figure 2.4: Solution loci for the 1&4 coil combination in the J space. 2.3 Current Trajectory Smoothing Technique To formulate a more detailed analysis of the inverse solution current coordinates using (2.21) and (2.22), a more precise representation of the desired force on the sphere is needed. As previously stated, the required force on the sphere may be regarded as the force required to overcome the inertial, drag, and buoyant forces acting on the sphere. From (2.12), the required force to be provided by the electromagnets for a given trajectory ~S is ~FMag,FM = m ¨~S + 6 rμ(˙~S − ˙~S fl) + ( fl − s)Vs~g. (2.23) In addition, a trajectory must be examined to fully understand the resultant current profiles calculated using the inverse current formulation presented. For this example, the trajectory to be used was chosen to be a circular path in the Y − Z plane, orbiting the origin with a frequency of 0.1Hz. This may be defined as ~S (t) = cos 2 10 t ˆj + sin 2 10 t ˆk , (2.24) where is the radius of the path and t represents time. As can be seen, the first and second time derivatives of (2.24) may be easily evaluated and substituted into (2.23) to calculate the necessary force ~FMag,FM. It will be assumed that the vector acceleration due to gravity, denoted as ~g, will be 20 normal to the plane of motion, hence offset by a constraint force and may be subsequently dropped from (2.23). Furthermore, it will also be assumed that this testing apparatus will move the sphere in a static fluid, allowing for the dropping of the term ˙~S fl in (2.23) as well. The physical parameters in (2.21) (or (2.22)) and (2.23) which need assignment consist of the mass of the sphere, the radius of the sphere, the viscosity of the fluid medium, the radius of the circular path, and the magnetization of the sphere. In an effort to formulate a more “realistic” example, parameters closer to that of the physical system will be utilized. For this system, m = 0.1257g, r = 1/8”, = 0.5”, and μ = 0.950. It should be noted that the fluid selected for initial testing is glycerine and the viscosity value is consistent with that found in many fluid property tables. The sphere used in this study is a permanent neodymium rareearth magnet of grade 40. Accounting for demagnetization effects of the spherical shape results in a permanent magnetization of M = 3Br 2μ0 , where Br is the remanent flux density in the magnetized sphere [19, 21]. For this study, the remanent flux density will be estimated to have the value of Br = 1T. A selection of this value is consistent with the material type and grade of magnet. Finally, the coil geometry used for this example is that of rectangular coils wound with 12 gauge magnet wire. The coils consist of 21 layers of wire with 18 turns per layer. The inner dimensions of the coils are 0.439” × 0.757”, while the distance from each coil face to the global origin is 2.5”. A more in depth analysis of this coil configuration will be provided in Chapter 4. The trajectory chosen for this study is that of a circle orbiting the origin with a constant angular velocity. Because of this, symmetry may be used and only one coil combination need be examined. In other words, the force requirements for each coil combination is the same, hence the current solutions are of the same form. Utilizing this, the combination of coils 1 and 3 will be examined initially. A plot of the four solutions given by (2.21) (see Table 2.1) for coil combination 1 and 3 is shown in Figure 2.5. As can be seen, the solution provided by (2.21) for coils 1 and 3 has two distinct regions in which [ 13] is semidefinite or indefinite. However, if this solution is used to calculate the resultant force generated and then compared to the required force it can be seen that one of the current regions is erroneous. This is demonstrated in Figure 2.6. The incorrect solution regions of both force components seen in Figure 2.6 are artifacts of the mathematical formulation of the inverse solution technique. As is evidenced by Figure 2.6, the requirement of [ ] to be semidefinite or indefinite is not “strong enough” to result in a true solution. The resulting inverse solution must also satisfy the desired direction of the applied magnetic force. In other words, calculated must 21 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J a (A) j 1 j 3 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J b (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J c (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J d (A) time (sec) Figure 2.5: Initial current solution for coil combination 1 and 3 during circular trajectory. 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,a (N) f req f mag 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,a (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,d (N) time (sec) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,d (N) time (sec) Figure 2.6: Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces. 22 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J a (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J b (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J c (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J d (A) time (sec) j 1 j 3 Figure 2.7: Correct current solution for coil combination 1 and 3 during circular trajectory. equal desired. Disregarding the solution regions which do not satisfy this second criteria, Figures 2.7 and 2.8 may be constructed which give the correct solution and calculated resultant forces for coil combination 1 and 3. As can be seen, the region of incorrect force results has vanished from the calculated solution. Now that a known correct solution with no erroneous sections has been formulated, a second facet of this result may be observed in Figure 2.7, that being the discontinuities present in the current trajectories. Indeed, this aspect when examined at face value may be construed as quite disturbing. While this current trajectory would in fact result in the necessary magnetic force being imparted on the sphere, the ability to achieve this trajectory in an inductive system may or may not be entirely possible. Because of this, it would be desirable for the current trajectories to be “smooth” or continuous. Upon examination of the discontinuity shown in Figure 2.7, it has been observed that the discontinuous point corresponds to the location at which the transformation [eQ ] changes from a rotation to a reflection, or vice versa. It is well known that the particular type of transformation (either rotation or reflection) may be ascertained by examining the sign of the determinant of the transformation [26, 25, 27]. Defining D = det[eQ ] and plotting sgn(D1,3) results in Figure 2.9. A comparison of Figures 2.7 and 2.9 shows that the location of the current solution 23 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,a (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,a (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,d (N) time (sec) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,d (N) time (sec) f req f mag Figure 2.8: Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces are now absent. 0 1 2 3 4 5 6 7 8 9 10 !1.5 !1 !0.5 0 0.5 1 1.5 ,./012/34 25617 183 4 Figure 2.9: A plot of sgn(D1,3), where D = det[eQ ]. A value of +1 indicates that transformation is a rotation, while a value of −1 is indicative of a reflection. 24 Table 2.2: Modified Sign Groupings for Smooth e J1,3 = ˜j 1 1 1 ,˜j 3 3 3 Solution Coordinates solution set 1 1 3 3 e Ja,1,3 +sgn(D1,3) +sgn(D1,3) + +sgn(D1,3) e Jb,1,3 −sgn(D1,3) +sgn(D1,3) − +sgn(D1,3) e Jc,1,3 +sgn(D1,3) −sgn(D1,3) − −sgn(D1,3) e Jd,1,3 −sgn(D1,3) −sgn(D1,3) + −sgn(D1,3) discontinuity corresponds to the point at which the transformation eQ 1,3 changes from a rotation to a reflection. This dependence may be accounted for by modifying the definitions of the correct solution signs given in Table 2.1. A “smoothed” version of the solution given in Figure 2.7 may be achieved by the solution signs given in Table 2.2. Note that the pattern of ±’s preceding the signum functions matches that seen in the sign pattern given in Table 2.1. Utilizing the signing pattern of Table 2.2, a smooth current solution for coil combination 1 and 3 may be formed. This is shown in Figure 2.10. As can be seen, there are two types of solution, one in which the coil currents of like sign and one in which the coil currents are of opposite sign. Additionally, for each of the two types of solution there are two forms in which the current polarity of each coil is reversed. As stated earlier, for the generic circular trajectory chosen here, the current trajectory solutions for the other coil combinations will be of similar form to that shown in Figure 2.10. However, the sign convention for providing a smooth current trajectory in each coil combination will differ from that given in Table 2.2. For the remaining coil combinations 2&3, 2&4, and 1&4, the resulting inverse solutions using equation (2.21) and the sign convention given in Table 2.1 results in the current trajectories shown in Figures 2.11, 2.12, and 2.13. In addition, the required sign convention needed for each of the combinations to achieve a “smooth” trajectory are given in Tables 2.3, 2.4, and 2.5, with the resulting smooth trajectories also shown in Figures 2.11, 2.12, and 2.13. As can be seen, the discontinuities have been removed by using the sign conventions provided in Tables 2.2, 2.3, 2.4, and 2.5. It should be noted that unlike the solution from equation (2.21), the sign convention given in the tables is not dependent on the solution form used. It has been previously stated that (2.21) is equally valid by substituting fMag,Y ! fMag,Z and [kFM,Y ] ! [kFM,Z]. In the case of sign convention, fMag,Y is utilized in the formulation regardless of the use of the Y direction solution versus the Z. Now that smooth current trajectories have been calculated, it is matter of 25 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J a (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J b (A) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 J c (A) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 J d (A) time (sec) j 1 j 3 Figure 2.10: Smoothed current solution for coil combination 1 and 3 during circular trajectory. selecting from the four different solutions provided by (2.21). Any of the four solutions presented will result in the desired force, however one may choose to prefer one solution over another based on the polarity or the total power utilized. For the purposes of experiment design and optimization of this design, symmetry may be leveraged and only analysis of the hardware configuration required to generate one coil combination’s solution is necessary. As will be seen in Chapter 4, only the solution for coil combination 1&3 shown in Figure 2.10 will be used in the optimization routine. 2.4 State Space Model and Closed Loop Control Considerations Implementation of closed loop control on this type of system is challenging and has not been explored in literature yet. Even with the reduced set and the exact calculation of the required currents, it has been shown in this work that for the two dimensional case there will be four possible solutions. For the three dimensional case utilizing six stationary coils, the reduced set would number three and the number of solution coordinates would be eight. Because of this, closed loop control employing the reduced set will be challenging as well. Even with that said, the benefits of being able to exactly calculate the possible solutions are obvious when compared with using optimization methods to search out solutions for an underdetermined system. Feedback control will be further complicated 26 0 2 4 6 8 10 !5 0 5 J a (A) 0 2 4 6 8 10 !5 0 5 J b (A) 0 2 4 6 8 10 !5 0 5 J c (A) 0 2 4 6 8 10 !5 0 5 J d (A) time (sec) 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !1 0 1 time (sec) j 2 j 3 j 2 j 3 Figure 2.11: Initial and smoothed current solution for coil combination 2 and 3 during circular trajectory. 0 2 4 6 8 10 !5 0 5 J a (A) 0 2 4 6 8 10 !5 0 5 J b (A) 0 2 4 6 8 10 !5 0 5 J c (A) 0 2 4 6 8 10 !5 0 5 J d (A) time (sec) 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !5 0 5 time (sec) j 2 j 4 j 2 j 4 Figure 2.12: Initial and smoothed current solution for coil combination 2 and 4 during circular trajectory. 27 0 2 4 6 8 10 !5 0 5 J a (A) j 1 j 4 0 2 4 6 8 10 !5 0 5 J b (A) 0 2 4 6 8 10 !5 0 5 J c (A) 0 2 4 6 8 10 !5 0 5 J d (A) time (sec) 0 2 4 6 8 10 !5 0 5 j 1 j 4 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !1 0 1 time (sec) Figure 2.13: Initial and smoothed current solution for coil combination 1 and 4 during circular trajectory. Table 2.3: Modified Sign Groupings for Smooth e J2,3 = ˜j 2 2 2 ,˜j 3 3 3 Solution Coordinates solution set 2 2 3 3 e Ja,2,3 +sgn(fMag,Y ) +sgn(D2,3fMag,Y ) +sgn(D2,3) +sgn(D2,3fMag,Y ) e Jb,2,3 −sgn(fMag,Y ) +sgn(D2,3fMag,Y ) −sgn(D2,3) +sgn(D2,3fMag,Y ) e Jc,2,3 + −sgn(D2,3fMag,Y ) −sgn(D2,3fMag,Y ) −sgn(D2,3fMag,Y ) e Jd,2,3 − −sgn(D2,3fMag,Y ) +sgn(D2,3fMag,Y ) −sgn(D2,3fMag,Y ) 28 Table 2.4: Modified Sign Groupings for Smooth e J2,4 = ˜j 2 2 2 ,˜j 4 4 4 Solution Coordinates solution set 2 2 4 4 e Ja,2,4 + +sgn(D2,4) +sgn(D2,4) +sgn(D2,4) e Jb,2,4 − +sgn(D2,4) −sgn(D2,4) +sgn(D2,4) e Jc,2,4 + −sgn(D2,4) −sgn(D2,4) −sgn(D2,4) e Jd,2,4 − −sgn(D2,4) +sgn(D2,4) −sgn(D2,4) Table 2.5: Modified Sign Groupings for Smooth e J1,4 = ˜j 1 1 1 ,˜j 4 4 4 Solution Coordinates solution set 1 1 4 4 e Ja,1,4 +sgn(D1,4) +sgn(D1,4fMag,Y ) +sgn(fMag,Y ) +sgn(D1,4fMag,Y ) e Jb,1,4 −sgn(D1,4) +sgn(D1,4fMag,Y ) −sgn(fMag,Y ) +sgn(D1,4fMag,Y ) e Jc,1,4 +sgn(D1,4fMag,Y ) −sgn(D1,4fMag,Y ) − −sgn(D1,4fMag,Y ) e Jd,1,4 −sgn(D1,4fMag,Y ) −sgn(D1,4fMag,Y ) + −sgn(D1,4fMag,Y ) 29 for this system due to the inductive effects of the coils. For a given (slow) trajectory, the current solutions to follow that path may be calculated from (2.21) or (2.22). The voltage solutions would then need to be evaluated using (2.13). The method of using a reduced coil set would inherently involve “switching” coils on and off. Inductive effects associated with this type of technique must be accounted and compensated for. One advantage of this system is that for the electrical dynamics, full state feedback is available. Therefore it is conceivable that the switching effects may be taken into account. In the case of large coils with significant inductance parameters, full state feedback of current may prove to be invaluable when dealing with the switching effects. To synthesize a controller, it is beneficial to first express the entire system model in a state space form. For this system, the state variables will be designated as the sphere’s Y and Z positions, the velocity components ˙Y and ˙Z , and currents in the four coils j1, j2, j3, and j4. The inputs for this system are the coil voltages v1, v2, v3, and v4. This results in the state vector x = Y ˙Y Z ˙Z j1 j2 j3 j4 T , (2.25) and the input vector u = v1 v2 v3 v4 T . (2.26) It will be assumed that the measured outputs are the sphere positions Y and Z, and the four coil currents. As can be seen from (2.12) and (2.13), this system may be represented as the two linear systems of sphere dynamics and electrical dynamics coupled by the nonlinear magnetic force terms. In addition any unknown drag force which deviates from the linear Stokes’ static drag may be treated as a disturbance. This leads to the state space representation of x˙ = [A]x + f(x) + [B]u + [D], (2.27a) y = [C]x, (2.27b) where [A] = 2 64 As [;]4×4 [;]4×4 −L−1R 3 75 , [B] = 2 64 [;]4×4 L−1 3 75 , As = 2 66666664 0 1 0 0 0 −6 rμ m 0 0 0 0 0 1 0 0 0 −6 rμ m 3 77777775 , [D] = 2 66666666664 0 dY 0 dZ [;]4×1 3 77777777775 , 30 [C] = 2 64 Cs [;]4×4 [;]4×4 [I]4×4 3 75 , Cs = 2 66666664 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 3 77777775 , f(x) = 2 66666666664 0 4 μ0r3 3m p M xT x xT [KFM,Y ]x 0 4 μ0r3 3m p M xT x xT [KFM,Z]x [;]4×1 3 77777777775 , [KFM,S] = 2 64 [;]4×4 [;]4×4 [;]4×4 [kFM,S] 3 75 , = 2 64 [;]4×4 [;]4×4 [;]4×4 [G] 3 75 . Note that in (2.27), [;]m×n represents an m×n matrix of zeros, [I]m×m represents an m×m identity matrix, and dS, S = X, Y,Z, represents the unknown component flow disturbances. Equation (2.27) is an eighth order nonlinear state space model which may be used for controller synthesis and design. As can be seen, the system is linear in the input u but contains rather formidable nonlinear terms describing the field forces. With the system now represented in a state space form, the formulation of the field components as they pertain to the coil geometries and the calculation of the mutual inductances must be examined in greater detail. This will be explored in the following chapter. 31 CHAPTER 3 Electromagnetic Field Description 3.1 The Single Filamentary Loop While a detailed examination of the planar equation of motion and the relation between all the fields of the individual coils is developed in Chapter 2, the functions describing the resulting field components and gradient field components generated by the rectangular electromagnets is not addressed. Greater insight into this physical system allows for the creation of a model which may then be used in a predictive capacity and allow for design and optimization of an experimental setup for evaluation of this concept. This chapter will expand on the work of Chapter 2 and provide the derivation of a model describing the field generated by elliptical coils wound on a rectangular core, represented as rectangular electromagnets. This geometric simplification allows for the analytic formulation of the descriptive field equations. In addition, the technique used will entail the approximation of the coils as assemblies of discrete filamentary loops [28]. This technique has been selected due to the less difficult analytical treatment when compared to the volumetric coil model [29, 30]. The law of BiotSavart will be used to calculate the field contribution of each loop. To develop a relation describing the resulting magnetic field from the quad coil configuration, an analysis of a single filamentary loop will initially be undertaken. The resulting magnetic field for a single (closed) filamentary loop may be ascertained by using the law of BiotSavart [31], h = j 4 I d~l ×~r r3 , (3.1) where h denotes the magnetic field at a point in 3D space, j denotes the applied current in the filamentary loop, d~l denotes a differential element of the loop, ~r denotes the vector from the differential element to the point at which the vector field is being evaluated, and r is the magnitude of ~r. In this study, we will first examine a single rectangular filamentary loop as represented in Figure 3.1. Application of (3.1) in a piecewise fashion to the rectangular conductor results in h = j 4 " Z b −b d~l 1 ×~r1 r3 1 + Z a −a d~l 2 ×~r2 r3 2 + Z b −b d~l 3 ×~r3 r3 3 + Z a −a d~l 4 ×~r4 r3 4 # , (3.2) where d~l m is a differential element of the m side of the rectangular filament, ~rm is the vector from 32 z y x b a 1 2 3 4 j Figure 3.1: Schematic of a single filamentary conducting loop. The local coordinate system is fixed to the center of the loop with the central axis aligned with the z coordinate direction. Positive current flow obeys the right hand rule. The numbers indicate sections of conductor while the letters denote the dimensions of the rectangle. the m side differential element to the point in 3D space at which the field is being evaluated, and m = 1, 2, 3, 4 indicates a given side of the rectangular filament as shown in Figure 3.1. Note that d~l m is in the direction of positive current flow, which in this case obeys the right hand rule. Evaluation of (3.2) in the rectangular coordinate system of Figure 3.1 results in h(x, y, z, j) = j[gx(x, y, z, a, b, c)ˆi + gy(x, y, z, a, b, c)ˆj + gz(x, y, z, a, b, c)ˆk] = jg(x, y, z, a, b, c), (3.3) where gx = z − c 4 1 '1 y + b '6 − y − b '5 − 1 '2 y + b '8 − y − b '7 , gy = z − c 4 1 '3 x + a '7 − x − a '5 − 1 '4 x + a '8 − x − a '6 , gz = 1 4 x + a '1 y + b '8 − y − b '7 − x − a '2 y + b '6 − y − b '5 + y + b '3 x + a '8 − x − a '6 − y − b '4 x + a '7 − x − a '5 , '1 = (x − a)2 + (z − c)2, '2 = (x + a)2 + (z − c)2, '3 = (y − b)2 + (z − c)2, 33 '4 = (y + b)2 + (z − c)2, '5 = [(x − a)2 + (y − b)2 + (z − c)2]1/2, '6 = [(x − a)2 + (y + b)2 + (z − c)2]1/2, '7 = [(x + a)2 + (y − b)2 + (z − c)2]1/2, '8 = [(x + a)2 + (y + b)2 + (z − c)2]1/2. In addition, (x, y, z) is the point in 3D space at which the field is being evaluated, ˆi , ˆj and ˆk are unit vectors in the x, y, and z directions, a and b indicate the dimensions of the rectangular loop (see Figure 3.1), and c denotes the offset of the loop from the x−y plane of Figure 3.1. It should be noted that there is a fundamental difference between parameters a, b and c. The parameters a and b denote dimensions, hence a, b > 0. However, c denotes an offset of the loop from a plane. In this formulation c > 0 indicates a positive offset, while c < 0 indicates a negative offset. Additionally it may be observed that if the loop resides in the x − y plane as seen in Figure 3.1, c = 0. As can be seen from (3.3), the resultant field from the single rectangular filament is a function of the applied current multiplied by the vector field g. The “g field” is also seen to be a nonlinear function of the 3dimensional location of the point at which the vector field h is being evaluated and the geometry of the filamentary loop. In addition to the vector field h (or g), the calculation of the resultant force on a particle requires additional quantities. It is known in [19] that the vector force on a magnetized sphere approximated as a dipole residing in free space is given as ~FMag(h) = 4 μ0r3 3 M· r h, (3.4) where μ0 is the permeability of free space, r is the radius of the sphere (not to be confused with the magnitude of the vector ~r given in (3.1)) , M is the magnetization vector field of the sphere, and h is the applied magnetic field from the filamentary loop at the location of the particle. It must also be noted that in most fluid media, the spherical object will be unconstrained in the rotational sense resulting in an alignment of the magnetization field with the applied field [32]. Use of this fact along with (3.3) and expansion of the vector quantity (M· r)h in (3.4) results in [33] M h h · r h = jM g " gx @gx @x + gy @gx @y + gz @gx @z ˆi + gx @gy @x + gy @gy @y + gz @gy @z ˆj + gx @gz @x + gy @gz @y + gz @gz @z ˆk # , (3.5) 34 where M is the magnitude of magnetization and g = q g2x + g2 y + g2 z . As can be seen in (3.5), to effectively describe the resultant vector force on the spherical target, not only are the g field components needed, but the partial derivatives are required as well. It should be noted that for a linear paramagnetic sphere, the force relation may be rewritten as ~FMag,PM(h) = 2 μ0r3( m) ( m + 3) r h · h , where m is the magnetic susceptibility of the paramagnetic material. The quantity r(h · h) will then take on the same form as (3.5) with the substitutions of jM/g ! 2j2 and @gm/@n ! @gn/@m, (m, n = x, y, z). The partial derivative terms in (3.5) may be directly evaluated from the g field components given in (3.3). While this has been performed, it proves to be quite tedious and will be omitted here for brevity, however these functions for the rectangular filamentary loop are available in Appendix A. 3.2 Multiturn Expansion To evaluate the resultant vector force on a spherical object from a coil with multiple turns, one may utilize the filamentary model of Section 3.1 and the principle of superposition. If each turn of a multiturn coil is approximated as a single filamentary loop, then the total magnetic field H at a given point in 3D space from a multiturn coil may be given as H = h1 + h2 + · · · + hq = Xq p=1 hp, (3.6) where q is the total number of turns for the coil. Equation (3.6) is only valid with adherence to a few important caveats. First, the individual fields h1, · · · , hq are given in the same coordinate system. For the principle of superposition to hold, the resultant field from each turn must be presented in terms of the same basis. Second, the point in 3D space where the field is being evaluated may not be coincident with a point on a filamentary loop. As can be seen from (3.2), application of the law of BiotSavart for the coincident loop results in an undefined integrand. This is not expected to be an issue due to the fact that no trajectory of the object can be coincident with the coil volume, hence points coincident with the filamentary loops are not needed in this analysis. Application of (3.3) to (3.6) results in H = jg1 + jg2 + · · · + jgq = j Xq p=1 gp = jG, (3.7) where gp is the geometry dependent vector field of the pth loop. Note that in (3.7) the total vector field H is a function of the current j in each turn multiplied by the vector field G which is a function 35 z y x t l d Figure 3.2: Schematic of a rectangular coil with a “perfect wind” packing pattern. The diameter of the wire is denoted as “d”, while the number of turns/layer and number of layers are denoted t and l, respectively. of the number of current loops, the geometry of the current loops, and the location of the loops in 3D space. To formulate the G field expression describing the multiturn coil, use may be made of the expressions in (3.3). As can be seen, (3.3) gives the single loop g field components in the coordinate system of Figure 3.1. This same description may be used for each of the the filamentary loops comprising the approximation of the multiturn coil, provided that the proper values of a, b, and c are used for the loop. Thus to describe the magnetic vector field from a multiturn coil by using the approximate technique of representing the turns as a finite number of current loops, one must only supply the proper values of a, b, and c in equation (3.3) for each pth loop and then apply equation (3.7). If it is assumed that each of the wound coil’s turns may be represented by one concentric loop and that the coil is wound with “perfect” packing, then a schematic of this system may be formed in Figure 3.2. The representation of a wound coil as a grouping of discrete loops is indeed an approximation of the true physical system. In reality, the coil consists of a single conductor wrapped multiple times about a given geometry and the previous “wrappings”. In fact, for a wound coil to be a continuous conductor there must be a region of wire “crossover” which is obviously being neglected in this technique. Additionally, the wound path is in actuality a spiral. This results in an angular 36 component of the current path being present that is not represented by the concentric current loops. This angle can be envisioned as one which would be between the current path and the x − y plane of Figure 3.2. Finally, the ability to fabricate a “rectangular” coil may be quite elusive. While the innermost layers of wire may be very close to a rectangle due to the shape of the bobbin on which the wire is wound, inevitably some rounding of the coil corners will occur. This fact is also being neglected for the analytical model formulation of this study. As is hinted by Figure 3.2, four parameters are necessary to establish the “perfectly” wound, concentric loop model. These consist of the number of wound layers l, the number of turns per layer t, the diameter d of the wire used, and the z direction in which the coil is built. For the representation in Figure 3.2, the coil is built in the positive z direction of the coil’s coordinate system. The parameter t may be thought of as the amount of coil buildup in the z direction while l is the amount of coil buildup in the transverse x − y directions. Note that with the cross sectional wind geometry chosen in Figure 3.2, all layers have the same number of turns per layer. To establish an expression for gp, p = 1, . . . , q, in (3.7), the following relations may be used: G = Xq p=1 gp = Xl =1 Xt =1 g(x, y, z, a , b , c± ), (3.8) where a = a0 + d p 3 2 ( − 1) b = b0 + d p 3 2 ( − 1) c± = 8>>< >>: ±[d( − 1)], for = 1, 3, 5, . . . (odd); ±[ d 2 + d( − 1)], for = 2, 4, 6, . . . (even). In equation (3.8), a0 and b0 denote the initial positions of the innermost layer of the rectangular coil, and the sign of c± indicates the build direction along the ±z axis as indicated in Figure 3.2. Note that q = tl and also that this formulation begins the coil in the local x − y plane, building out from there. An analogous expression may be formed for the partial derivative terms needed in equation (3.5), available here as @Gm @n = Xq p=1 @gm,p @n = Xl =1 Xt =1 @gm(x, y, z, a , b , c± ) @n , (3.9) where m, n = x, y, z. Expressions (3.8) and (3.9) may now be used to calculate the G field and partial derivatives of the G field for a given coil of t turns per layer and l layers (q = tl total turns) at a given point in 3D space in the local coil coordinate system of Figure 3.2. The configuration 37 proposed to facilitate the particle motion consists of four coils, hence an expansion of the single coil model to the quadcoil model is necessary. 3.3 Quad Coil Model The expressions derived in Section 3.2 may also be used to describe the multicoil configuration. A schematic of this system is available in Figure 2.2. Once again using the principle of superposition, the total field resulting from an array of multiturn coils may be given as H = H1 +H2 + · · · +HQ = XQ P=1 HP , where Q represents the number of coils in the array. Note that in system of Figure 2.2, Q = 4. Representing the total field in terms of the G field and current of each coil in the array results in H = XQ P=1 jPGP = j1G1 + · · · + jQGQ = [J]T 2 66664 G1 ... GQ 3 77775 = [J]T [G], where the individual fields GP must be given in a common basis. As has been shown in [32] and Chapter 2, the vector field (M·r)H proportional to the force resulting on a magnetized sphere for an array of coils is derived as M H H· r H = M[J]T H " [GX] @GX @X T + [GY ] @GX @Y T + [GZ] @GX @Z T ˆi + [GX] @GY @X T + [GY ] @GY @Y T + [GZ] @GY @Z T ˆj + [GX] @GZ @X T + [GY ] @GZ @Y T + [GZ] @GZ @Z T ˆk # [J], (3.10) where: H = r [J]T h [GX][GX]T + [GY ][GY ]T + [GZ][GZ]T i [J], [GS] = 2 66664 GS,1 ... GS,Q 3 77775 , @GS @S T = @GS,1 @S · · · @GS,Q @S . As previously stated, the field components and partial derivatives of field components in (3.10) are given in the global coordinate system S = X, Y,Z. To express the fields in the global coordinate 38 system, homogeneous coordinate transformations may be employed on equations (3.8) and (3.9) to transform the calculated field results from the local coil coordinates to the global system [26]. To utilize the representation (3.10) and transform to the proper location in the global coordinate system, the individual local coordinate systems for each coil must be established. The local coil coordinate locations are available in Figure 2.2. It is important to note that the local coordinates given in Figure 2.2 are positioned in the plane of the turns closest to the global origin. Additionally, the correct build directions in the local systems must be chosen, and hence the proper sign of c± in (3.8) and (3.9) selected. For this system it is determined by the orientation of the local coordinate systems; coils one and two utilize c− and coils three and four utilize c+ . As can be seen, local coordinate system 2 is translated and not rotated relative to the global system. Therefore the homogenous coordinate transformation may be given as [A2] = 2 64 [R2] [d2] ; 1 3 75 , (3.11) where [R2] = 2 66664 1 0 0 0 1 0 0 0 1 3 77775 , [d2]T = 0 0 −c0 , the quantity [R2] denotes the rotation component of the transform, [d2] denotes the translation component of the transform, and c0 is the distance from the global origin to the x − y plane of the closest turns of the coil to the global origin. Note that because coil 2 is only translated in relation to the global coordinates, the rotation matrix is effectively an identity matrix. For the other three coils this will not be the case. The coordinate transforms for coils 1, 3, and 4 may be given as [A1] = 2 64 [R1] [d1] ; 1 3 75 , (3.12) [A3] = 2 64 [R3] [d3] ; 1 3 75 , (3.13) and [A4] = 2 64 [R4] [d4] ; 1 3 75 , (3.14) 39 where: [R1] = 2 66664 1 0 0 0 −1 0 0 0 −1 3 77775 , [d1]T = 0 0 c0 = −[d2]T , [R3] = 2 66664 1 0 0 0 0 −1 0 1 0 3 77775 , [d3]T = 0 −c0 0 = −[d4]T , [R4] = 2 66664 1 0 0 0 0 1 0 −1 0 3 77775 . Equation (3.8) may now be used with the relations in (3.11), (3.12), (3.13), and (3.14) to formulate the G field components (given in the global coordinate system) of equation (3.10). This may be accomplished by first representing global point at which the sphere is located in the local basis of each coil, calculating each coil’s local G field and @G/@s contributions (where s = x, y, z), and finally transforming those results back into the global coordinate system for use in (3.10). To transform the point S = (X, Y,Z) to a point sP = (xP , yP , zP ) in the Pth coil’s local coordinate system, it is beneficial to first present the point in its augmented homogeneous representation [26] [ eS]T = X Y Z 1 , [esP ]T = xP yP zP 1 . Transformation of the global point may then be achieved by [esP ] = [AP ]−1[ eS]. (3.15) Once the global coordinate is transformed into the individual coordinate systems, equations (3.3), (3.8), and (3.9) may be used to calculate the local field and spatial partial derivatives of the field. Following this, the local coil field contributions may be transformed back to the global coordinate system by utilizing [GP ] = [RP ][GP ], (3.16) where [GP ]T = GX,P GY,P GZ,P , [GP ]T = Gx,P Gy,P Gz,P . 40 The transform of the spatial partial derivatives of GP is not as elegant as (3.16). This vector field may be transformed back into the global coordinates using [rSGX,1] = [R1][rs,1Gx,1] [rSGY,1] = −[R1][rs,1Gy,1] (3.17) [rSGZ,1] = −[R1][rs,1Gz,1], [rSGX,2] = [R2][rs,1Gx,2] [rSGY,2] = [R2][rs,1Gy,2] (3.18) [rSGZ,2] = [R2][rs,1Gz,2], [rSGX,3] = [R3][rs,1Gx,3] [rSGY,3] = −[R3][rs,1Gy,3] (3.19) [rSGZ,3] = [R3][rs,1Gz,3], and [rSGX,4] = [R4][rs,1Gx,4] [rSGY,4] = [R4][rs,1Gy,4] (3.20) [rSGZ,4] = −[R4][rs,1Gz,4], where rS and rs,P denote the gradient operation in the global coordinates S and local coordinates s, P (P = 1, 2, 3, 4), respectively. Once the G field and partial spatial derivatives of G have been calculated, (3.10) may be utilized to calculate the vector field proportional to the force on the magnetized sphere, provided that the applied current vector [J] is given as well. 3.4 Coil Resistance and Inductance As has been shown in [32] and Chapter 2, the resistance, self inductance, and mutual inductance are needed for an accurate representation of the electrical dynamics of the system shown in Figure (2.2). For a volumetric coil residing in free space with the resulting vector magnetic potential denoted as ~A, the inductance L may be calculated analytically by L = 1 j2 Z V ~A ·~jdV, (3.21) where ~A = 1 4 Z V ~j R dV, 41 ~j is the current density in the coil, V is the volume of the coil, and R is the distance from the volume element dV to the point at which the vector potential is being evaluated [19]. As can be seen, (3.21) is the result of two triple integrals and may prove to be quite challenging. Keeping with the spirit of using a finite loop representation of the field result, the inductance may be formulated in a like manner. Inductance of a coil is defined as L = q j , where is the total magnetic flux linked by the each turn, q is the total number of turns, and j is the current in the coil [31]. Expanding the definition to include the effect between two coils, an examination on a turn by turn basis yields a new definition LM,N = M,N jN = Pq p=1 M,N,p jN , (3.22) where M,N is the total magnetic flux generated by the Nth coil linked by the Mth coil, jN is the current flowing in the Nth coil, and M,N,p is the total magnetic flux generated by the Nth coil linked by the pth filamentary loop of the Mth coil, summed up for q total loops. In this system, the self inductance of a coil may be defined as LM,M and the mutual inductance between two coils as LM,N where M,N = 1, 2, 3, 4 and M 6= N. It is well known that LM,N = LN,M [31, 19, 20]. Through symmetry, it may also be seen that for the quad coil system of Figure 2.2, L1,3 = L1,4 = L2,3 = L2,4 and L1,2 = L3,4. It is obvious that for four identical coils L1,1 = L2,2 = L3,3 = L4,4. From these symmetry observations, it is evident that only three formulations are needed to summarize all of the self and mutual inductance values for the system in Figure 2.2. Calculation of L1,1, L3,1, and L2,1 will facilitate a full description of the system inductance. For the first treatment, the self inductance (M = N) will be examined and the indices M and N will be temporarily dropped. In a general sense, the total flux linked by the pth loop of a coil may be defined as the surface integral p = Z @Sp B · d ~ Sp, (3.23) where B is the total flux density vector field and @Sp is the surface outlined by the pth loop. For the single filamentary loop of Figure 3.1, the total flux linked may be rewritten as = Z @A BzdA, (3.24) where Bz is the component of the total flux density in the z direction linked by the loop and @A is the area outlined by the loop. One may suggest now that (3.24) be solved directly using the relation 42 B = μ0H =) Bz = μ0Hz. For the multiturn coil with q loops, the z component of the total flux density may be formulated from (3.6) and (3.7) as Bz = μ0Hz = jμ0Gz = jμ0 Xq p=1 gz,p. This may be used with (3.24) and (3.8) to formulate the total flux (produced by a coil having q = tl turns) linked by the same coil’s ˜pth loop as ˜p = jμ0 Xl =1 Xt =1 Zbp˜ −b˜p Zap˜ −a˜p gz,p(x, y, z˜p, a , b , c+ )dxdy, (3.25) where z˜p is the z location of the ˜pth loop in the local coordinate system of Figure 3.2. It must also be noted that for the self inductance formulation, a , b , and c+ obey the definitions provided in (3.8). Substitution of the relation for gz from (3.3) into the integrand of (3.25) makes for a demanding exercise in calculus. The evaluation of the total flux can be simplified by utilizing the magnetic vector potential. For the single filamentary loop of Figure 3.1, the magnetic vector potential may be formulated as ~A = j 4 I d~l r = jA, (3.26) where A may be thought of as a geometric vector potential. The magnetic vector potential is related to the magnetic field through the relation H = r×~A, therefore the flux density may be related to the geometric vector potential through B = jμ0(r×A). This may be substituted into (3.23), resulting in a formula describing the total flux linked by the pth loop in terms of total vector potential p = jμ0 Z @Sp (r × A) · d ~ Sp. (3.27) The integration of (3.27) is no less complicated than (3.23), however (3.27) may now be reduced from a surface integral to a path integral by use of Stokes’ theorem [33]. Application of this theorem results in the total flux linkage being defined as p = jμ0 I @sp A · d ~ sp, (3.28) where @sp is the closed path around the pth loop. Clearly, this is a more mathematically tractable problem. Before (3.28) may be used, the quantity A must be evaluated. Initially, the vector potential of the single filamentary loop, denoted by the vector field a, will be examined and then superposed with all loops comprising a coil to form A. The vector field a may be evaluated by using (3.26) and 43 a piecewise integration similar to the technique used in (3.2), a = 1 4 " Z b −b d~l 1 r1 + Z −a a d~l 2 r2 + Z −b b d~l 3 r3 + Z a −a d~l 4 r4 # . (3.29) Note that the vector field a exists only in the ˆi and ˆj directions. Evaluation of (3.29) in the rectangular coordinate system of Figure 3.1 yields a(x, y, z, a, b, c) = ax(x, y, z, a, b, c)ˆi + ay(x, y, z, a, b, c)ˆj , (3.30) where ax = 1 4 ln #2#3 #1#4 , ay = 1 4 ln #6#7 #5#8 , #1 = (x + a) + '7, #2 = (x − a) + '5, #3 = (x + a) + '8, #4 = (x − a) + '6, #5 = (y − b) + '5, #6 = (y + b) + '6, #7 = (y − b) + '7, #8 = (y + b) + '8. In a fashion similar to that of (3.8), the total geometric vector potential for the coil in Figure 3.2 with q turns consisting of l layers and t turns per layer is A = Xq p=1 ap = Xl =1 Xt =1 a(x, y, z, a , b , c+ ). (3.31) Equation (3.31) may be substituted into (3.28) to calculate the total flux linked by the path @sp. The flux linked by the ˜tth turn in the˜l th layer may now be ascertained through the piecewise integration ˜t,˜l = jμ0 Xl =1 Xt =1 " Z b˜l −b˜l ay(a˜l , y, c˜t, a , b , c+ )dy + Z −a˜l a˜l ax(x, b˜l , c˜t, a , b , c+ )dx+ Z −b˜l b˜l ay(−a˜l , y, c˜t, a , b , c+ )dy + Z a˜l −a˜l ax(x,−b˜l , c˜t, a , b , c+ )dx # . (3.32) Solution of (3.32) yields ˜t,˜l = jμ0 4 Xl =1 Xt =1 ˜t,˜l , , , (3.33) 44 where ˜t,˜l , , = b˜l ln Y8 m=5 #m(I1)#m(I3) #m(I2)#m(I4) + a˜l ln Y4 m=1 #m(I3)#m(I4) #m(I1)#m(I2) + b ln Y8 m=5 #m(I1)#m(I4) #m(I2)#m(I3) (−1)m + a ln Y4 m=1 #m(I1)#m(I4) #m(I2)#m(I3) (−1)m , I1 = (a˜l , b˜l , c+ ˜t , a , b , c+ ), I2 = (−a˜l , b˜l , c+ ˜t , a , b , c+ ), I3 = (a˜l ,−b˜l , c+ ˜t , a , b , c+ ), I4 = (−a˜l ,−b˜l , c+ ˜t , a , b , c+ ). Note that in (3.33), I denotes the variable assignments in the functions # and '. Additionally, a˜l = a , b˜l = b , and c+ ˜t = c+ , for ˜l = , and ˜t = . Use of (3.32) and (3.22) allows for the formulation of the self inductance L1,1 as L1,1 = 1 j ˜l X ˜ =1 Xt˜ ˜ =1 ˜ ˜ = μ0 4 ˜l X ˜ =1 Xt˜ ˜ =1 Xl =1 Xt =1 ˜ ,˜ , , . While (3.33) is convenient, it does have limitations. The integrands of (3.32) are undefined when ˜t = and ˜l = , and hence must be ignored when calculating (3.33), or in other words e ˜t,˜l = 8>>< >>: 0, for ˜t = and ˜l = ; jμ0 4 Pl =1 Pt =1 ˜t,˜l , , , for ˜t 6= or ˜l 6= . (3.34) A modified self inductance may now be formed which neglects the undefined integrands, eL 1,1 = 1 j ˜l X ˜ =1 Xt˜ ˜ =1 e ˜ ˜ . (3.35) This will result in error being introduced to the calculation, however if the coil windings are assumed to be closely packed, an estimate of the error may be obtained. Assuming that the coil windings are closely packed allows for the approximations A = a1 + a2 + · · · + aq qa, p = jμ0 I @sp A · d ~ sp qjμ0 I @s a · d~s. The total flux linked by the coil would approximately be = Xq p=1 p q2jμ0 I @s a · d~s. 45 If the indefinite path integral is neglected, then eA (q − 1)a, e p (q − 1)jμ0 I @s a · d~s. The total flux linked by the coil when neglecting the undefined integrands is then e = Xq p=1 e p q(q − 1)jμ0 I @s a · d~s. Using (3.22), the error estimate Es of the self inductance resulting from neglecting the undefined integrands may be formed as Es L −eL L = 1 q . As can be seen, the error in the self inductance calculation utilizing finite loops and neglecting the undefined integrands is proportional to the inverse of the number of turns in the coil. A corrected self inductance calculation may then be formed as bL 1,1 = eL 1,1 1 + 1 q . To calculate the mutual inductances L3,1 and L2,1, a similar approach may be used. Additionally, for the mutual inductance formulations it will not be necessary to neglect some of the integrations, as there will no longer be any overlapping integration paths and hence no undefined integrands. Performing a similar derivation to that of the self inductance formulation yields the mutual inductance of two adjacent coils as L3,1 = μ0 4 ˜l X ˜ =1 Xt˜ ˜ =1 Xl =1 Xt =1 ˜ ,˜ , , , (3.36) where ˜ ,˜ , , = (a˜ )3,1 ln Y4 m=1 #m(L3)#m(L4) #m(L1)#m(L2) + a ln Y4 m=1 #m(L1)#m(L4) #m(L2)#m(L3) (−1)m + X8 m=5 ['m(L4) − 'm(L3) + 'm(L1) − 'm(L2)], L1 = ((a˜ )3,1, (b˜ )3,1, (c˜ )− 3,1, a , b , c− ), L2 = (−(a˜ )3,1, (b˜ )3,1, (c˜ )− 3,1, a , b , c− ), L3 = ((a˜ )3,1,−(b˜ )3,1, (c˜ )+3 ,1, a , b , c− ), L4 = (−(a˜ )3,1,−(b˜ )3,1, (c˜ )+3 ,1, a , b , c− ), (a˜ )3,1 = a0 + d p 3 2 (˜ − 1), 46 (b˜ )3,1 = 8>>< >>: c0 + (˜ − 1)d; if ˜ is odd, c0 + d 2 + (˜ − 1)d; if ˜ is even, (c˜ )± 3,1 = c0 ± b0 + d p 3 2 (˜ − 1) , and the mutual inductance of two opposite coils as L2,1 = μ0 4 ˜l X ˜ =1 Xt˜ ˜ =1 Xl =1 Xt =1 ˜ ,˜ , , , (3.37) where ˜ ,˜ , , = b˜ ln Y8 m=5 #m(J2)#m(J4) #m(J1)#m(J3) + a˜ ln Y4 m=1 #m(J1)#m(J2) #m(J3)#m(J4) + b ln Y8 m=5 #m(J2)#m(J3) #m(J1)#m(J4) (−1)m + a ln Y4 m=1 #m(J2)#m(J3) #m(J1)#m(J4) (−1)m , J1 = (a˜ , b˜ , (c˜ )2,1, a , b , c− ), J2 = (−a˜ , b˜ , (c˜ )2,1, a , b , c− ), J3 = (a˜ ,−b˜ , (c˜ )2,1, a , b , c− ), J4 = (−a˜ −, b˜ , (c˜ )2,1, a , b , c− ), (c˜ )2,1 = 8>>< >>: 2c0 + [d(˜ − 1)], for ˜ = 1, 3, 5, . . . (odd); 2c0 + [ d 2 + d(˜ − 1)], for ˜ = 2, 4, 6, . . . (even). The coil resistance calculation is quite straightforward. If the total length of the wire in the coil is known, then this length may be multiplied by the resistance per length Rlength of the wire used to wind the coil. As in the previous field calculations and inductance formulations, the coil is approximated as a series of finite loops. This results in the calculation of resistance as R = 4tlRlength a0 + b0 + d p 3 2 (l − 1) . (3.38) Alternatively, coil resistance may be calculated on a weight W basis, where the parameter of resistance per unit weight Rwgt is utilized to calculate the resistance R = RwgtW. While (3.38) is a useful formulation for use in a predictive model or optimization routine, the weight based calculation is much more practical when winding the coil. 47 CHAPTER 4 Design and Optimization of Experiment 4.1 Functional for Minimization and Algorithm In a system such as the one proposed in Figure 2.2, there are many design parameters which must be determined before fabrication of an experimental setup can proceed. It may be assumed for the purposes of this study that certain aspects of the system are “fixed” and are not considered variables to be changed or adjusted in any kind of optimization procedure. The fixed characteristics are listed as follows: sphere radius, sphere magnetization, fluid viscosity, and sphere 2D trajectory. Additionally, a large amount of 12 gauge heavy insulation magnet wire was already in hand, resulting in wire diameter being fixed as well. In contrast, certain parameters of this system will be used as variables in an optimization process. These parameters are given in Table 4.1. Upon first observation, one may reason that the coil resistance should be among the variables listed in Table 4.1, however this is not the case. Fabrication of the experimental apparatus will utilize a magnet wire already in hand (12 gauge), thus establishing the wire diameter. The coil resistance is a function of wire diameter, electrical conductivity of the wire utilized, and length of the wire wound into coil form. Of these three criteria, only the wire length will be varied during the optimization and it is a function of the first four parameters of Table 4.1. In devising a “bench top” test apparatus, it is desirable to have a hardware setup which will move the sphere with the least amount of effort. Hence a minimization of some sort is in order. From the standpoint of the physical system, the attribute which must be minimized is the power applied to the electromagnet coils. The instantaneous power supplied to any reactive electrical component is defined as P = vj, where v is the potential across the device and j is current flowing in the device [34]. Adopting the approach used in [32] and representing the coil as an RL circuit, the instantaneous power in the mth coil may be stated as Pm = Rmj2m + jm XQ n=1 Lm,n djn dt , (4.1) 48 Table 4.1: System Parameters to be Used for Optimization parameter symbol unit number of wire turns/layer t    number of wire layers l    coil inner “x” dimension ain inch coil inner “y” dimension bin inch distance of coil face from global origin c0 inch where Rm and jm are the resistance and current for the mth coil, Lm,n is the mutual (m = n) or self (m 6= n) inductance of the mth coil on the nth coil, and Q is the total number of coils. Consistent with Chapter 2, the optimization problem will be based the ability of the hardware system to move a sphere about a circular orbit. Hence, symmetry may be leveraged and only coils 1 and 3 examined. To calculate the instantaneous power of the coil combination depicted in Figure 2.10, denoted here as 1,3 = (P1,P3), the formulations for coil resistance and inductance (both self and mutual) developed in Chapter 3 may be employed. The analytical evaluation of the time derivative of current would prove to be quite tedious, if even possible, due to the required formulation of the partial derivatives of the vector field term ~kFM. A complete mathematical description of the vector field term is available in Chapter 2. Because of this, the time derivative of the required currents in coils 1 and 3 will be calculated numerically from discrete points comprising the inverse current solution of a given system configuration. Obviously in a reactive system such as an RL ciruit, achieving the instantaneous current change at the beginning and end of the inverse current solution as depicted in Figure 2.10 would be impossible. Indeed an inductive “switching effect” would be inherent in the proposed technique, however for the purposes of optimizing an experimental design this will be neglected. For the set of discrete points making up the inverse current solution J ,1,3, where = a, b, c, d, a forward difference approximation was used to calculate the derivative at the initial point while a backward difference approximation was used to calculate the derivative at the last point. For all other points making up the set, a central difference technique was used. This will provide a reasonable approximation of the function J˙ ,1,3, provided that the time step spacing is small [35]. The time derivative of current was calculated using points spaced at 10 msec intervals. An 49 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ a /dt (A/sec) dj 1 /dt dj 3 /dt 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ b /dt (A/sec) 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ c /dt (A/sec) 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ d /dt (A/sec) time (seconds) Figure 4.1: Calculated current derivative from inverse solution of coil combination 1&3. Displayed results correspond to system with R1,3 = 0.547 , L1,1 = 6.79mH, and L1,3 = 0.58mH (calculated). example calculation on the current trajectory depicted in Figure 2.10 is shown in Figure 4.1. It should be noted that the trace shown in Figure 2.10 corresponds to a system where the resistance, inductance and mutual inductance have been calculated to be R1,3 = 0.547 , L1,1 = 6.79mH, and L1,3 = 0.58mH. Utilizing equation (4.1) with the discrete traces of j1,3 and dj1,3/dt along with the calculated values for resistance and inductance allows for the estimation of the instantaneous power required by the coil combination which results in the inverse current solution. The instantaneous power calculation for this configuration is shown in Figure 4.2. As can be seen, the instantaneous power has two distinct forms rather than four. This is indicative of the same power being needed for the solutions corresponding to opposite polarities. With the proposed minimal coil set technique, two coils will be powered at any one given time. Hence, it would be advantageous to utilize an expression as the functional for minimization which incorporates the power applied to both members of the coil set. The Euclidean norm was employed to satisfy this convenience and is defined here as k m,nk = p P2m + P2n , where m, n = 1, 2, 3, 4. Obvioulsy k m,nk = k n,mk and for the case being examined m = 1 and n = 3. The norms of the two instantaneous power traces are shown in Figure 4.3. From 50 0 1 2 3 4 5 6 7 8 9 10 0 2 4 !a (Watts) P 1 P 3 0 1 2 3 4 5 6 7 8 9 10 0 2 4 !b (Watts) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 !c (Watts) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 ! d (Watts) time (seconds) Figure 4.2: Instantaneous power required for inverse current solution of coil combination 1&3. 0 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 time (seconds) !1,3  (Watts) set "a" set "c" Figure 4.3: Plot of k 1,3k for solution sets a and c. The norm of solution sets b and d are identical to sets a and c, respectively. 51 the power norm, a single metric may be defined to describe the performance of a coil system with the parameters seen in Table 4.1. If the maximum value of the current norm, defined here as M= max k m,nk, is used as criteria for minimization, then the minimized configuration Smin would be the set of parameters S = {t, l, a0, b0, c0} which results in the minimum M when compared to all parameter values possible. More precisely, {t, l, ain, bin, c0} = S = Smin if 9{t, l, ain, bin, c0} : k m,n(t, l, ain, bin, c0)k inf{M(S)}8S 2 R. In an effort to reduce the number of parameters utilized in the search routine, it was decided to make c0, the distance from the face of the coils to the origin (see Figure 2.2), a function of the parameters l and bin. To achieve this, a spacing criteria csp is defined as csp = bin + d 2 + (l − 1) p 3 2 + sep, (4.2) where d is the wound wire diameter, bin is the inner Z dimension of the rectangular coil, sep is a separation constant, and l is the number of layers of wire wound on the coil. It should be noted that the formulation of (4.2) is based upon a coil having a “perfect” winding pattern as defined in Chapter 3. The minimum coil to coil spacing for this experiment was chosen to be 5”, resulting in a minimum coil spacing from the origin of cmin = 2.5”. In addition, the separation constant in (4.2) was chosen to be sep = 0.25”. To reduce the number of parameters for the optimization problem, c0 may then be defined as c0 = 8>>< >>: cmin if csp cmin, csp otherwise. This relationship sets the coil spacing as 2.5 inches unless the coil physical dimensions do not allow for them to be positioned this closely. In that case, the coil spacing from the origin is dictated by (4.2). The technique selected to find the minimal parameter set is a basic gradient method: steepest descent [36, 37]. More specifically, a discrete steepest descent algorithm was used because the system model was evaluated at specific points in the parameter space. This work will utilize a rather straigtforward version of this gradient method, which will be presented in two parts. The first component may be summarized as a method by which, for a given ain and bin, the number of turns per layer t and the number of layers l resulting in the minimumMis found. To achieve this, we may look at M as a function of the minimizing variables t, l, ain, and bin, hence M = M(t, l, ain, bin). For a given initial guess (ti, li) and a given (ain, bin), the initial value of M = M(ti, li, ain, bin) is calculated. The points “around” the initial guess are utilized to then calculate the values of M around the initial guess. These points make up the ordered set Si u = {(ti+1, li), (ti+1, li−1), (ti, li− 52 yes no Figure 4.4: A flowchart illustrating the algorithm used to find the number of turns tf and number of layers lf which results in minimal power used for a given ain and bin. 1), (ti − 1, li − 1), (ti − 1, li), (ti − 1, li + 1), (ti, li + 1), (ti + 1, li + 1)}, where u = 1, 2, . . . , 8 denotes a particular ordered set member. This set may then be used to calculate the values of M at each set point, {M(Si u, ain, bin)} = {M(Si 1, ain, bin), . . . ,M(Si 8, ain, bin)}. The “next step” (ti+1, li+1) may then be determined as the vth element of Si u which results in G = min{M(Si u, ain, bin) − M(ti, li, ain, bin)} and G 0. This process may be repeated until a point (tf , lf ) is found such that G > 0. At which time, a minimizing t and l have been found. This process is illustrated in the flow chart shown in Figure 4.4. The second component of the complete minimization algorithm may be viewed as an analogy of the first component discussed in the previous paragraph and shown in Figure 4.4. It may be summarized as a method by which, for a given t and l, the inner x coil dimension ain and the inner y coil dimension bin resulting in the minimumMis found. As before, this may be achieved by regarding M as a function of the variables which are being “varied”. For a given initial guess (ain,j , bin,j) and a known value (t, l), the initial value ofM=M(t, l, ain,j , bin,j) is calculated. Once again, the points “around” the initial guess are utilized to calculate the values of M around the initial guess. These points make up the ordered set Tj u = {(ain,j + , bin,j), (ain,j + , bin,j − ), (ain,j , bin,j − ), (ain,j − , bin,j − ), (ain,j − , bin,j), (ain,j − , bin,j + ), (ain,j , bin,j + ), (ain,j + , bin,j + )}, where u = 53 yes no Figure 4.5: A flowchart illustrating the algorithm used to find the inner x coil dimension ain,f and inner y coil dimension bin,f which results in minimal power used for a given t, l, and . 1, 2, . . . , 8 denotes a particular ordered set member and denotes a search length. This set may then be used to calculate the values ofMat each set point, {M(t, l,Tj u)} = {M(t, l,Tj 1), . . . ,M(t, l,Tj 8)}. The “next step” (ain,j+1, bin,j+1) may then be determined as the wth element of Tj u which results in G = min{M(t, l,Tj u) −M(t, l, ain,j , bin,j)} and G 0. This process may be repeated until a point (ain,f , bin,f ) is found such that G > 0. At which time, a minimizing ain and bin have been found. This process is illustrated in the flow chart shown in Figure 4.5. As can be seen, Figures 4.4 and 4.5 are similar in structure. To achieve a single algorithm which finds a minimizing point (tf , lf , ain,f , bin,f ), one algorithm will be “nested” in the other. For this study, it was decided to nest the algorithm depicted in Figure 4.4 into the algorithm depicted in Figure 4.5. The flow of the complete process is similar to that of the procedures previously outlined. For an initial guess (ti, li, ain,j , bin,j ), the initial value of M = M(ti, li, ain,j , bin,j) must be calculated. Subsequently, M = M(t , l , ain,j , bin,j) which has been minimized on (t, l) may be calculated using a nested algorithm similar to that in Figure 4.4. In this case, represents an arbitrary number of iterations required for convergence. Similarly, the minimal {(t , l )u} for each set member of Tj u may be calculated as well by a nested (t, l) minimization algorithm, where represents an arbitrary number of iterations required for convergence. These results are then used to check 54 yes no Figure 4.6: A flowchart illustrating the complete algorithm used to find the parameters (tf , lf , ain,f , bin,f ) resulting in minimal power. The nested algorithms are indicated by “broken line” boxes. for a descent direction in the variable space. This procedure is then repeated until a minimizing configuration (tf , lf , ain,f , bin,f ) is found. The complete algorithm is shown in Figure 4.6. It must be noted that while not explicitly depicted in Figure 4.6, a value for the search length must be selected before using the algorithm. The final aspect to be resolved before the minimization routine may be implemented is which norm set will be used (see Figure 4.3). For this study, “a” was selected as the set utilized in the functional evaluation. With the minimization routine developed, and the system model in place, a minimizing parameter set may be calculated. 4.2 Minimization Results and Experimental Setup A minimizing parameter set (tf , lf , ain,f , bin,f ) was indeed found from utilizing the routines developed in Section 4.1. Three different search efforts were performed with search lengths valued as = 0.1”, 0.01”, and 0.001”. This method may be viewed as one in which an initial “coarse” search is performed to establish a general area of the extremum and then each successive search utilizes a finer search grid. The optimized design parameters for the quad coil assembly are given in Table 4.2. As can be seen, the optimal coil design corresponds to one in which the coil faces are placed as close 55 Table 4.2: Optimized System Parameters parameter symbol value unit number of wire turns/layer t 18    number of wire layers l 21    coil inner “x” dimension ain 0.439 inch coil inner “y” dimension bin 0.757 inch distance of coil face from global origin c0 2.501 inch to the origin as allowed, in this case 2.5”. This is indicative of the proximity of the particle to the electromagnet as being a major contributor to the amount of power required for particle movement. A coil assembly along with the tray providing the plane of motion for the sphere was constructed to these specifications and is shown in Figure 4.10. As expected, the geometry of the wound coi
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Modeling and Control of a Permanently Magnetized Sphere's Motion by Field Manipulation 
Date  20061201 
Author  Duvall, Matthew Glenn 
Document Type  
Full Text Type  Open Access 
Abstract  A model describing the motion of a magnetized sphere in a fluid medium subject to an applied field is developed. The inspiration for this system is rooted in technologies and procedures being explored in the biomedical community. In these systems, the accepted approximation of Stokes' flow about the sphere is discussed and limitations to this are identified with drag force relations accounting for wall effects presented. For this system, the sphere is approximated as a point dipole and the resulting electromagnetic force model is presented. For motion of the sphere using static electromagnets, a reduced coil set technique is proposed. The technique consists of utilizing the minimum number of coils for generation of the desired force on the sphere. This approach is in contrast to the existing solutions which use the full coil set for parallel and antiparallel motion and single out a solution for the underdetermined system. A technique for determining the proper coil combination is developed from examining the definiteness of the geometric field functions. A coil array configuration consisting of a four coil assembly is investigated to facilitate planar motion of a magnetized spherical particle in a fluid medium. For this configuration, the reduced coil set consists of only two adjacent coils being active at any given time. The exact inverse current solution for the minimized coil set is derived and presented. The system dynamic model is formulated and represented in a nonlinear state space system. This study also derives analytical expressions describing the magnetic field and field gradient terms generated by an assembly of four coils grouped about a coordinate origin with the central axes of the coils all residing in the same plane. This model will prove to be quite useful for optimal design of an apparatus to serve as an experimental platform for such a system. The model predicts the magnetic field produced by an elliptically shaped coil, approximated as a finite number of rectangularly shaped loops. The law of BiotSavart is utilized to calculate the magnetic field and field gradient terms for each loop. The loops are then used to approximate each turn of a wound coil and are superposed to calculate the total field generated by a single coil. Homogenous coordinate transforms are utilized to describe the field generated by the quad coil configuration. An analytical formulation based on finite loops for calculation of self inductance, mutual inductance, and resistance is also derived and presented. Magnetic field measurements of wound coils are performed and compared to the predictive model, with satisfactory results being demonstrated. The technique of representing the wound coils as an assemblage of finite loops is found sufficient to formulate a model which may be used for predictive modeling of a system to facilitate two dimensional motion of a magnetized sphere. A minimization routine based on discrete steepest descent gradient method is developed and presented. The functional used for the minimization routine has been developed and is established as the norm of instantaneous power applied to an adjacent set of coils. This allows for summary of the power in two coils by one calculation. The routine calculates the required coil spacing and coil geometry necessary to achieve motion of the magnetized sphere along a given circular trajectory with minimal power dissipation. Control of a small magnetized spherical particle in 2D space using a square arrangement of electromagnets is presented. A classical controller with and without feedforward is implemented in concert with a novel switching approach. A sample of the experimental results for motion of the particle along a desired trajectory in a fluid medium are presented and discussed. Further, the experimental results are compared with computer simulation results to get insights into the accuracy of the dynamic model. 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  MODELING AND CONTROL OF A PERMANENTLY MAGNETIZED SPHERE’S MOTION BY FIELD MANIPULATION By MATTHEW GLENN DUVALL Bachelor of Science Mechanical Engineering Oklahoma State University Stillwater, Oklahoma, USA 1994 Master of Science Mechanical Engineering Oklahoma State University Stillwater, Oklahoma, USA 1997 Submitted to the Faculty of the Graduate College of Oklahoma State University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy December, 2006 COPYRIGHT c By MATTHEW GLENN DUVALL December, 2006 MODELING AND CONTROL OF A PERMANENTLY MAGNETIZED SPHERE’S MOTION BY FIELD MANIPULATION Dissertation Approved: Eduardo A. Misawa Dissertation Advisor Prabhakar R. Pagilla Frank W. Chambers Martin Hagan Dean of the Graduate College iii ACKNOWLEDGMENTS The author would like to sincerely thank the following individuals: Firstly, my academic advisor Dr. Eduardo A. Misawa whose support and advice greatly facilitated this project. Secondly, my academic committee, specifically Dr. Prabhakar R. Pagilla. His guidance, mentorship, and not to mention the many hours of time sacrificed to examining draft manuscripts proved to be invaluable to my educational experience. I would also like to thank Dr. Gary Young for his support during my teaching assistantship which provided a crucial aspect of the doctoral candidacy experience. In addition I would like to thank Dr. Kenneth L. Pottebaum of Seagate Technology, LLC for his professional guidance and counseling regarding the balance of school and academia. Moreover, I would like to thank Dr. Ryan T. Ratliff of the Boeing Company for his assistance in, and clarification of, various technical difficulties. Furthermore, I would like to thank Leon E. Boomershine andWenjie Li. Without assistance from these two individuals it is doubtful that the experimental hardware setup for this project would have been completed in a timely fashion. Finally, I must offer an extreme amount of gratitude to my wife, Heather D. Duvall, for her unflinching support in my academic endeavors and the many sacrifices for which she has had to endure. Without the support of these individuals, it is inconceivable that this body of work would be completed and this particular scholastic goal accomplished. iv TABLE OF CONTENTS Chapter Page 1 Introductory Work 1 1.1 Impetus for the Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Overview of Relevant Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Contributions of this Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Overview of this Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Dynamic System Model 8 2.1 Dynamic Model Development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Inverse Current Solution for Constrained 2D Motion with a Reduced Coil Set . . . 14 2.3 Current Trajectory Smoothing Technique . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 State Space Model and Closed Loop Control Considerations . . . . . . . . . . . . . . 26 3 Electromagnetic Field Description 32 3.1 The Single Filamentary Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2 Multiturn Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Quad Coil Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.4 Coil Resistance and Inductance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 Design and Optimization of Experiment 48 4.1 Functional for Minimization and Algorithm . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Minimization Results and Experimental Setup . . . . . . . . . . . . . . . . . . . . . 55 4.3 Experimental Platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.4 Experimental Validation of Platform Design . . . . . . . . . . . . . . . . . . . . . . . 59 5 Closed Loop Control of Sphere Motion 65 5.1 Coil Switching Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2 Closed Loop Control Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.3 Closed Loop Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 v 6 Conclusions and Future Work 89 6.1 Conclusions from this Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.2 Recommended Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 BIBLIOGRAPHY 95 A Partial Spatial Derivatives of g Field Components Resulting From a Single Rectangular Filamentary Loop 99 B Power Circuit Diagram 103 C Closed Loop Matlab/Simulink Simulation Source Code 104 C.1 Model Parameters Script Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 C.1.1 Coil Resistance and Weight Calculation Script . . . . . . . . . . . . . . . . . 108 C.1.2 Coil Self Inductance Calculation MEX Script . . . . . . . . . . . . . . . . . . 109 C.1.3 Coil Mutual Inductance MEX Script . . . . . . . . . . . . . . . . . . . . . . . 116 C.1.4 G Field Calculation MEX Script . . . . . . . . . . . . . . . . . . . . . . . . . 121 C.2 Simulink Top level Model and Subsystems . . . . . . . . . . . . . . . . . . . . . . . . 134 C.3 Embedded Matlab Function Block Scripts . . . . . . . . . . . . . . . . . . . . . . . . 147 C.3.1 Square path embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . 147 C.3.2 Point path embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . . 148 C.3.3 Inverse current calculation embedded function. . . . . . . . . . . . . . . . . . 149 C.3.4 Current selector embedded function. . . . . . . . . . . . . . . . . . . . . . . . 209 C.3.5 K matrices embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . . 212 C.3.6 Wall forces embedded function. . . . . . . . . . . . . . . . . . . . . . . . . . . 217 vi LIST OF FIGURES Figure Page 2.1 Typical M − H curve for a ferromagnetic material. . . . . . . . . . . . . . . . . . . . 10 2.2 Schematic of a four coil arrangement to provide forces for planar motion of a sphere. The individual coil coordinates indicate direction of positive current polarity (right hand rule). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Graphical representation of the applied planar magnetic force. . . . . . . . . . . . . . 15 2.4 Solution loci for the 1&4 coil combination in the J space. . . . . . . . . . . . . . . . 20 2.5 Initial current solution for coil combination 1 and 3 during circular trajectory. . . . . 22 2.6 Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces. . . . . . . . . . . . . . . . . . . . . . 22 2.7 Correct current solution for coil combination 1 and 3 during circular trajectory. . . . 23 2.8 Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces are now absent. . . . . . . . . . . . . 24 2.9 A plot of sgn(D1,3), where D = det[eQ ]. A value of +1 indicates that transformation is a rotation, while a value of −1 is indicative of a reflection. . . . . . . . . . . . . . . 24 2.10 Smoothed current solution for coil combination 1 and 3 during circular trajectory. . 26 2.11 Initial and smoothed current solution for coil combination 2 and 3 during circular trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.12 Initial and smoothed current solution for coil combination 2 and 4 during circular trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.13 Initial and smoothed current solution for coil combination 1 and 4 during circular trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1 Schematic of a single filamentary conducting loop. The local coordinate system is fixed to the center of the loop with the central axis aligned with the z coordinate direction. Positive current flow obeys the right hand rule. The numbers indicate sections of conductor while the letters denote the dimensions of the rectangle. . . . . 33 vii 3.2 Schematic of a rectangular coil with a “perfect wind” packing pattern. The diameter of the wire is denoted as “d”, while the number of turns/layer and number of layers are denoted t and l, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.1 Calculated current derivative from inverse solution of coil combination 1&3. Displayed results correspond to system with R1,3 = 0.547 , L1,1 = 6.79mH, and L1,3 = 0.58mH (calculated). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Instantaneous power required for inverse current solution of coil combination 1&3. . 51 4.3 Plot of k 1,3k for solution sets a and c. The norm of solution sets b and d are identical to sets a and c, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 A flowchart illustrating the algorithm used to find the number of turns tf and number of layers lf which results in minimal power used for a given ain and bin. . . . . . . . 53 4.5 A flowchart illustrating the algorithm used to find the inner x coil dimension ain,f and inner y coil dimension bin,f which results in minimal power used for a given t, l, and . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.6 A flowchart illustrating the complete algorithm used to find the parameters (tf , lf , ain,f , bin,f ) resulting in minimal power. The nested algorithms are indicated by “broken line” boxes. 55 4.7 Experimental apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.8 Schematic of experimental setup and signal flow. . . . . . . . . . . . . . . . . . . . . 57 4.9 Schematic of tray cross section. The sphere resides below the centerline of the fluid gap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.10 Experimental apparatus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.11 Experimental setup for measurement of the flux density produced by a single coil in the y direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.12 Experimental setup for measurement of the flux density produced by a single coil in the z direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.13 A comparison of the predictive model for the field component Gy to measured data points. The brackets around the points indicate ±3 distributions of the measurement technique at that point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.14 A comparison of the predictive model for the field component Gz to measured data points. The brackets around the points indicate ±3 distributions of the measurement technique at that point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.1 Calculated current trajectories for different coil combinations. . . . . . . . . . . . . . 67 viii 5.2 Block diagram of closed loop control simulation with a classical “proportional” controller is implemented. The upper diagram is the closed loop system while the lower diagram is a detailed view of the controller. . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 Block diagram of closed loop control simulation with a feedforward block. A classical “proportional” controller is implemented in concert with a feedforward block which estimates the desired forces on the spherical particle. The upper diagram is the closed loop system while the lower diagram is a detailed view of the controller. . . . . . . . 71 5.4 Simulation of closed loop PD controller resultant trajectory following a spiral reference. 72 5.5 Simulation of closed loop PD controller resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.6 Simulation of closed loop feedforward controller resultant trajectory following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.7 Simulation of closed loop feedforward controller resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.8 Experimental closed loop controller resultant trajectory following a spiral reference. . 75 5.9 Experimental closed loop controller resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.10 Experimental closed loop controller with feedforward block resultant trajectory following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 5.11 Experimental closed loop controller with feed forward block resultant current and position error following a spiral reference. . . . . . . . . . . . . . . . . . . . . . . . . 77 5.12 Experimental closed loop controller resultant trajectory following a square reference. 77 5.13 Experimental closed loop controller resultant current and position error following a square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.14 Experimental closed loop controller with feedforward resultant trajectory following a square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.15 Experimental closed loop controller with feedforward resultant current and position error following a square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.16 Experimental closed loop controller resultant trajectory attempting a point regulation. 80 5.17 Experimental closed loop controller resultant current and position error attempting a point regulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.18 Experimental closed loop controller resultant trajectory attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 ix 5.19 Experimental closed loop controller resultant current and position error attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.20 Experimental closed loop feedforward controller resultant trajectory attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.21 Experimental closed loop feedforward controller resultant current and position error attempting to track the fast spiral trajectory. . . . . . . . . . . . . . . . . . . . . . . 83 5.22 Experimental closed loop feedforward controller resultant trajectory attempting to track the fast spiral trajectory. (kP = 0.5) . . . . . . . . . . . . . . . . . . . . . . . . 84 5.23 Experimental closed loop feedforward controller resultant current and position error attempting to track the fast spiral trajectory. (kP = 0.5) . . . . . . . . . . . . . . . . 84 5.24 Experimental closed loop controller resultant trajectory following the fast square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.25 Experimental closed loop controller resultant current and position error following the fast square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.26 Experimental closed loop feedforward controller resultant trajectory following the fast square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.27 Experimental closed loop feedforward controller resultant current and position error following the square reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.28 Experimental open loop controller resultant trajectory attempting to follow a spiral. 87 5.29 Experimental open loop controller resultant current and position error attempting to follow a spiral. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 B.1 Circuit diagram for coil power op amps and signal conditional op amps. Note that these are inverting amplifiers and that each coil has a power circuit. . . . . . . . . . 103 C.1 Top level Simulink Model (“closed loop”). . . . . . . . . . . . . . . . . . . . . . . . . 134 C.2 Spiral trajectory subsystem (closed loop!spiral trajectory). . . . . . . . . . . . . . . 135 C.3 Spiral trajectory path definition (closed loop!spiral trajectory!spiral). . . . . . . . 136 C.4 Square trajectory subsystem (closed loop!square trajectory). . . . . . . . . . . . . . 137 C.5 Point regulator subsystem (closed loop!point regulator). . . . . . . . . . . . . . . . 138 C.6 Controller subsystem (closed loop!Controller). . . . . . . . . . . . . . . . . . . . . . 139 C.7 Inverse current calculation (closed loop!Controller!inverse current calculation). . . 140 C.8 Inverse voltage conversion (closed loop!Controller!inverse voltage conversion). . . 141 C.9 Required force subsystem (closed loop!Controller!required force). . . . . . . . . . 142 x C.10 Plant model subsystem (closed loop!plant model). . . . . . . . . . . . . . . . . . . 143 C.11 Electrical dynamics subsystem (closed loop!plant model!electrical dynamics). . . 144 C.12 Magnetic force subsystem (closed loop!plant model!magnetic force). . . . . . . . . 145 C.13 Sphere dynamics subsystem (closed loop!plant model!sphere dynamics). . . . . . 146 xi LIST OF TABLES Table Page 2.1 Correct Line Loci Groupings for e J Space Solution Coordinates . . . . . . . . . . . . 19 2.2 Modified Sign Groupings for Smooth e J1,3 = ˜j 1 1 1 ,˜j 3 3 3 Solution Coordinates . . 25 2.3 Modified Sign Groupings for Smooth e J2,3 = ˜j 2 2 2 ,˜j 3 3 3 Solution Coordinates . . 28 2.4 Modified Sign Groupings for Smooth e J2,4 = ˜j 2 2 2 ,˜j 4 4 4 Solution Coordinates . . 29 2.5 Modified Sign Groupings for Smooth e J1,4 = ˜j 1 1 1 ,˜j 4 4 4 Solution Coordinates . . 29 4.1 System Parameters to be Used for Optimization . . . . . . . . . . . . . . . . . . . . 49 4.2 Optimized System Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3 Targeted wound coil parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 xii CHAPTER 1 Introductory Work 1.1 Impetus for the Research In recent decades many new technologies have been proposed and investigated in the biomedical sciences. Several studies have suggested the use of magnetic nanoparticles in vivo to the human body. The benefits of such methods are believed to be invaluable because of the minimally invasive capabilities of anticipated procedures. Envisioned applications of nanoparticles include, but are not limited to, cell labeling, magnetic separation, targeted drug delivery, hyperthermic cell treatment, and MRI contrast enhancement [1]. Significant study has been performed regarding these concepts and is available in references [2, 3, 4, 5]. A second novel application of a magnetically reactive device being introduced to the internal workings of a living body and being manipulated by an applied ex vivo magnetic field is magnetic implant and catheter guidance. The practice of this procedure on the human brain is referred to as stereotactic neurosurgery and a significant effort has been put forth on this technique [6, 7]. Similarly, a process utilizing nano sized particles and applied external magnetic fields for intraocular retinal repair is being explored as well [8]. Finally, the fascinating and emerging science of nanoscaled machines has offered the potential employment of microdevices in medical procedures. To reduce device complexity, an external magnetic propulsion technique is envisioned [9, 10]. The common attribute present in all of the aforementioned systems is the use of an applied magnetic field. This field is utilized to exert a desired force and/or torque on an object susceptible to such fields. Upon that realization, one may attempt to first recognize the analogy and similarities to the classic magnetic levitation problem. Indeed a parallel may be observed that manipulation of a particle or microdevice along a unidirectional trajectory is comparable to the wellknown technique of suspending or manipulating objects by either the application of alternating fields or utilizing a direct field to maintain the static position of the target object [11, 12]. However, in the aforementioned biomedical processes, the ability to either sustain a static position in a moving fluid medium possessing varying components of velocity in multiple coordinate directions, or propel 1 the target object(s) along a predetermined spatial trajectory elevates this problem to one beyond classical levitation. Achieving this type of dynamics requires application of forces with components in multiple coordinate directions. The resulting fields, field gradients, and the means by which to generate them are far more complex than that of the levitation problem, which often requires force generation only along a single direction. A system such as this presents an interesting control problem for which a suitable descriptive model is necessary before a sufficient treatment and analysis may occur. The technique of moving an implant or microdevice by fields from static electromagnets is not a new concept [6, 7, 13]. Indeed a substantial amount of work has been performed using a cubic arrangement of coils to facilitate three dimensional motion. A recent technique as described in [13] utilizes six active coils to impart a predetermined force and torque on a target object. The studied approach powers all coils simultaneously with the three opposite the primary attractive coils being operated with a lower current of reverse polarity. This would result in an increase in the field gradients at the location of the target and in effect raise the magnetic force on the object. While this technique has been explored, calculation of the required currents is quite tedious. The device operates in an open loop control scheme and employs minimization techniques to numerically solve for the six coil currents. The difficulty in the current solution for this system rests with the fact that it is underdetermined (i.e., there are more equations than unknowns). For the scenario of moving a spherical object where the orientation is unconstrained, a solution of six unknowns (the currents) from three equations is needed. Hence, many solutions exist, requiring the use of optimization methods to single out a minimal current set. This work explores the technique of reducing the number of coils to only those which are necessary to generate the force. This results in a determined system, allowing for an exact calculation of the current set. The analysis will be limited to the less complicated problem of moving a sphere two dimensionally in a fluid medium using static electromagnets. The system consists of four electromagnets, of which only two are used at any time to provide planar motion to the particle. The pair of active coils which will be utilized is determined by the force requirement on the sphere. Before exploring the dynamics of this system in detail, a synopsis of the work available in the literature will be presented. 1.2 Overview of Relevant Research The use of magnetic particles in a biomedical application is not new and the available work is quite extensive. Early concepts consisted of using magnetic fluids to perform magnetic tracing and 2 mapping of blood flow systems [14]. This technique is “passive” in nature relying on the high permeability characteristic of the introduced fluid for detection by an electromagnetic pickup array. A second concept explored by [14] involves actively guiding and controlling an unmixed and undiluted material in the bloodstream. The envisioned technique would be helpful in sealing off defects in the circulatory system allowing for the initiation of clotting and scarring. Work on this has been limited with very little success due to the nature of ferrofluidic material not being easily controlled thru magnetic propulsion methods. More recent proposals involve manipulation of particles rather than a fluid based material. The envisioned techniques by which an active guidance system utilizing magnetic propulsion may be employed consist of targeted drug delivery and guidance of micro devices in blood vessels. For the technique of targeted drug delivery, the vast majority of present methods utilize uniform magnetic fields in concert with fluid flow dynamics. In this method, the magnetic particles are moved to the desired location by way of fluid system and applied uniform magnetic fields are then applied to hold the particles in the desired location [1, 2, 5]. While these techniques are promising, a more robust solution would be the employment of an active system to guide the particle to the desired location, independent of the flow dynamics. While this would be the ideal system, present technology is somewhat limited due to the fact that extremely small particles require large field gradients to produce a significant degree of body force on the object. While this does not preclude the use of a magnetic propulsion system, the large power requirements necessary for implementation can be a significant engineering challenge. Guidance of larger scale particles and microdevices is a concept on which magnetic propulsion and guidance techniques are readily applicable. Previous investigations of magnetic guidance of microdevices have been undertaken. Once such study has been given in [9, 10]. In this work, an initial investigation to evaluate the feasibility of utilizing a conventional magnetic resonance imaging machine to impart large enough forces on a ferromagnetic sphere was performed. It was the intent of this study to assess the potential for employment of a microdevice in a human circulatory system, therefore the size of the sphere was chosen accordingly. The investigators utilized a 1/8” 1010 carbon steel sphere immersed in water and subjected to uniaxial flow in a 1/4” diameter tube. The MRI was utilized to keep the global position of the sphere static. It was determined that with a field gradient of 0.018T/m, the sphere could withstatnd a maximum flow of 0.370 ± 0.0064l/min. This value is comparable to that of the human circulatory system a good distance away from the heart, but still significantly larger than the flow rates in regions furthest from the heart. As previously stated, the fluid flow only was present in a single coordinate direction therefore evaluation of device guidance in multiple directions was 3 not achieved in this work. In addition, the static position control was achieved in an open loop fashion with the current in the electromagnet coils adjusted accordingly to arrive at equilibrium of the sphere in the flow. Feedback position control of a microdevice was not demonstrated in this work. A second novel concept of a microdevice guided by applied magnetic fields investigated in recent years is available in [15, 16]. This work studies motion of a “screw” type device propelled through both fluid and organic tissue by application of a rotating magnetic field. The device is essentially comprised of a cylindrical magnet with a spiral blade residing on its surface, analogous to a mechanical screw or auger. The cylindrical magnet is magnetized in a direction perpendicular to its central axis (a radial direction). The magnetization in this coordinate direction allows for a torque about the central axis of the cylinder upon application of a magnetic field perpendicular to the central magnet axis. The torque about the cylinder will be present until the magnetic moment of the cylinder aligns itself with the applied field. Subsequent rotation of the applied field results in a rotation of the cylinder about its central axis. The forward (or backward) motion of the device then occurs in the direction of the central cylindrical axis and is a product of the mechanical interaction of the spiral blade with the medium in which the device is immersed. The mechanism may be viewed much in the same way as that of an auger boring its way through a material body. Directional changes of the device motion are facilitated by a change in direction of the axis of rotation of the rotating magnetic field. Experimental investigations of this concept utilized a 15mm long device with a diameter of 1.2mm and a 0.15mm spiral blade with a pitch of 45 . Results of this study were quite interesting with maximum velocities achieved valued at approximately 20mm/s and 0.022mm/s in fluids with kinematic viscosities of 500mm2/s and 5 × 105mm2/s, respectively. Additionally, maneuvering of the auger device was accomplished in a mazed patterned waterway with 15mm wide passage and filled with a silicone oil having a viscosity of 1000mm2/s. Finally, additional experimental work has been performed evaluating the ability of a device such as this to “bore” its way through gel or a soft solid. Results in this area were quite satisfactory for the initial study into this technique. While this is a novel technique, it differs subtly from the method proposed in this work. The method proposed in [15, 16] utilizes an applied torque and a mechanical “propeller” effect to achieve device motion. This is different than the technique to be investigated by this work in which the direction of motion is determined by the applied field and the device is effectively pulled through the medium in which it resides. The body of work most closely related to the study presented in this investigation may be found in [17, 6, 7, 13, 18]. The cited research investigates a system known as the magnetic stereotaxis 4 system. This system is capable in theory of moving a small intracranial implant along a path of minimum neurological damage in the brain by utilizing applied field gradients. The insert is then positioned inside a deepseated tumor and inductively heated, hence heating and eliminating only the tumor. The envisioned implant would have a dimensional scale < 3mm. The initial investigations of this system incorporated a single electromagnet situated on the end of a robotic manipulator. The imaging method employed by this technique consist of a biplanar fluoroscope which was chosen due to its insensitivity to the applied magnetic fields. Further work into this apparatus emphasized the shortcomings of the moveable electromagnet technique and a system employing an array of stationary electromagnets was the selected path of continued study. More specifically, the stationary electromagnet assembly is a cubic arrangement on which the central axis of each electromagnet is collinear with the global coordinate system. This system is similar to the electromagnetic actuator technique explored in this work. Indeed the mechanism by which the magnetized seed is moved is the method studied. The target object is essentially pulled through the medium in which it resides. For this system however, all six electromagnets are energized. This operation technique results in an underdetermined system in current in which more variables are present than equations available to solve. As a result, optimization techniques must be employed to evaluate and search for a current solution satisfying the desired force on the target object. One advantage to this approach of utilizing all electromagnets is the ability to increase the field gradient at the target location, however it must be recognized that the overall system power is greatly increased due to utilization of all coils. In contrast, this work will examine a technique which will employ only the minimum number of coils required to generate the desired force on the target object. Another aspect of the magnetic stereotaxis system is that it is only capable of operation in an “open loop” mode. The target object is propelled in an impulsive fashion with the location being reached by just “nudging” the seed into position. This is presumably due to the difficulties associated with the optimal calculation of the current vector for the movement instance. By utilizing the minimal coil set method, it can be shown that a direct formulation giving the required coil currents can be made. This will be demonstrated with more detail in a later passage. 1.3 Contributions of this Study There are several contributions produced by this body of work. First, a detailed first principle model of the particle dynamics and electromagnetics will be presented. The electromagnetic model will be one formulated from a simplified coil representation consisting of finite current loops. Use 5 of this method allows for incorporation of the model into a controller in which closed loop control may be achieved. Secondly, a new technique based on the minimal coil set will be presented and explored. Advantages to this technique lie in the fact that a direct formulation may be made to calculate the required coil currents for producing the necessary or desired force on the target object. In addition, the coil switching algorithm will be rooted in a method by which the knowledge of the actual coil currents is unnecessary. Only the geometric field relations of the coils and the target object position is required. Thirdly, an experimental apparatus for investigation of this concept is designed and presented. The test configuration is of “bench top” scale and is optimally designed for minimal use of power. Finally, extensive closed loop experiments in which a small magnetized sphere is guided along 2D trajectories are performed and the results presented. As this is the initial work for achieving closed loop control in 2 dimensions utilizing electromagnetic actuators, a classical control technique is employed along with the minimal coil set current solution method. 1.4 Overview of this Study As has been previously stated, this study will limit itself to the simplified system of a magnetized sphere residing in a fluid medium. Additionally the sphere will be constrained to have only planar motion. Therefore, this study will examine the details associated with a system of four coil which are engineered to provide controlled motion of the sphere in the plane. This body of work will be presented in six separate chapters. Chapter 1 provides the motivation for the work and any existing relevant research. These details are available in the previous passages. The theoretical derivation of the system model will be presented in Chapter 2. In this chapter, the dynamic model will be derived from first principles. It describes the particle dynamics as well as the electromagnetics with the geometric field vector functions presented in a general form. Additionally, a state space representation will be presented. Also, the inverse current solution will be formulated. This will correspond to the system being studied which consists of the particle being constrained to a 2D plane. Subsequently, the technique by which the coil switching will be determined and the resulting calculated current trajectories smoothed will be discussed. A detailed electromagnetic field description will be presented in Chapter 3. This will be a more rigorous treatment of the general geometric vector field functions associated with the electromagnets which were introduced in Chapter 2. The coil field model will be based on a discrete filamentary loop description. Initially, the field resultant from a single filamentary loop will be presented, followed by a multiturn expansion. Subsequently, a formulation describing the net field and field gradient 6 generated by the quad coil assembly will be presented and discussed. Following that, relationships for calculating the electromagnet resistance, self inductance, and mutual inductance will be derived. Finally, experimental verification of the predicted field will be presented and discussed. The design and optimization of the experimental electromagnet assembly will be presented in Chapter 4. Calculation of the electromagnetic characteristics used in the optimization will employ the results of Chapters 2 and 3. A functional for minimization based on minimal power required for moving a magnetized sphere in a circular orbit is derived and presented. Additionally, a minimization technique based on steepest descent is discussed and employed. Finally, the details of the optimally designed electromagnetic coil set are presented. The results of closed loop control of the sphere’s 2D motion are presented in Chapter 5. A closed loop control law utilizing a classical PID controller with and without feedforward in concert with the inverse current solution is presented in Chapter 2. A numerical simulation is formulated and the results utilized to discuss further considerations about the coil switching technique. Additionally, three different path trajectories are examined in simulation. These consist of a spiral beginning at the origin and approaching a circular orbit, a linear path following the outline of a square, and the point regulation problem in which a single coordinate is provided and the path to approach it is not specified. Simulation results for both controllers with all three trajectories are presented and discussed. Finally, the control method was implemented in experiment and the results are presented and discussed. 7 CHAPTER 2 Dynamic System Model 2.1 Dynamic Model Development The equation of motion describing the movement of the target object (sphere) may be devised through a Newtonian formulation given here as m ¨~S (t) = ~FBias + ~FMag + ~FDrag, (2.1) where ~FBias denotes any bias force, ~FMag denotes the applied magnetic force, ~FDrag denotes a fluid drag force from the motion of the object relative to the fluid medium it resides in, m denotes the mass of the object, and ¨~S (t) denotes the vector acceleration of the object in an inertial reference frame. The corresponding trajectory ~S (t) describing the object’s motion may be ascertained by solving the differential equation (2.1). For the system consisting of a spherical target residing in a fluid medium, the first term on the right hand side of (2.1) may be defined to represent the resulting force of gravity on the sphere and the buoyant force of the fluid medium on the sphere ~FBias = ( sph − fl)Vsph~g, (2.2) where sph denotes the density of the sphere, fl denotes the density of the fluid medium in which the sphere resides, Vsph denotes the volume of the sphere, and ~g denotes the acceleration of the sphere due to gravitational attraction. The second term on the right hand side of (2.1) may be expressed by approximating the sphere as a magnetic dipole. This approximation for small objects is widely accepted and has seen some experimental correlation [2, 8]. When a magnetic field is applied, a resulting force and torque act on the dipole. In this case, assuming that the fluid medium has a permeability close to that of free space, the force and torque on the sphere are ~FMag(H) = 4 μ0r3 3 M· r H (2.3a) ~ Mag(H) = 4 μ0r3 3 (M×H), (2.3b) where M is the magnetization of the body, μ0 is the permeability of free space (valued at 4 × 10−7T · m/A), r is the radius of the sphere, and H is the applied magnetic field strength at the 8 location of the sphere [19, 20]. Note that boldface indicates a vector field. In the case of the particle residing in a medium where the relative permeability μfl 1, then μ0 ! μ0μfl. Recognizing that rotational resistance of the spherical object will be small in most fluid media, it is apparent from (2.3b) that the body will have a tendency to align with the applied field H. It may be noted that in some applications, particularly those involving an elongated magnetic object as suggested in [13], the alignment of the object relative to the applied force ~FMag may be considered as an additional constraint on the target. In this study, however, this requirement is disregarded due to the symmetry of the object. At this point in the model formulation, the type of material utilized for the spherical target must be addressed, as this affects the magnetization M achieved when the field H is applied. Primarily two types of materials are proposed for use, paramagnetic or ferromagnetic. The paramagnetic class of material is characterized by a linear relationship between the resulting magnetization and the internal magnetic field of the object. This linear constant is referred to as the magnetic susceptibility and is typically denoted as m. The resulting magnetic force on the paramagnetic sphere is [19] ~FMag,PM(H) = 2 μ0r3( m) ( m + 3) r H·H . The second class of material, ferromagnetic, is characterized by a nonlinear relationship between the applied magnetic field and the magnetization of the object. More specifically, the M −H curve is sigmoidal in shape indicative of these materials reaching a magnetic saturation. Additionally, the sigmoidal M − H curve may also possess a region of hysteresis about the origin depending on the classification of the ferromagnetic material as “hard” or “soft”. An example of a typical ferromagnetic M − H curve is available in Figure 2.1. Hard ferromagnetic materials are typified as those in which the M − H hysteresis loop is large and upon removal of the applied H, the sample is permanently magnetized (e.g. a permanent magnet). This is often referred to as magnetic remanence. In contrast, soft ferromagnetic materials have a small hysteresis loop and very little residual magnetization upon removal of any applied magnetic field (e.g. transformer iron) [21]. This work will address only the ferromagnetic material type. For the ferromagnetic sphere residing in a medium with permeability close to that of free space, the resultant magnetic force for a given amount of magnetization is ~FMag,FM(H) = 4 μ0r3 3 M H H· r H, (2.4) where M is the magnitude of magnetization and H is the magnitude of the applied field H at the sphere’s location [19]. For the permanently magnetized sphere, M may be treated as a constant. In 9 M H indicates initial magnetization Figure 2.1: Typical M − H curve for a ferromagnetic material. the Cartesian coordinate system, taking H = HXˆi + HYˆj + HZˆk whereˆi , ˆj , and ˆk are unit vectors in the respective coordinate directions, the magnitude of the applied field is given by H = H = q H2X+ H2 Y + H2Z . For the final force on the right hand side of (2.1), it has been readily accepted in the literature [1, 2, 3, 4, 5, 8, 10] that the flow about a spherical particle or object in these applications corresponds to a low Reynolds’ number condition (i.e. Re 1) and that the drag force may be formulated by applying the Stokes’ flow unbounded free stream drag force relation ~FDrag = 6 rμ(˙~S fl − ˙~S ), (2.5) where μ represents the viscosity of the fluid medium, ˙~S fl represents the free stream velocity and the quantity (˙~S fl − ˙~S ) denotes the velocity of the moving object relative to the free stream velocity. Indeed there is experimental evidence that, at least for the cases where particles are manipulated in a stagnant or slowly moving medium, this approximation correlates with the observed phenomena [2, 8]. It must be acknowledged, however, that flow scenarios may exist in which the Stokes approximation (2.5) may lose some accuracy. For instance, the Reynolds’ number for flow over a sphere 10 is Re = 2 fl ˙S flr μ , where fl is the density of the fluid. The extreme smallness of the particle may still allow for a calculation of Re 1 even with a significant magnitude of flow velocity around the sphere. If the object is of a larger size (e.g. a microdevice or catheter guide), the result of Re 1 may no longer hold, pushing the drag coefficient into a nonlinear function of Re. Additionally, the relative size of the sphere to that of the channel in which it is moving and the geometry of that channel will have an impact on the calculation of the drag coefficient. As an example, consider the system of a sphere moving along the central axis of a cylindrical duct or a square duct. For the cylindrical duct, the drag force on the sphere is given as FDrag,cyl = 6 μUr 1 + k r R0 , (2.6) where U is the velocity of the sphere, R0 is the radius of the cylindrical duct, and k is a correction coefficient valued at k = 2.10444 [22]. For the square duct, the drag force on the sphere is given by FDrag,sqr = 6 μUr 1 + k r l , (2.7) where l is the half width of the square duct and k = 1.903266 [22]. Likewise, the wall effects on an object in close proximity to the channel surface have been neglected. These will influence the drag force and often impart lift and moment components on the sphere. Depending on the flow boundary conditions, often the most convenient way to evaluate these effects is through numerical simulation [23, 24]. Additionally, it is well known that fluids in living bodies do not behave in a “Newtonian” fashion. Additional forces may be present, i.e. collisions with other small bodies. Depending on the size of the target object, these dynamics may need to be accounted for. Finally, the condition of turbulent flow in the medium is neglected. Turbulent flow of the fluid medium will lead to significant deviation from the Stokes’ approximation. The wall and turbulent flow effects will be ignored for the purposes of this study and left for analysis at a later time. For now, equations (2.2), (2.4), and (2.5) may be used in (2.1) to solve for the motion of the sphere. A component which remains to be addressed is the method to generate the required field H in (2.4). The technique utilized in this study for generation of the magnetic field on the sphere is that of an air core electromagnet. It will be assumed that the required motion of the target and hence the required field changes will be slow enough that the magnetic quasistatic approximation may be used to ignore any transient effects (e.g. eddy currents, etc.). In this system, stationary electromagnets will be employed. Because of this, multiple coils must be used to facilitate parallel and antiparallel 11 forces along multiple coordinate directions. This is due to the fact that when the sphere is unrestrained in the rotational sense and allowed to freely align with the field, a single electromagnet will act in an attractive capacity alone. For three dimensional motion an obvious configuration is to make use of six coils arranged in a cubic fashion similar to that of the magnetic implant guidance system by Stereotaxis, Inc. [6, 7, 13]. The resulting field from a closed path current density conforms to a solution of Maxwell’s field equations which are readily available in [19, 20]. Application of these equations and subsequent reduction of this set to the law of BiotSavart yields a relation expressing the field H of the Pth individual coil as a linear function of the current applied to the coil HP = jP [GXˆi + GYˆj + GZˆk] = jPGP (X, Y,Z), (2.8) where jP is the applied current,ˆi , ˆj , and ˆk are unit vectors in the respective coordinate directions, and Gp is a nonlinear vector field multiplier, which is a function of the location of the sphere relative to the coil and the geometry of the coil itself. The combined field from the employment of Q electromagnets may be ascertained through superposition of the individual coil fields H = [J]T 2 66664 G1 ... GQ 3 77775 = [J]T [G], (2.9) where [J] is a Q × 1 column vector of the applied currents ([J] 2 RQ) and GP is the geometry and position dependent vector field from the Pth coil evaluated in an inertial coordinate system. Using (2.9) the following quantity in (2.4) is derived M H H· r H =M H [J]T [kFM,X]ˆi + [kFM,Y ]ˆj + [kFM,Z]ˆk [J] =M H [J]T [~kFM][J], (2.10) where kFM,X = [GX] @GX @X T + [GY ] @GX @Y T + [GZ] @GX @Z T , kFM,Y = [GX] @GY @X T + [GY ] @GY @Y T + [GZ] @GY @Z T , kFM,Z = [GX] @GZ @X T + [GY ] @GZ @Y T + [GZ] @GZ @Z T , H = r [J]T h [GX][GX]T + [GY ][GY ]T + [GZ][GZ]T i [J], 12 [GS] = 2 66664 GS,1 ... GS,Q 3 77775 , (2.11a) @GS @S T = @GS,1 @S , · · · , @GS,Q @S . (2.11b) Note that in equations (2.11a) and (2.11b), S = X, Y,Z. All field components and partial derivatives of the field components are evaluated in a global inertial coordinate system. Substituting (2.2), (2.4), (2.5), and (2.10) into (2.1) yields the equation of motion m ¨~S = ( s − fl)Vs~g + 4 μ0r3 3 M H [J]T [~kFM][J] + 6 rμ(˙~S fl − ˙~S ). (2.12) Observe that (2.12) relates a given column vector of input current to the resultant motion of the sphere. Because this system utilizes electromagnet coils, the inductive effects of the electromagnet array must be taken into account. Hence, the equation of motion must be augmented by a suitable description of the electrical dynamics for a complete description of the system. Approximating each coil as an RL circuit (a resistor in series with an inductor) results in the following relationship between the input voltage and current flow of the electromagnets vm = Rmjm + XQ n=1 Lm,n djn dt , (2.13) where Q is the number of coils, m = 1, 2, . . . ,Q, n = 1, 2, . . . ,Q, vm denotes the input voltage of coil m, Rm denotes the resistance of coil m, jm denotes the current in coil m, and djn/dt denotes the time rate of change of the current in coil n. Additionally, Lm,n signifies the mutual magnetic inductance of coil n on coil m or the self inductance of coil m when n = m [19]. It should be noted that while the self inductance is always positively valued, the sign of the mutual inductance is dependent on the orientation of coil m to n. Observing that (2.13) is a linear system, it may be rewritten in state space form [J˙] = −[L−1][R][J] + [L−1][V ] (2.14) where [V ] and [J] are Q × 1 column vectors of the voltages and currents in Q coils, [R] is a Q × Q diagonal matrix of the coil resistances, and [L] is a Q×Q symmetric matrix of inductance with the diagonal elements being the self inductance terms and the off diagonal elements being the mutual inductance terms. Note that [L] is symmetric because Ln,m = Lm,n and also note that [L] is invertible [19]. Equation (2.12) is used in concert with (2.14) to provide a full description of the sphere subject to Stokes’ flow and a magnetic field generated by Q discrete electromagnets with input voltages [V ]. 13 Z X Y x1 x3 x2 x4 y1 z1 y2 z2 z y 4 4 y3 z3 plane of motion (YZ) j1 j2 j3 j4 1 4 3 2 c0 sphere Figure 2.2: Schematic of a four coil arrangement to provide forces for planar motion of a sphere. The individual coil coordinates indicate direction of positive current polarity (right hand rule). 2.2 Inverse Current Solution for Constrained 2D Motion with a Reduced Coil Set In an effort to further explore a system using stationary coils to manipulate the position of a magnetized sphere in a fluid medium, this study will examine the two dimensional planar motion case with the sphere constrained to the Y −Z plane (see Fig. 2.2). As previously mentioned, a cubic arrangement of coils can facilitate three dimensional motion of the target sphere [6, 7, 13]. Similarly, a set of four coils arranged on the edges of a plane with the central axes of the coils lying in the plane will facilitate forced motion in that plane. A schematic of this type of system is shown in Fig. 2.2. The inverse current solution entails determining the current trajectories [J] which would cause the sphere to move along a desired motion trajectory ~S . Indeed, the inverse current solution may be considered as a more thorough examination of (2.4), acknowledging in this case that ~FMag is the force required to “offset” the inertial, bias, and drag forces incurred by the sphere moving along ~S . It is obvious from (2.4) that ~FMag,FM is a function linear in [J] but is nonlinear in the components of H, so a direct algebraic solution is not feasible. With this noted, it is evident that in the configuration where the sphere is constrained to a central plane and when the coil configuration to provide both parallel and antiparallel motion is employed, the system is underdetermined. In other words, a 14 θ FZ FY fMag,FM,Y fMag,FM,Z FMag,FM Figure 2.3: Graphical representation of the applied planar magnetic force. solution in four variables (the coil currents) is to be ascertained from just two equations (the vector components of (2.4) in the Y and Z directions). Due to the fact that electromagnets will act in an attractive capacity alone, one may reduce the number of utilized coils to only those which are necessary for the required force. In the scenario of planar motion utilizing the system in Fig. 2.2, this entails reducing the active coils to that of “corner” combinations (i.e. coils 1&4, 1&3, 2&4, and 2&3), commensurate with utilizing only adjacent coils. This will reduce the number of solution variables from four to two and make the system determined. For example, suppose that there is a known required magnetic force ~FMag,FM in the Y − Z plane. This force may be represented in a cylindrical coordinate system as ~FMag,FM = fMag,FM[cos( )ˆj + sin( )ˆk] where fMag,FM is the force magnitude and gives the angular direction of the force measured from the (global) Y axis with the counter clockwise direction taken as positive (see Fig. 2.3). Recognizing that the coils will only attract the sphere, for a required force in the direction of = 45 , the necessary coils to achieve this will be numbers 1 and 4 (see Fig. 2.2) with the actual current amplitudes and polarity dependent on the location of the sphere. In other words, the coil combination for a force pointing into a given quadrant will be the coils which bound that quadrant. One may postulate that the corner combinations of 1&4, 1&3, 2&3, and 2&4 correspond to force directions of 0 90 , 90 180 , 180 270 , and 270 360 , respectively. However, this quadrant 15 assignment to pairs of coils is approximate due to the nonlinear nature of the attractive magnetic fields. The components of the required force ~FMag,FM may be expressed as a function of one another fMag,Z = tan( )fMag,Y , (2.15) where once again is the angular direction of the required force measured from the Y axis with the counter clockwise direction taken as positive. As can be seen from (2.4) and (2.12), the Y and Z components of magnetic force for the ferromagnetic sphere are fMag,Y = 4 μ0r3 3 M H [J]T [kFM,Y ][J], (2.16a) fMag,Z = 4 μ0r3 3 M H [J]T [kFM,Z][J]. (2.16b) Substituting (2.16a) and (2.16b) into (2.15), rearranging, and then simplifying yields [J]T [ ][J] = 0, (2.17) where [ ] = 8>>< >>: [kFM,Y ], for = 2 ± n , n = 0, 1, 2. . . . ; [kFM,Z] − tan( )[kFM,Y ], otherwise. The quantity [J] is a 2×1 vector of the currents in the reduced coil set and the vector fields [kFM,Y ] and [kFM,Z] are 2×2 matrices containing the G and @G/@S terms of the reduced set in accordance with (2.11a) and (2.11b). The matrices [kFM,Y ] and [kFM,Z] are given by the quantities in (2.10), however GS and @GS/@S are now 2 × 1 matrices with the first element corresponding to coil 1 or 2 and the second element corresponding to coil 3 or 4, depending on the coil combination being examined. Once again note that all field components and partial derivatives of the field components are evaluated in the global inertial coordinate system. Equation (2.17) corresponds to a quadratic surface in the J (current) space, centered about the origin. For a current solution other than J = [0 0]T to exist, [ ] must be indefinite or semidefinite [25]. In other words, for a nontrivial solution to exist [ ] cannot be positive or negative definite. Several techniques are available to ascertain the definiteness of [ ]. One method involves examining the signs of the principle minors of [ ]. If the principle minors are not exclusively < 0 or > 0, then [ ] is semidefinite or indefinite [25]. Alternatively, for the 2×2 matrix [ ] one may examine the product of eigenvalues. If the eigenvalue product is 0, then [ ] is semidefinite or indefinite. To facilitate the inverse current solution, it is necessary to require the eigenvalues to be real. This may be achieved by forcing symmetry through the relation [J]T [ ][J] = [J]T [( + T )/2][J] = [J]T [ ][J]. 16 Recognizing that the 2×2 matrix [ ] of the reduced coil set will yield two eigenvalues, the eigenvalues of [ ] will be denoted as 1 and 2. For a sphere located in the Y −Z plane, a force may be effected on the sphere by any coil combination with an eigenvalue product of 1 2 0. As an example, the sign of the eigenvalue products for coil combination 1&4 was examined over a range of −45 135 for three different locations in the Y −Z plane. The initial system analyzed and demonstrated here consists of four rectangular filamentary current loops with a major dimension of 1 unit, a minor dimension of 0.5 units, and an offset “c” from the global origin (see Fig. 2.2) of 1.5 units. The effective force angle range evaluated at the Y − Z locations of (0.5, 0.5), (0.0), and (−0.5,−0.5) are −30.51 120.51 , −2.43 82.43 , and 10.89 79.11 , respectively. These angular ranges indicate the directions of force that coil combination 1&4 is capable of achieving at the given sphere locations. As alluded to earlier, the range of force directions a particular coil combination may impart on the sphere is dependent on the nonlinear field expressions which are functions of the sphere’s position as well as coil geometry. The effective range of force direction decreases as the sphere moves away from the electromagnet pair. To solve for the required currents to achieve the desired (planar) force ~FMag,FM, the orthonormal basis [eQ ] composed of the normalized eigenvectors of [ ] may be used to form the orthogonal transformation [J] = [eQ ][ e J]. Application of this to [J]T [ ][J] = 0 results now in the diagonal quadratic [ e J]T [eD ][ e J] = 0 where [eD ] is a diagonal matrix of the eigenvalues of [ ]. In the e J space, two line loci may now be formed from the diagonal quadratic ˜j 2 i 1 +˜j 2 ii 2 = 0 8>>< >>: ˜j ii = ±˜j i r − 1 2 ; 2 6= 0, (2.18a) ˜j i = ±˜j ii r − 2 1 ; 1 6= 0, (2.18b) where the subscripts i = 1, 2 and ii = 3, 4 denote which coil is used in the calculation. As a reminder, recall that [ ] is also a function of which coils are included in the reduced coil set. Note that (2.18) actually provides two separate formulations for the line loci with a caveat on the existence of a 17 nonzero eigenvalue. Utilizing (2.18), the [ e J] space current vector may now be rewritten as [ e J] = 8>>>>>>>>< >>>>>>>>: ˜j i 2 64 1 ± q − 1 2 3 75 =˜j i[ ± 1 ]; 2 6= 0, (2.19a) ˜j ii 2 64 ± q − 2 1 1 3 75 =˜j ii[ ± 2 ]; 1 6= 0. (2.19b) Applying the orthogonal transformation to (2.16a) yields fMag,Y = 4 μ0r3 3 M eH [ e J]T [eQ ]T [kFM,Y ][eQ ][ e J], (2.20) where eH = q [ e J]T [eQ ]T [G][eQ ][ e J], [G] = GXGT X + GY GTY + GZGT Z. Substitution of (2.19a) into (2.20), simplifying, and utilizing (2.18a) results in the [ e J] solution of ˜j i = ± 3fMag,Y 4 μ0r3M q [ ± 1 ]T [eQ ]T [G][eQ ][ ± 1 ] [ ± 1 ]T [eQ ]T [kFM,Y ][eQ ][ ± 1 ] , (2.21a) ˜j ii = ±˜j i r − 1 2 . (2.21b) Note that (2.21) must still conform to the caveat of 2 6= 0. For the condition that 2 = 0, a similar process may be employed utilizing (2.18b) and (2.19b) to provide the solution ˜j i = ±˜j ii r − 2 1 , (2.22a) ˜j ii = ± 3fMag,Y 4 μ0r3M q [ ± 2 ]T [eQ ]T [G][eQ ][ ± 2 ] [ ± 2 ]T [eQ ]T [kFM,Y ][eQ ][ ± 2 ] , (2.22b) where 1 6= 0. Additionally, solutions may be formulated using (2.16b). This will result in forms similar to (2.21) and (2.22) with fMag,Y ! fMag,Z and [kFM,Y ] ! [kFM,Z]. To select the proper signs for the solution coordinates, the correct groupings on loci intersections must be observed. These are available in Table 2.1, where in the solution value˜j ( , = ±, = i, ii), = ± denotes the sign just to the right of the “=” symbol in (2.21) or (2.22) and = ± denotes the sign used for [ ], as defined in (2.19). It can be seen from (2.21), (2.22), and Table 2.1 that for the planar motion reduced set, there are four different solution coordinates. The coordinate system e J is a transformed system and therefore has no physical meaning. To complete the calculation of the current solution, the e J coordinate solutions must be transformed back to the J space via [J] = [eQ ][ e J]. 18 Table 2.1: Correct Line Loci Groupings for e J Space Solution Coordinates solution set sign e Ja (˜j ++ i ,˜j ++ ii ) e Jb (˜j −+ i ,˜j −+ ii ) e Jc (˜j +− i ,˜j −− ii ) e Jd (˜j −− i ,˜j +− ii ) A plot of the intersecting loci indicating the solution coordinates in the J space for the previously mentioned filamentary loop configuration, a sphere with a unit parameter multiplier (4 μ0r3M/3 = 1) located at the global origin, and a required magnetic force of ~FMag = 1ˆj + 1ˆk is shown in Fig. 2.4. As can be seen, the four solution coordinates may be ascertained by the intersection of the component loci of equation(s) (2.16) or the intersection of the line loci for equation (2.17) with just one of the component loci of (2.16). Indeed it is the latter technique which yields the results of (2.21) and (2.22). The calculated solution coordinates using (2.21) for solutions Ja, Jb, Jc, and Jd are (11.208, 11.208), (−11.208,−11.208), (−25.519, 25.519), and (25.519,−25.519), respectively. As can be seen, they correlate with the loci intersections. The inverse current solution presented in this work was derived for the ferromagnetic material type. A similar derivation may be used for the paramagnetic material and will give the same number of solutions, four. However, the solution is of a simpler nature because the shapes of the component loci for a given force are analytically known. In this case the current solution coordinates correspond to intersecting conic sections. Extending this line of thought would result in the inverse current solution for the paramagnetic material being manipulated by a cubic arrangement of coils being the intersection of quadric surfaces. While the inverse current solution given in relations (2.21) and (2.22) do indeed result in coil currents providing the desired magnetic force on the sphere, as will be shown the sign convention given in Table 2.1 is not sufficient to provide a smooth current trajectory as a function of time. It is advantageous however for the current profile to be smooth to minimize inductive effects incurred by discontinuities in the desired current profile. Further analysis into this formulation is necessary to achieve this and will be provided in the next section. 19 −50 −40 −30 −20 −10 0 10 20 30 40 50 −50 −40 −30 −20 −10 0 10 20 30 40 50 j 1 j 4 line loci f Mag,Z locus f Mag,Y locus Figure 2.4: Solution loci for the 1&4 coil combination in the J space. 2.3 Current Trajectory Smoothing Technique To formulate a more detailed analysis of the inverse solution current coordinates using (2.21) and (2.22), a more precise representation of the desired force on the sphere is needed. As previously stated, the required force on the sphere may be regarded as the force required to overcome the inertial, drag, and buoyant forces acting on the sphere. From (2.12), the required force to be provided by the electromagnets for a given trajectory ~S is ~FMag,FM = m ¨~S + 6 rμ(˙~S − ˙~S fl) + ( fl − s)Vs~g. (2.23) In addition, a trajectory must be examined to fully understand the resultant current profiles calculated using the inverse current formulation presented. For this example, the trajectory to be used was chosen to be a circular path in the Y − Z plane, orbiting the origin with a frequency of 0.1Hz. This may be defined as ~S (t) = cos 2 10 t ˆj + sin 2 10 t ˆk , (2.24) where is the radius of the path and t represents time. As can be seen, the first and second time derivatives of (2.24) may be easily evaluated and substituted into (2.23) to calculate the necessary force ~FMag,FM. It will be assumed that the vector acceleration due to gravity, denoted as ~g, will be 20 normal to the plane of motion, hence offset by a constraint force and may be subsequently dropped from (2.23). Furthermore, it will also be assumed that this testing apparatus will move the sphere in a static fluid, allowing for the dropping of the term ˙~S fl in (2.23) as well. The physical parameters in (2.21) (or (2.22)) and (2.23) which need assignment consist of the mass of the sphere, the radius of the sphere, the viscosity of the fluid medium, the radius of the circular path, and the magnetization of the sphere. In an effort to formulate a more “realistic” example, parameters closer to that of the physical system will be utilized. For this system, m = 0.1257g, r = 1/8”, = 0.5”, and μ = 0.950. It should be noted that the fluid selected for initial testing is glycerine and the viscosity value is consistent with that found in many fluid property tables. The sphere used in this study is a permanent neodymium rareearth magnet of grade 40. Accounting for demagnetization effects of the spherical shape results in a permanent magnetization of M = 3Br 2μ0 , where Br is the remanent flux density in the magnetized sphere [19, 21]. For this study, the remanent flux density will be estimated to have the value of Br = 1T. A selection of this value is consistent with the material type and grade of magnet. Finally, the coil geometry used for this example is that of rectangular coils wound with 12 gauge magnet wire. The coils consist of 21 layers of wire with 18 turns per layer. The inner dimensions of the coils are 0.439” × 0.757”, while the distance from each coil face to the global origin is 2.5”. A more in depth analysis of this coil configuration will be provided in Chapter 4. The trajectory chosen for this study is that of a circle orbiting the origin with a constant angular velocity. Because of this, symmetry may be used and only one coil combination need be examined. In other words, the force requirements for each coil combination is the same, hence the current solutions are of the same form. Utilizing this, the combination of coils 1 and 3 will be examined initially. A plot of the four solutions given by (2.21) (see Table 2.1) for coil combination 1 and 3 is shown in Figure 2.5. As can be seen, the solution provided by (2.21) for coils 1 and 3 has two distinct regions in which [ 13] is semidefinite or indefinite. However, if this solution is used to calculate the resultant force generated and then compared to the required force it can be seen that one of the current regions is erroneous. This is demonstrated in Figure 2.6. The incorrect solution regions of both force components seen in Figure 2.6 are artifacts of the mathematical formulation of the inverse solution technique. As is evidenced by Figure 2.6, the requirement of [ ] to be semidefinite or indefinite is not “strong enough” to result in a true solution. The resulting inverse solution must also satisfy the desired direction of the applied magnetic force. In other words, calculated must 21 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J a (A) j 1 j 3 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J b (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J c (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J d (A) time (sec) Figure 2.5: Initial current solution for coil combination 1 and 3 during circular trajectory. 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,a (N) f req f mag 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,a (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,d (N) time (sec) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,d (N) time (sec) Figure 2.6: Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces. 22 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J a (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J b (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J c (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J d (A) time (sec) j 1 j 3 Figure 2.7: Correct current solution for coil combination 1 and 3 during circular trajectory. equal desired. Disregarding the solution regions which do not satisfy this second criteria, Figures 2.7 and 2.8 may be constructed which give the correct solution and calculated resultant forces for coil combination 1 and 3. As can be seen, the region of incorrect force results has vanished from the calculated solution. Now that a known correct solution with no erroneous sections has been formulated, a second facet of this result may be observed in Figure 2.7, that being the discontinuities present in the current trajectories. Indeed, this aspect when examined at face value may be construed as quite disturbing. While this current trajectory would in fact result in the necessary magnetic force being imparted on the sphere, the ability to achieve this trajectory in an inductive system may or may not be entirely possible. Because of this, it would be desirable for the current trajectories to be “smooth” or continuous. Upon examination of the discontinuity shown in Figure 2.7, it has been observed that the discontinuous point corresponds to the location at which the transformation [eQ ] changes from a rotation to a reflection, or vice versa. It is well known that the particular type of transformation (either rotation or reflection) may be ascertained by examining the sign of the determinant of the transformation [26, 25, 27]. Defining D = det[eQ ] and plotting sgn(D1,3) results in Figure 2.9. A comparison of Figures 2.7 and 2.9 shows that the location of the current solution 23 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,a (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,a (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,b (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,c (N) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Y,d (N) time (sec) 0 2 4 6 8 10 !5 0 5 x 10 !4 F Z,d (N) time (sec) f req f mag Figure 2.8: Calculated forces using initial current solution compared to required force. Note the regions of erroneous calculated magnetic forces are now absent. 0 1 2 3 4 5 6 7 8 9 10 !1.5 !1 !0.5 0 0.5 1 1.5 ,./012/34 25617 183 4 Figure 2.9: A plot of sgn(D1,3), where D = det[eQ ]. A value of +1 indicates that transformation is a rotation, while a value of −1 is indicative of a reflection. 24 Table 2.2: Modified Sign Groupings for Smooth e J1,3 = ˜j 1 1 1 ,˜j 3 3 3 Solution Coordinates solution set 1 1 3 3 e Ja,1,3 +sgn(D1,3) +sgn(D1,3) + +sgn(D1,3) e Jb,1,3 −sgn(D1,3) +sgn(D1,3) − +sgn(D1,3) e Jc,1,3 +sgn(D1,3) −sgn(D1,3) − −sgn(D1,3) e Jd,1,3 −sgn(D1,3) −sgn(D1,3) + −sgn(D1,3) discontinuity corresponds to the point at which the transformation eQ 1,3 changes from a rotation to a reflection. This dependence may be accounted for by modifying the definitions of the correct solution signs given in Table 2.1. A “smoothed” version of the solution given in Figure 2.7 may be achieved by the solution signs given in Table 2.2. Note that the pattern of ±’s preceding the signum functions matches that seen in the sign pattern given in Table 2.1. Utilizing the signing pattern of Table 2.2, a smooth current solution for coil combination 1 and 3 may be formed. This is shown in Figure 2.10. As can be seen, there are two types of solution, one in which the coil currents of like sign and one in which the coil currents are of opposite sign. Additionally, for each of the two types of solution there are two forms in which the current polarity of each coil is reversed. As stated earlier, for the generic circular trajectory chosen here, the current trajectory solutions for the other coil combinations will be of similar form to that shown in Figure 2.10. However, the sign convention for providing a smooth current trajectory in each coil combination will differ from that given in Table 2.2. For the remaining coil combinations 2&3, 2&4, and 1&4, the resulting inverse solutions using equation (2.21) and the sign convention given in Table 2.1 results in the current trajectories shown in Figures 2.11, 2.12, and 2.13. In addition, the required sign convention needed for each of the combinations to achieve a “smooth” trajectory are given in Tables 2.3, 2.4, and 2.5, with the resulting smooth trajectories also shown in Figures 2.11, 2.12, and 2.13. As can be seen, the discontinuities have been removed by using the sign conventions provided in Tables 2.2, 2.3, 2.4, and 2.5. It should be noted that unlike the solution from equation (2.21), the sign convention given in the tables is not dependent on the solution form used. It has been previously stated that (2.21) is equally valid by substituting fMag,Y ! fMag,Z and [kFM,Y ] ! [kFM,Z]. In the case of sign convention, fMag,Y is utilized in the formulation regardless of the use of the Y direction solution versus the Z. Now that smooth current trajectories have been calculated, it is matter of 25 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J a (A) 0 1 2 3 4 5 6 7 8 9 10 !5 0 5 J b (A) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 J c (A) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 J d (A) time (sec) j 1 j 3 Figure 2.10: Smoothed current solution for coil combination 1 and 3 during circular trajectory. selecting from the four different solutions provided by (2.21). Any of the four solutions presented will result in the desired force, however one may choose to prefer one solution over another based on the polarity or the total power utilized. For the purposes of experiment design and optimization of this design, symmetry may be leveraged and only analysis of the hardware configuration required to generate one coil combination’s solution is necessary. As will be seen in Chapter 4, only the solution for coil combination 1&3 shown in Figure 2.10 will be used in the optimization routine. 2.4 State Space Model and Closed Loop Control Considerations Implementation of closed loop control on this type of system is challenging and has not been explored in literature yet. Even with the reduced set and the exact calculation of the required currents, it has been shown in this work that for the two dimensional case there will be four possible solutions. For the three dimensional case utilizing six stationary coils, the reduced set would number three and the number of solution coordinates would be eight. Because of this, closed loop control employing the reduced set will be challenging as well. Even with that said, the benefits of being able to exactly calculate the possible solutions are obvious when compared with using optimization methods to search out solutions for an underdetermined system. Feedback control will be further complicated 26 0 2 4 6 8 10 !5 0 5 J a (A) 0 2 4 6 8 10 !5 0 5 J b (A) 0 2 4 6 8 10 !5 0 5 J c (A) 0 2 4 6 8 10 !5 0 5 J d (A) time (sec) 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !1 0 1 time (sec) j 2 j 3 j 2 j 3 Figure 2.11: Initial and smoothed current solution for coil combination 2 and 3 during circular trajectory. 0 2 4 6 8 10 !5 0 5 J a (A) 0 2 4 6 8 10 !5 0 5 J b (A) 0 2 4 6 8 10 !5 0 5 J c (A) 0 2 4 6 8 10 !5 0 5 J d (A) time (sec) 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !5 0 5 time (sec) j 2 j 4 j 2 j 4 Figure 2.12: Initial and smoothed current solution for coil combination 2 and 4 during circular trajectory. 27 0 2 4 6 8 10 !5 0 5 J a (A) j 1 j 4 0 2 4 6 8 10 !5 0 5 J b (A) 0 2 4 6 8 10 !5 0 5 J c (A) 0 2 4 6 8 10 !5 0 5 J d (A) time (sec) 0 2 4 6 8 10 !5 0 5 j 1 j 4 0 2 4 6 8 10 !5 0 5 0 2 4 6 8 10 !1 0 1 0 2 4 6 8 10 !1 0 1 time (sec) Figure 2.13: Initial and smoothed current solution for coil combination 1 and 4 during circular trajectory. Table 2.3: Modified Sign Groupings for Smooth e J2,3 = ˜j 2 2 2 ,˜j 3 3 3 Solution Coordinates solution set 2 2 3 3 e Ja,2,3 +sgn(fMag,Y ) +sgn(D2,3fMag,Y ) +sgn(D2,3) +sgn(D2,3fMag,Y ) e Jb,2,3 −sgn(fMag,Y ) +sgn(D2,3fMag,Y ) −sgn(D2,3) +sgn(D2,3fMag,Y ) e Jc,2,3 + −sgn(D2,3fMag,Y ) −sgn(D2,3fMag,Y ) −sgn(D2,3fMag,Y ) e Jd,2,3 − −sgn(D2,3fMag,Y ) +sgn(D2,3fMag,Y ) −sgn(D2,3fMag,Y ) 28 Table 2.4: Modified Sign Groupings for Smooth e J2,4 = ˜j 2 2 2 ,˜j 4 4 4 Solution Coordinates solution set 2 2 4 4 e Ja,2,4 + +sgn(D2,4) +sgn(D2,4) +sgn(D2,4) e Jb,2,4 − +sgn(D2,4) −sgn(D2,4) +sgn(D2,4) e Jc,2,4 + −sgn(D2,4) −sgn(D2,4) −sgn(D2,4) e Jd,2,4 − −sgn(D2,4) +sgn(D2,4) −sgn(D2,4) Table 2.5: Modified Sign Groupings for Smooth e J1,4 = ˜j 1 1 1 ,˜j 4 4 4 Solution Coordinates solution set 1 1 4 4 e Ja,1,4 +sgn(D1,4) +sgn(D1,4fMag,Y ) +sgn(fMag,Y ) +sgn(D1,4fMag,Y ) e Jb,1,4 −sgn(D1,4) +sgn(D1,4fMag,Y ) −sgn(fMag,Y ) +sgn(D1,4fMag,Y ) e Jc,1,4 +sgn(D1,4fMag,Y ) −sgn(D1,4fMag,Y ) − −sgn(D1,4fMag,Y ) e Jd,1,4 −sgn(D1,4fMag,Y ) −sgn(D1,4fMag,Y ) + −sgn(D1,4fMag,Y ) 29 for this system due to the inductive effects of the coils. For a given (slow) trajectory, the current solutions to follow that path may be calculated from (2.21) or (2.22). The voltage solutions would then need to be evaluated using (2.13). The method of using a reduced coil set would inherently involve “switching” coils on and off. Inductive effects associated with this type of technique must be accounted and compensated for. One advantage of this system is that for the electrical dynamics, full state feedback is available. Therefore it is conceivable that the switching effects may be taken into account. In the case of large coils with significant inductance parameters, full state feedback of current may prove to be invaluable when dealing with the switching effects. To synthesize a controller, it is beneficial to first express the entire system model in a state space form. For this system, the state variables will be designated as the sphere’s Y and Z positions, the velocity components ˙Y and ˙Z , and currents in the four coils j1, j2, j3, and j4. The inputs for this system are the coil voltages v1, v2, v3, and v4. This results in the state vector x = Y ˙Y Z ˙Z j1 j2 j3 j4 T , (2.25) and the input vector u = v1 v2 v3 v4 T . (2.26) It will be assumed that the measured outputs are the sphere positions Y and Z, and the four coil currents. As can be seen from (2.12) and (2.13), this system may be represented as the two linear systems of sphere dynamics and electrical dynamics coupled by the nonlinear magnetic force terms. In addition any unknown drag force which deviates from the linear Stokes’ static drag may be treated as a disturbance. This leads to the state space representation of x˙ = [A]x + f(x) + [B]u + [D], (2.27a) y = [C]x, (2.27b) where [A] = 2 64 As [;]4×4 [;]4×4 −L−1R 3 75 , [B] = 2 64 [;]4×4 L−1 3 75 , As = 2 66666664 0 1 0 0 0 −6 rμ m 0 0 0 0 0 1 0 0 0 −6 rμ m 3 77777775 , [D] = 2 66666666664 0 dY 0 dZ [;]4×1 3 77777777775 , 30 [C] = 2 64 Cs [;]4×4 [;]4×4 [I]4×4 3 75 , Cs = 2 66666664 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 3 77777775 , f(x) = 2 66666666664 0 4 μ0r3 3m p M xT x xT [KFM,Y ]x 0 4 μ0r3 3m p M xT x xT [KFM,Z]x [;]4×1 3 77777777775 , [KFM,S] = 2 64 [;]4×4 [;]4×4 [;]4×4 [kFM,S] 3 75 , = 2 64 [;]4×4 [;]4×4 [;]4×4 [G] 3 75 . Note that in (2.27), [;]m×n represents an m×n matrix of zeros, [I]m×m represents an m×m identity matrix, and dS, S = X, Y,Z, represents the unknown component flow disturbances. Equation (2.27) is an eighth order nonlinear state space model which may be used for controller synthesis and design. As can be seen, the system is linear in the input u but contains rather formidable nonlinear terms describing the field forces. With the system now represented in a state space form, the formulation of the field components as they pertain to the coil geometries and the calculation of the mutual inductances must be examined in greater detail. This will be explored in the following chapter. 31 CHAPTER 3 Electromagnetic Field Description 3.1 The Single Filamentary Loop While a detailed examination of the planar equation of motion and the relation between all the fields of the individual coils is developed in Chapter 2, the functions describing the resulting field components and gradient field components generated by the rectangular electromagnets is not addressed. Greater insight into this physical system allows for the creation of a model which may then be used in a predictive capacity and allow for design and optimization of an experimental setup for evaluation of this concept. This chapter will expand on the work of Chapter 2 and provide the derivation of a model describing the field generated by elliptical coils wound on a rectangular core, represented as rectangular electromagnets. This geometric simplification allows for the analytic formulation of the descriptive field equations. In addition, the technique used will entail the approximation of the coils as assemblies of discrete filamentary loops [28]. This technique has been selected due to the less difficult analytical treatment when compared to the volumetric coil model [29, 30]. The law of BiotSavart will be used to calculate the field contribution of each loop. To develop a relation describing the resulting magnetic field from the quad coil configuration, an analysis of a single filamentary loop will initially be undertaken. The resulting magnetic field for a single (closed) filamentary loop may be ascertained by using the law of BiotSavart [31], h = j 4 I d~l ×~r r3 , (3.1) where h denotes the magnetic field at a point in 3D space, j denotes the applied current in the filamentary loop, d~l denotes a differential element of the loop, ~r denotes the vector from the differential element to the point at which the vector field is being evaluated, and r is the magnitude of ~r. In this study, we will first examine a single rectangular filamentary loop as represented in Figure 3.1. Application of (3.1) in a piecewise fashion to the rectangular conductor results in h = j 4 " Z b −b d~l 1 ×~r1 r3 1 + Z a −a d~l 2 ×~r2 r3 2 + Z b −b d~l 3 ×~r3 r3 3 + Z a −a d~l 4 ×~r4 r3 4 # , (3.2) where d~l m is a differential element of the m side of the rectangular filament, ~rm is the vector from 32 z y x b a 1 2 3 4 j Figure 3.1: Schematic of a single filamentary conducting loop. The local coordinate system is fixed to the center of the loop with the central axis aligned with the z coordinate direction. Positive current flow obeys the right hand rule. The numbers indicate sections of conductor while the letters denote the dimensions of the rectangle. the m side differential element to the point in 3D space at which the field is being evaluated, and m = 1, 2, 3, 4 indicates a given side of the rectangular filament as shown in Figure 3.1. Note that d~l m is in the direction of positive current flow, which in this case obeys the right hand rule. Evaluation of (3.2) in the rectangular coordinate system of Figure 3.1 results in h(x, y, z, j) = j[gx(x, y, z, a, b, c)ˆi + gy(x, y, z, a, b, c)ˆj + gz(x, y, z, a, b, c)ˆk] = jg(x, y, z, a, b, c), (3.3) where gx = z − c 4 1 '1 y + b '6 − y − b '5 − 1 '2 y + b '8 − y − b '7 , gy = z − c 4 1 '3 x + a '7 − x − a '5 − 1 '4 x + a '8 − x − a '6 , gz = 1 4 x + a '1 y + b '8 − y − b '7 − x − a '2 y + b '6 − y − b '5 + y + b '3 x + a '8 − x − a '6 − y − b '4 x + a '7 − x − a '5 , '1 = (x − a)2 + (z − c)2, '2 = (x + a)2 + (z − c)2, '3 = (y − b)2 + (z − c)2, 33 '4 = (y + b)2 + (z − c)2, '5 = [(x − a)2 + (y − b)2 + (z − c)2]1/2, '6 = [(x − a)2 + (y + b)2 + (z − c)2]1/2, '7 = [(x + a)2 + (y − b)2 + (z − c)2]1/2, '8 = [(x + a)2 + (y + b)2 + (z − c)2]1/2. In addition, (x, y, z) is the point in 3D space at which the field is being evaluated, ˆi , ˆj and ˆk are unit vectors in the x, y, and z directions, a and b indicate the dimensions of the rectangular loop (see Figure 3.1), and c denotes the offset of the loop from the x−y plane of Figure 3.1. It should be noted that there is a fundamental difference between parameters a, b and c. The parameters a and b denote dimensions, hence a, b > 0. However, c denotes an offset of the loop from a plane. In this formulation c > 0 indicates a positive offset, while c < 0 indicates a negative offset. Additionally it may be observed that if the loop resides in the x − y plane as seen in Figure 3.1, c = 0. As can be seen from (3.3), the resultant field from the single rectangular filament is a function of the applied current multiplied by the vector field g. The “g field” is also seen to be a nonlinear function of the 3dimensional location of the point at which the vector field h is being evaluated and the geometry of the filamentary loop. In addition to the vector field h (or g), the calculation of the resultant force on a particle requires additional quantities. It is known in [19] that the vector force on a magnetized sphere approximated as a dipole residing in free space is given as ~FMag(h) = 4 μ0r3 3 M· r h, (3.4) where μ0 is the permeability of free space, r is the radius of the sphere (not to be confused with the magnitude of the vector ~r given in (3.1)) , M is the magnetization vector field of the sphere, and h is the applied magnetic field from the filamentary loop at the location of the particle. It must also be noted that in most fluid media, the spherical object will be unconstrained in the rotational sense resulting in an alignment of the magnetization field with the applied field [32]. Use of this fact along with (3.3) and expansion of the vector quantity (M· r)h in (3.4) results in [33] M h h · r h = jM g " gx @gx @x + gy @gx @y + gz @gx @z ˆi + gx @gy @x + gy @gy @y + gz @gy @z ˆj + gx @gz @x + gy @gz @y + gz @gz @z ˆk # , (3.5) 34 where M is the magnitude of magnetization and g = q g2x + g2 y + g2 z . As can be seen in (3.5), to effectively describe the resultant vector force on the spherical target, not only are the g field components needed, but the partial derivatives are required as well. It should be noted that for a linear paramagnetic sphere, the force relation may be rewritten as ~FMag,PM(h) = 2 μ0r3( m) ( m + 3) r h · h , where m is the magnetic susceptibility of the paramagnetic material. The quantity r(h · h) will then take on the same form as (3.5) with the substitutions of jM/g ! 2j2 and @gm/@n ! @gn/@m, (m, n = x, y, z). The partial derivative terms in (3.5) may be directly evaluated from the g field components given in (3.3). While this has been performed, it proves to be quite tedious and will be omitted here for brevity, however these functions for the rectangular filamentary loop are available in Appendix A. 3.2 Multiturn Expansion To evaluate the resultant vector force on a spherical object from a coil with multiple turns, one may utilize the filamentary model of Section 3.1 and the principle of superposition. If each turn of a multiturn coil is approximated as a single filamentary loop, then the total magnetic field H at a given point in 3D space from a multiturn coil may be given as H = h1 + h2 + · · · + hq = Xq p=1 hp, (3.6) where q is the total number of turns for the coil. Equation (3.6) is only valid with adherence to a few important caveats. First, the individual fields h1, · · · , hq are given in the same coordinate system. For the principle of superposition to hold, the resultant field from each turn must be presented in terms of the same basis. Second, the point in 3D space where the field is being evaluated may not be coincident with a point on a filamentary loop. As can be seen from (3.2), application of the law of BiotSavart for the coincident loop results in an undefined integrand. This is not expected to be an issue due to the fact that no trajectory of the object can be coincident with the coil volume, hence points coincident with the filamentary loops are not needed in this analysis. Application of (3.3) to (3.6) results in H = jg1 + jg2 + · · · + jgq = j Xq p=1 gp = jG, (3.7) where gp is the geometry dependent vector field of the pth loop. Note that in (3.7) the total vector field H is a function of the current j in each turn multiplied by the vector field G which is a function 35 z y x t l d Figure 3.2: Schematic of a rectangular coil with a “perfect wind” packing pattern. The diameter of the wire is denoted as “d”, while the number of turns/layer and number of layers are denoted t and l, respectively. of the number of current loops, the geometry of the current loops, and the location of the loops in 3D space. To formulate the G field expression describing the multiturn coil, use may be made of the expressions in (3.3). As can be seen, (3.3) gives the single loop g field components in the coordinate system of Figure 3.1. This same description may be used for each of the the filamentary loops comprising the approximation of the multiturn coil, provided that the proper values of a, b, and c are used for the loop. Thus to describe the magnetic vector field from a multiturn coil by using the approximate technique of representing the turns as a finite number of current loops, one must only supply the proper values of a, b, and c in equation (3.3) for each pth loop and then apply equation (3.7). If it is assumed that each of the wound coil’s turns may be represented by one concentric loop and that the coil is wound with “perfect” packing, then a schematic of this system may be formed in Figure 3.2. The representation of a wound coil as a grouping of discrete loops is indeed an approximation of the true physical system. In reality, the coil consists of a single conductor wrapped multiple times about a given geometry and the previous “wrappings”. In fact, for a wound coil to be a continuous conductor there must be a region of wire “crossover” which is obviously being neglected in this technique. Additionally, the wound path is in actuality a spiral. This results in an angular 36 component of the current path being present that is not represented by the concentric current loops. This angle can be envisioned as one which would be between the current path and the x − y plane of Figure 3.2. Finally, the ability to fabricate a “rectangular” coil may be quite elusive. While the innermost layers of wire may be very close to a rectangle due to the shape of the bobbin on which the wire is wound, inevitably some rounding of the coil corners will occur. This fact is also being neglected for the analytical model formulation of this study. As is hinted by Figure 3.2, four parameters are necessary to establish the “perfectly” wound, concentric loop model. These consist of the number of wound layers l, the number of turns per layer t, the diameter d of the wire used, and the z direction in which the coil is built. For the representation in Figure 3.2, the coil is built in the positive z direction of the coil’s coordinate system. The parameter t may be thought of as the amount of coil buildup in the z direction while l is the amount of coil buildup in the transverse x − y directions. Note that with the cross sectional wind geometry chosen in Figure 3.2, all layers have the same number of turns per layer. To establish an expression for gp, p = 1, . . . , q, in (3.7), the following relations may be used: G = Xq p=1 gp = Xl =1 Xt =1 g(x, y, z, a , b , c± ), (3.8) where a = a0 + d p 3 2 ( − 1) b = b0 + d p 3 2 ( − 1) c± = 8>>< >>: ±[d( − 1)], for = 1, 3, 5, . . . (odd); ±[ d 2 + d( − 1)], for = 2, 4, 6, . . . (even). In equation (3.8), a0 and b0 denote the initial positions of the innermost layer of the rectangular coil, and the sign of c± indicates the build direction along the ±z axis as indicated in Figure 3.2. Note that q = tl and also that this formulation begins the coil in the local x − y plane, building out from there. An analogous expression may be formed for the partial derivative terms needed in equation (3.5), available here as @Gm @n = Xq p=1 @gm,p @n = Xl =1 Xt =1 @gm(x, y, z, a , b , c± ) @n , (3.9) where m, n = x, y, z. Expressions (3.8) and (3.9) may now be used to calculate the G field and partial derivatives of the G field for a given coil of t turns per layer and l layers (q = tl total turns) at a given point in 3D space in the local coil coordinate system of Figure 3.2. The configuration 37 proposed to facilitate the particle motion consists of four coils, hence an expansion of the single coil model to the quadcoil model is necessary. 3.3 Quad Coil Model The expressions derived in Section 3.2 may also be used to describe the multicoil configuration. A schematic of this system is available in Figure 2.2. Once again using the principle of superposition, the total field resulting from an array of multiturn coils may be given as H = H1 +H2 + · · · +HQ = XQ P=1 HP , where Q represents the number of coils in the array. Note that in system of Figure 2.2, Q = 4. Representing the total field in terms of the G field and current of each coil in the array results in H = XQ P=1 jPGP = j1G1 + · · · + jQGQ = [J]T 2 66664 G1 ... GQ 3 77775 = [J]T [G], where the individual fields GP must be given in a common basis. As has been shown in [32] and Chapter 2, the vector field (M·r)H proportional to the force resulting on a magnetized sphere for an array of coils is derived as M H H· r H = M[J]T H " [GX] @GX @X T + [GY ] @GX @Y T + [GZ] @GX @Z T ˆi + [GX] @GY @X T + [GY ] @GY @Y T + [GZ] @GY @Z T ˆj + [GX] @GZ @X T + [GY ] @GZ @Y T + [GZ] @GZ @Z T ˆk # [J], (3.10) where: H = r [J]T h [GX][GX]T + [GY ][GY ]T + [GZ][GZ]T i [J], [GS] = 2 66664 GS,1 ... GS,Q 3 77775 , @GS @S T = @GS,1 @S · · · @GS,Q @S . As previously stated, the field components and partial derivatives of field components in (3.10) are given in the global coordinate system S = X, Y,Z. To express the fields in the global coordinate 38 system, homogeneous coordinate transformations may be employed on equations (3.8) and (3.9) to transform the calculated field results from the local coil coordinates to the global system [26]. To utilize the representation (3.10) and transform to the proper location in the global coordinate system, the individual local coordinate systems for each coil must be established. The local coil coordinate locations are available in Figure 2.2. It is important to note that the local coordinates given in Figure 2.2 are positioned in the plane of the turns closest to the global origin. Additionally, the correct build directions in the local systems must be chosen, and hence the proper sign of c± in (3.8) and (3.9) selected. For this system it is determined by the orientation of the local coordinate systems; coils one and two utilize c− and coils three and four utilize c+ . As can be seen, local coordinate system 2 is translated and not rotated relative to the global system. Therefore the homogenous coordinate transformation may be given as [A2] = 2 64 [R2] [d2] ; 1 3 75 , (3.11) where [R2] = 2 66664 1 0 0 0 1 0 0 0 1 3 77775 , [d2]T = 0 0 −c0 , the quantity [R2] denotes the rotation component of the transform, [d2] denotes the translation component of the transform, and c0 is the distance from the global origin to the x − y plane of the closest turns of the coil to the global origin. Note that because coil 2 is only translated in relation to the global coordinates, the rotation matrix is effectively an identity matrix. For the other three coils this will not be the case. The coordinate transforms for coils 1, 3, and 4 may be given as [A1] = 2 64 [R1] [d1] ; 1 3 75 , (3.12) [A3] = 2 64 [R3] [d3] ; 1 3 75 , (3.13) and [A4] = 2 64 [R4] [d4] ; 1 3 75 , (3.14) 39 where: [R1] = 2 66664 1 0 0 0 −1 0 0 0 −1 3 77775 , [d1]T = 0 0 c0 = −[d2]T , [R3] = 2 66664 1 0 0 0 0 −1 0 1 0 3 77775 , [d3]T = 0 −c0 0 = −[d4]T , [R4] = 2 66664 1 0 0 0 0 1 0 −1 0 3 77775 . Equation (3.8) may now be used with the relations in (3.11), (3.12), (3.13), and (3.14) to formulate the G field components (given in the global coordinate system) of equation (3.10). This may be accomplished by first representing global point at which the sphere is located in the local basis of each coil, calculating each coil’s local G field and @G/@s contributions (where s = x, y, z), and finally transforming those results back into the global coordinate system for use in (3.10). To transform the point S = (X, Y,Z) to a point sP = (xP , yP , zP ) in the Pth coil’s local coordinate system, it is beneficial to first present the point in its augmented homogeneous representation [26] [ eS]T = X Y Z 1 , [esP ]T = xP yP zP 1 . Transformation of the global point may then be achieved by [esP ] = [AP ]−1[ eS]. (3.15) Once the global coordinate is transformed into the individual coordinate systems, equations (3.3), (3.8), and (3.9) may be used to calculate the local field and spatial partial derivatives of the field. Following this, the local coil field contributions may be transformed back to the global coordinate system by utilizing [GP ] = [RP ][GP ], (3.16) where [GP ]T = GX,P GY,P GZ,P , [GP ]T = Gx,P Gy,P Gz,P . 40 The transform of the spatial partial derivatives of GP is not as elegant as (3.16). This vector field may be transformed back into the global coordinates using [rSGX,1] = [R1][rs,1Gx,1] [rSGY,1] = −[R1][rs,1Gy,1] (3.17) [rSGZ,1] = −[R1][rs,1Gz,1], [rSGX,2] = [R2][rs,1Gx,2] [rSGY,2] = [R2][rs,1Gy,2] (3.18) [rSGZ,2] = [R2][rs,1Gz,2], [rSGX,3] = [R3][rs,1Gx,3] [rSGY,3] = −[R3][rs,1Gy,3] (3.19) [rSGZ,3] = [R3][rs,1Gz,3], and [rSGX,4] = [R4][rs,1Gx,4] [rSGY,4] = [R4][rs,1Gy,4] (3.20) [rSGZ,4] = −[R4][rs,1Gz,4], where rS and rs,P denote the gradient operation in the global coordinates S and local coordinates s, P (P = 1, 2, 3, 4), respectively. Once the G field and partial spatial derivatives of G have been calculated, (3.10) may be utilized to calculate the vector field proportional to the force on the magnetized sphere, provided that the applied current vector [J] is given as well. 3.4 Coil Resistance and Inductance As has been shown in [32] and Chapter 2, the resistance, self inductance, and mutual inductance are needed for an accurate representation of the electrical dynamics of the system shown in Figure (2.2). For a volumetric coil residing in free space with the resulting vector magnetic potential denoted as ~A, the inductance L may be calculated analytically by L = 1 j2 Z V ~A ·~jdV, (3.21) where ~A = 1 4 Z V ~j R dV, 41 ~j is the current density in the coil, V is the volume of the coil, and R is the distance from the volume element dV to the point at which the vector potential is being evaluated [19]. As can be seen, (3.21) is the result of two triple integrals and may prove to be quite challenging. Keeping with the spirit of using a finite loop representation of the field result, the inductance may be formulated in a like manner. Inductance of a coil is defined as L = q j , where is the total magnetic flux linked by the each turn, q is the total number of turns, and j is the current in the coil [31]. Expanding the definition to include the effect between two coils, an examination on a turn by turn basis yields a new definition LM,N = M,N jN = Pq p=1 M,N,p jN , (3.22) where M,N is the total magnetic flux generated by the Nth coil linked by the Mth coil, jN is the current flowing in the Nth coil, and M,N,p is the total magnetic flux generated by the Nth coil linked by the pth filamentary loop of the Mth coil, summed up for q total loops. In this system, the self inductance of a coil may be defined as LM,M and the mutual inductance between two coils as LM,N where M,N = 1, 2, 3, 4 and M 6= N. It is well known that LM,N = LN,M [31, 19, 20]. Through symmetry, it may also be seen that for the quad coil system of Figure 2.2, L1,3 = L1,4 = L2,3 = L2,4 and L1,2 = L3,4. It is obvious that for four identical coils L1,1 = L2,2 = L3,3 = L4,4. From these symmetry observations, it is evident that only three formulations are needed to summarize all of the self and mutual inductance values for the system in Figure 2.2. Calculation of L1,1, L3,1, and L2,1 will facilitate a full description of the system inductance. For the first treatment, the self inductance (M = N) will be examined and the indices M and N will be temporarily dropped. In a general sense, the total flux linked by the pth loop of a coil may be defined as the surface integral p = Z @Sp B · d ~ Sp, (3.23) where B is the total flux density vector field and @Sp is the surface outlined by the pth loop. For the single filamentary loop of Figure 3.1, the total flux linked may be rewritten as = Z @A BzdA, (3.24) where Bz is the component of the total flux density in the z direction linked by the loop and @A is the area outlined by the loop. One may suggest now that (3.24) be solved directly using the relation 42 B = μ0H =) Bz = μ0Hz. For the multiturn coil with q loops, the z component of the total flux density may be formulated from (3.6) and (3.7) as Bz = μ0Hz = jμ0Gz = jμ0 Xq p=1 gz,p. This may be used with (3.24) and (3.8) to formulate the total flux (produced by a coil having q = tl turns) linked by the same coil’s ˜pth loop as ˜p = jμ0 Xl =1 Xt =1 Zbp˜ −b˜p Zap˜ −a˜p gz,p(x, y, z˜p, a , b , c+ )dxdy, (3.25) where z˜p is the z location of the ˜pth loop in the local coordinate system of Figure 3.2. It must also be noted that for the self inductance formulation, a , b , and c+ obey the definitions provided in (3.8). Substitution of the relation for gz from (3.3) into the integrand of (3.25) makes for a demanding exercise in calculus. The evaluation of the total flux can be simplified by utilizing the magnetic vector potential. For the single filamentary loop of Figure 3.1, the magnetic vector potential may be formulated as ~A = j 4 I d~l r = jA, (3.26) where A may be thought of as a geometric vector potential. The magnetic vector potential is related to the magnetic field through the relation H = r×~A, therefore the flux density may be related to the geometric vector potential through B = jμ0(r×A). This may be substituted into (3.23), resulting in a formula describing the total flux linked by the pth loop in terms of total vector potential p = jμ0 Z @Sp (r × A) · d ~ Sp. (3.27) The integration of (3.27) is no less complicated than (3.23), however (3.27) may now be reduced from a surface integral to a path integral by use of Stokes’ theorem [33]. Application of this theorem results in the total flux linkage being defined as p = jμ0 I @sp A · d ~ sp, (3.28) where @sp is the closed path around the pth loop. Clearly, this is a more mathematically tractable problem. Before (3.28) may be used, the quantity A must be evaluated. Initially, the vector potential of the single filamentary loop, denoted by the vector field a, will be examined and then superposed with all loops comprising a coil to form A. The vector field a may be evaluated by using (3.26) and 43 a piecewise integration similar to the technique used in (3.2), a = 1 4 " Z b −b d~l 1 r1 + Z −a a d~l 2 r2 + Z −b b d~l 3 r3 + Z a −a d~l 4 r4 # . (3.29) Note that the vector field a exists only in the ˆi and ˆj directions. Evaluation of (3.29) in the rectangular coordinate system of Figure 3.1 yields a(x, y, z, a, b, c) = ax(x, y, z, a, b, c)ˆi + ay(x, y, z, a, b, c)ˆj , (3.30) where ax = 1 4 ln #2#3 #1#4 , ay = 1 4 ln #6#7 #5#8 , #1 = (x + a) + '7, #2 = (x − a) + '5, #3 = (x + a) + '8, #4 = (x − a) + '6, #5 = (y − b) + '5, #6 = (y + b) + '6, #7 = (y − b) + '7, #8 = (y + b) + '8. In a fashion similar to that of (3.8), the total geometric vector potential for the coil in Figure 3.2 with q turns consisting of l layers and t turns per layer is A = Xq p=1 ap = Xl =1 Xt =1 a(x, y, z, a , b , c+ ). (3.31) Equation (3.31) may be substituted into (3.28) to calculate the total flux linked by the path @sp. The flux linked by the ˜tth turn in the˜l th layer may now be ascertained through the piecewise integration ˜t,˜l = jμ0 Xl =1 Xt =1 " Z b˜l −b˜l ay(a˜l , y, c˜t, a , b , c+ )dy + Z −a˜l a˜l ax(x, b˜l , c˜t, a , b , c+ )dx+ Z −b˜l b˜l ay(−a˜l , y, c˜t, a , b , c+ )dy + Z a˜l −a˜l ax(x,−b˜l , c˜t, a , b , c+ )dx # . (3.32) Solution of (3.32) yields ˜t,˜l = jμ0 4 Xl =1 Xt =1 ˜t,˜l , , , (3.33) 44 where ˜t,˜l , , = b˜l ln Y8 m=5 #m(I1)#m(I3) #m(I2)#m(I4) + a˜l ln Y4 m=1 #m(I3)#m(I4) #m(I1)#m(I2) + b ln Y8 m=5 #m(I1)#m(I4) #m(I2)#m(I3) (−1)m + a ln Y4 m=1 #m(I1)#m(I4) #m(I2)#m(I3) (−1)m , I1 = (a˜l , b˜l , c+ ˜t , a , b , c+ ), I2 = (−a˜l , b˜l , c+ ˜t , a , b , c+ ), I3 = (a˜l ,−b˜l , c+ ˜t , a , b , c+ ), I4 = (−a˜l ,−b˜l , c+ ˜t , a , b , c+ ). Note that in (3.33), I denotes the variable assignments in the functions # and '. Additionally, a˜l = a , b˜l = b , and c+ ˜t = c+ , for ˜l = , and ˜t = . Use of (3.32) and (3.22) allows for the formulation of the self inductance L1,1 as L1,1 = 1 j ˜l X ˜ =1 Xt˜ ˜ =1 ˜ ˜ = μ0 4 ˜l X ˜ =1 Xt˜ ˜ =1 Xl =1 Xt =1 ˜ ,˜ , , . While (3.33) is convenient, it does have limitations. The integrands of (3.32) are undefined when ˜t = and ˜l = , and hence must be ignored when calculating (3.33), or in other words e ˜t,˜l = 8>>< >>: 0, for ˜t = and ˜l = ; jμ0 4 Pl =1 Pt =1 ˜t,˜l , , , for ˜t 6= or ˜l 6= . (3.34) A modified self inductance may now be formed which neglects the undefined integrands, eL 1,1 = 1 j ˜l X ˜ =1 Xt˜ ˜ =1 e ˜ ˜ . (3.35) This will result in error being introduced to the calculation, however if the coil windings are assumed to be closely packed, an estimate of the error may be obtained. Assuming that the coil windings are closely packed allows for the approximations A = a1 + a2 + · · · + aq qa, p = jμ0 I @sp A · d ~ sp qjμ0 I @s a · d~s. The total flux linked by the coil would approximately be = Xq p=1 p q2jμ0 I @s a · d~s. 45 If the indefinite path integral is neglected, then eA (q − 1)a, e p (q − 1)jμ0 I @s a · d~s. The total flux linked by the coil when neglecting the undefined integrands is then e = Xq p=1 e p q(q − 1)jμ0 I @s a · d~s. Using (3.22), the error estimate Es of the self inductance resulting from neglecting the undefined integrands may be formed as Es L −eL L = 1 q . As can be seen, the error in the self inductance calculation utilizing finite loops and neglecting the undefined integrands is proportional to the inverse of the number of turns in the coil. A corrected self inductance calculation may then be formed as bL 1,1 = eL 1,1 1 + 1 q . To calculate the mutual inductances L3,1 and L2,1, a similar approach may be used. Additionally, for the mutual inductance formulations it will not be necessary to neglect some of the integrations, as there will no longer be any overlapping integration paths and hence no undefined integrands. Performing a similar derivation to that of the self inductance formulation yields the mutual inductance of two adjacent coils as L3,1 = μ0 4 ˜l X ˜ =1 Xt˜ ˜ =1 Xl =1 Xt =1 ˜ ,˜ , , , (3.36) where ˜ ,˜ , , = (a˜ )3,1 ln Y4 m=1 #m(L3)#m(L4) #m(L1)#m(L2) + a ln Y4 m=1 #m(L1)#m(L4) #m(L2)#m(L3) (−1)m + X8 m=5 ['m(L4) − 'm(L3) + 'm(L1) − 'm(L2)], L1 = ((a˜ )3,1, (b˜ )3,1, (c˜ )− 3,1, a , b , c− ), L2 = (−(a˜ )3,1, (b˜ )3,1, (c˜ )− 3,1, a , b , c− ), L3 = ((a˜ )3,1,−(b˜ )3,1, (c˜ )+3 ,1, a , b , c− ), L4 = (−(a˜ )3,1,−(b˜ )3,1, (c˜ )+3 ,1, a , b , c− ), (a˜ )3,1 = a0 + d p 3 2 (˜ − 1), 46 (b˜ )3,1 = 8>>< >>: c0 + (˜ − 1)d; if ˜ is odd, c0 + d 2 + (˜ − 1)d; if ˜ is even, (c˜ )± 3,1 = c0 ± b0 + d p 3 2 (˜ − 1) , and the mutual inductance of two opposite coils as L2,1 = μ0 4 ˜l X ˜ =1 Xt˜ ˜ =1 Xl =1 Xt =1 ˜ ,˜ , , , (3.37) where ˜ ,˜ , , = b˜ ln Y8 m=5 #m(J2)#m(J4) #m(J1)#m(J3) + a˜ ln Y4 m=1 #m(J1)#m(J2) #m(J3)#m(J4) + b ln Y8 m=5 #m(J2)#m(J3) #m(J1)#m(J4) (−1)m + a ln Y4 m=1 #m(J2)#m(J3) #m(J1)#m(J4) (−1)m , J1 = (a˜ , b˜ , (c˜ )2,1, a , b , c− ), J2 = (−a˜ , b˜ , (c˜ )2,1, a , b , c− ), J3 = (a˜ ,−b˜ , (c˜ )2,1, a , b , c− ), J4 = (−a˜ −, b˜ , (c˜ )2,1, a , b , c− ), (c˜ )2,1 = 8>>< >>: 2c0 + [d(˜ − 1)], for ˜ = 1, 3, 5, . . . (odd); 2c0 + [ d 2 + d(˜ − 1)], for ˜ = 2, 4, 6, . . . (even). The coil resistance calculation is quite straightforward. If the total length of the wire in the coil is known, then this length may be multiplied by the resistance per length Rlength of the wire used to wind the coil. As in the previous field calculations and inductance formulations, the coil is approximated as a series of finite loops. This results in the calculation of resistance as R = 4tlRlength a0 + b0 + d p 3 2 (l − 1) . (3.38) Alternatively, coil resistance may be calculated on a weight W basis, where the parameter of resistance per unit weight Rwgt is utilized to calculate the resistance R = RwgtW. While (3.38) is a useful formulation for use in a predictive model or optimization routine, the weight based calculation is much more practical when winding the coil. 47 CHAPTER 4 Design and Optimization of Experiment 4.1 Functional for Minimization and Algorithm In a system such as the one proposed in Figure 2.2, there are many design parameters which must be determined before fabrication of an experimental setup can proceed. It may be assumed for the purposes of this study that certain aspects of the system are “fixed” and are not considered variables to be changed or adjusted in any kind of optimization procedure. The fixed characteristics are listed as follows: sphere radius, sphere magnetization, fluid viscosity, and sphere 2D trajectory. Additionally, a large amount of 12 gauge heavy insulation magnet wire was already in hand, resulting in wire diameter being fixed as well. In contrast, certain parameters of this system will be used as variables in an optimization process. These parameters are given in Table 4.1. Upon first observation, one may reason that the coil resistance should be among the variables listed in Table 4.1, however this is not the case. Fabrication of the experimental apparatus will utilize a magnet wire already in hand (12 gauge), thus establishing the wire diameter. The coil resistance is a function of wire diameter, electrical conductivity of the wire utilized, and length of the wire wound into coil form. Of these three criteria, only the wire length will be varied during the optimization and it is a function of the first four parameters of Table 4.1. In devising a “bench top” test apparatus, it is desirable to have a hardware setup which will move the sphere with the least amount of effort. Hence a minimization of some sort is in order. From the standpoint of the physical system, the attribute which must be minimized is the power applied to the electromagnet coils. The instantaneous power supplied to any reactive electrical component is defined as P = vj, where v is the potential across the device and j is current flowing in the device [34]. Adopting the approach used in [32] and representing the coil as an RL circuit, the instantaneous power in the mth coil may be stated as Pm = Rmj2m + jm XQ n=1 Lm,n djn dt , (4.1) 48 Table 4.1: System Parameters to be Used for Optimization parameter symbol unit number of wire turns/layer t    number of wire layers l    coil inner “x” dimension ain inch coil inner “y” dimension bin inch distance of coil face from global origin c0 inch where Rm and jm are the resistance and current for the mth coil, Lm,n is the mutual (m = n) or self (m 6= n) inductance of the mth coil on the nth coil, and Q is the total number of coils. Consistent with Chapter 2, the optimization problem will be based the ability of the hardware system to move a sphere about a circular orbit. Hence, symmetry may be leveraged and only coils 1 and 3 examined. To calculate the instantaneous power of the coil combination depicted in Figure 2.10, denoted here as 1,3 = (P1,P3), the formulations for coil resistance and inductance (both self and mutual) developed in Chapter 3 may be employed. The analytical evaluation of the time derivative of current would prove to be quite tedious, if even possible, due to the required formulation of the partial derivatives of the vector field term ~kFM. A complete mathematical description of the vector field term is available in Chapter 2. Because of this, the time derivative of the required currents in coils 1 and 3 will be calculated numerically from discrete points comprising the inverse current solution of a given system configuration. Obviously in a reactive system such as an RL ciruit, achieving the instantaneous current change at the beginning and end of the inverse current solution as depicted in Figure 2.10 would be impossible. Indeed an inductive “switching effect” would be inherent in the proposed technique, however for the purposes of optimizing an experimental design this will be neglected. For the set of discrete points making up the inverse current solution J ,1,3, where = a, b, c, d, a forward difference approximation was used to calculate the derivative at the initial point while a backward difference approximation was used to calculate the derivative at the last point. For all other points making up the set, a central difference technique was used. This will provide a reasonable approximation of the function J˙ ,1,3, provided that the time step spacing is small [35]. The time derivative of current was calculated using points spaced at 10 msec intervals. An 49 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ a /dt (A/sec) dj 1 /dt dj 3 /dt 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ b /dt (A/sec) 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ c /dt (A/sec) 0 1 2 3 4 & 6 ( 8 * 10 !10 0 10 dJ d /dt (A/sec) time (seconds) Figure 4.1: Calculated current derivative from inverse solution of coil combination 1&3. Displayed results correspond to system with R1,3 = 0.547 , L1,1 = 6.79mH, and L1,3 = 0.58mH (calculated). example calculation on the current trajectory depicted in Figure 2.10 is shown in Figure 4.1. It should be noted that the trace shown in Figure 2.10 corresponds to a system where the resistance, inductance and mutual inductance have been calculated to be R1,3 = 0.547 , L1,1 = 6.79mH, and L1,3 = 0.58mH. Utilizing equation (4.1) with the discrete traces of j1,3 and dj1,3/dt along with the calculated values for resistance and inductance allows for the estimation of the instantaneous power required by the coil combination which results in the inverse current solution. The instantaneous power calculation for this configuration is shown in Figure 4.2. As can be seen, the instantaneous power has two distinct forms rather than four. This is indicative of the same power being needed for the solutions corresponding to opposite polarities. With the proposed minimal coil set technique, two coils will be powered at any one given time. Hence, it would be advantageous to utilize an expression as the functional for minimization which incorporates the power applied to both members of the coil set. The Euclidean norm was employed to satisfy this convenience and is defined here as k m,nk = p P2m + P2n , where m, n = 1, 2, 3, 4. Obvioulsy k m,nk = k n,mk and for the case being examined m = 1 and n = 3. The norms of the two instantaneous power traces are shown in Figure 4.3. From 50 0 1 2 3 4 5 6 7 8 9 10 0 2 4 !a (Watts) P 1 P 3 0 1 2 3 4 5 6 7 8 9 10 0 2 4 !b (Watts) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 !c (Watts) 0 1 2 3 4 5 6 7 8 9 10 !1 0 1 ! d (Watts) time (seconds) Figure 4.2: Instantaneous power required for inverse current solution of coil combination 1&3. 0 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 time (seconds) !1,3  (Watts) set "a" set "c" Figure 4.3: Plot of k 1,3k for solution sets a and c. The norm of solution sets b and d are identical to sets a and c, respectively. 51 the power norm, a single metric may be defined to describe the performance of a coil system with the parameters seen in Table 4.1. If the maximum value of the current norm, defined here as M= max k m,nk, is used as criteria for minimization, then the minimized configuration Smin would be the set of parameters S = {t, l, a0, b0, c0} which results in the minimum M when compared to all parameter values possible. More precisely, {t, l, ain, bin, c0} = S = Smin if 9{t, l, ain, bin, c0} : k m,n(t, l, ain, bin, c0)k inf{M(S)}8S 2 R. In an effort to reduce the number of parameters utilized in the search routine, it was decided to make c0, the distance from the face of the coils to the origin (see Figure 2.2), a function of the parameters l and bin. To achieve this, a spacing criteria csp is defined as csp = bin + d 2 + (l − 1) p 3 2 + sep, (4.2) where d is the wound wire diameter, bin is the inner Z dimension of the rectangular coil, sep is a separation constant, and l is the number of layers of wire wound on the coil. It should be noted that the formulation of (4.2) is based upon a coil having a “perfect” winding pattern as defined in Chapter 3. The minimum coil to coil spacing for this experiment was chosen to be 5”, resulting in a minimum coil spacing from the origin of cmin = 2.5”. In addition, the separation constant in (4.2) was chosen to be sep = 0.25”. To reduce the number of parameters for the optimization problem, c0 may then be defined as c0 = 8>>< >>: cmin if csp cmin, csp otherwise. This relationship sets the coil spacing as 2.5 inches unless the coil physical dimensions do not allow for them to be positioned this closely. In that case, the coil spacing from the origin is dictated by (4.2). The technique selected to find the minimal parameter set is a basic gradient method: steepest descent [36, 37]. More specifically, a discrete steepest descent algorithm was used because the system model was evaluated at specific points in the parameter space. This work will utilize a rather straigtforward version of this gradient method, which will be presented in two parts. The first component may be summarized as a method by which, for a given ain and bin, the number of turns per layer t and the number of layers l resulting in the minimumMis found. To achieve this, we may look at M as a function of the minimizing variables t, l, ain, and bin, hence M = M(t, l, ain, bin). For a given initial guess (ti, li) and a given (ain, bin), the initial value of M = M(ti, li, ain, bin) is calculated. The points “around” the initial guess are utilized to then calculate the values of M around the initial guess. These points make up the ordered set Si u = {(ti+1, li), (ti+1, li−1), (ti, li− 52 yes no Figure 4.4: A flowchart illustrating the algorithm used to find the number of turns tf and number of layers lf which results in minimal power used for a given ain and bin. 1), (ti − 1, li − 1), (ti − 1, li), (ti − 1, li + 1), (ti, li + 1), (ti + 1, li + 1)}, where u = 1, 2, . . . , 8 denotes a particular ordered set member. This set may then be used to calculate the values of M at each set point, {M(Si u, ain, bin)} = {M(Si 1, ain, bin), . . . ,M(Si 8, ain, bin)}. The “next step” (ti+1, li+1) may then be determined as the vth element of Si u which results in G = min{M(Si u, ain, bin) − M(ti, li, ain, bin)} and G 0. This process may be repeated until a point (tf , lf ) is found such that G > 0. At which time, a minimizing t and l have been found. This process is illustrated in the flow chart shown in Figure 4.4. The second component of the complete minimization algorithm may be viewed as an analogy of the first component discussed in the previous paragraph and shown in Figure 4.4. It may be summarized as a method by which, for a given t and l, the inner x coil dimension ain and the inner y coil dimension bin resulting in the minimumMis found. As before, this may be achieved by regarding M as a function of the variables which are being “varied”. For a given initial guess (ain,j , bin,j) and a known value (t, l), the initial value ofM=M(t, l, ain,j , bin,j) is calculated. Once again, the points “around” the initial guess are utilized to calculate the values of M around the initial guess. These points make up the ordered set Tj u = {(ain,j + , bin,j), (ain,j + , bin,j − ), (ain,j , bin,j − ), (ain,j − , bin,j − ), (ain,j − , bin,j), (ain,j − , bin,j + ), (ain,j , bin,j + ), (ain,j + , bin,j + )}, where u = 53 yes no Figure 4.5: A flowchart illustrating the algorithm used to find the inner x coil dimension ain,f and inner y coil dimension bin,f which results in minimal power used for a given t, l, and . 1, 2, . . . , 8 denotes a particular ordered set member and denotes a search length. This set may then be used to calculate the values ofMat each set point, {M(t, l,Tj u)} = {M(t, l,Tj 1), . . . ,M(t, l,Tj 8)}. The “next step” (ain,j+1, bin,j+1) may then be determined as the wth element of Tj u which results in G = min{M(t, l,Tj u) −M(t, l, ain,j , bin,j)} and G 0. This process may be repeated until a point (ain,f , bin,f ) is found such that G > 0. At which time, a minimizing ain and bin have been found. This process is illustrated in the flow chart shown in Figure 4.5. As can be seen, Figures 4.4 and 4.5 are similar in structure. To achieve a single algorithm which finds a minimizing point (tf , lf , ain,f , bin,f ), one algorithm will be “nested” in the other. For this study, it was decided to nest the algorithm depicted in Figure 4.4 into the algorithm depicted in Figure 4.5. The flow of the complete process is similar to that of the procedures previously outlined. For an initial guess (ti, li, ain,j , bin,j ), the initial value of M = M(ti, li, ain,j , bin,j) must be calculated. Subsequently, M = M(t , l , ain,j , bin,j) which has been minimized on (t, l) may be calculated using a nested algorithm similar to that in Figure 4.4. In this case, represents an arbitrary number of iterations required for convergence. Similarly, the minimal {(t , l )u} for each set member of Tj u may be calculated as well by a nested (t, l) minimization algorithm, where represents an arbitrary number of iterations required for convergence. These results are then used to check 54 yes no Figure 4.6: A flowchart illustrating the complete algorithm used to find the parameters (tf , lf , ain,f , bin,f ) resulting in minimal power. The nested algorithms are indicated by “broken line” boxes. for a descent direction in the variable space. This procedure is then repeated until a minimizing configuration (tf , lf , ain,f , bin,f ) is found. The complete algorithm is shown in Figure 4.6. It must be noted that while not explicitly depicted in Figure 4.6, a value for the search length must be selected before using the algorithm. The final aspect to be resolved before the minimization routine may be implemented is which norm set will be used (see Figure 4.3). For this study, “a” was selected as the set utilized in the functional evaluation. With the minimization routine developed, and the system model in place, a minimizing parameter set may be calculated. 4.2 Minimization Results and Experimental Setup A minimizing parameter set (tf , lf , ain,f , bin,f ) was indeed found from utilizing the routines developed in Section 4.1. Three different search efforts were performed with search lengths valued as = 0.1”, 0.01”, and 0.001”. This method may be viewed as one in which an initial “coarse” search is performed to establish a general area of the extremum and then each successive search utilizes a finer search grid. The optimized design parameters for the quad coil assembly are given in Table 4.2. As can be seen, the optimal coil design corresponds to one in which the coil faces are placed as close 55 Table 4.2: Optimized System Parameters parameter symbol value unit number of wire turns/layer t 18    number of wire layers l 21    coil inner “x” dimension ain 0.439 inch coil inner “y” dimension bin 0.757 inch distance of coil face from global origin c0 2.501 inch to the origin as allowed, in this case 2.5”. This is indicative of the proximity of the particle to the electromagnet as being a major contributor to the amount of power required for particle movement. A coil assembly along with the tray providing the plane of motion for the sphere was constructed to these specifications and is shown in Figure 4.10. As expected, the geometry of the wound coi 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


