UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 COUPLED MAGNETOELASTOSTATIC AN ALYSIS USING IMPLICIT BOUNDARY FINITE ELEMENT METHOD By SUNGUK ZHANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORID A IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010 PAGE 2 2 2010 SungUk Zhang PAGE 3 3 To God and my family PAGE 4 4 ACKNOWLEDGMENTS I would like to express my sincere gratit ude to my advisor, Dr. Ashok V. Kumar for his guidance, enthusiasm, and constant suppor t throughout my doctoral research. His recommendations and suggestions have been inval uable for the research. I would also like to thank the members of my advisory committee, Dr. Bhav ani V. Sankar, Dr. John K. Schueller, Dr. Peter G. Ifju, and Dr. Toshi Nishi da. I am grateful fo r their willingness to serve on my committee, for valuable suggesti ons in my oral qualifying examination and for reviewing this dissertation. I also wish to thank my colleagues at Design and Rapid Protot yping Laboratory at the University of Florida for their help and thei r dissertations. The dissertations of Ravi K. Burla and Sanjeev Padmanabhan we re valuable. The help of Prem Dheepak and Nitin Chandola improved my research. I would like to thank my family. Without their constant love and support, this woul d not have been possible. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDG MENTS..................................................................................................4LIST OF TABLES............................................................................................................8LIST OF FI GURES..........................................................................................................9ABSTRACT ...................................................................................................................14 CHAPTER 1 INTRODUC TION....................................................................................................16Overvi ew.................................................................................................................16Goals and Obje ctives..............................................................................................21Outlin e....................................................................................................................222 MAGNETOSTATI C ANALSI S.................................................................................24Computational El ectromagnet ics............................................................................24Maxwells E quations ...............................................................................................25Overview of Magnet ic Act uator...............................................................................27Reluctance Method.................................................................................................303 IMPLICIT BOUNDARY FINI TE ELEMENT METHOD.............................................33Overview of Finite Element Method........................................................................33IBFEM for Elas tostatics...........................................................................................35Finite Element Formulati on for Elasto statics...........................................................38Bspline el ement .....................................................................................................39Timoshenko Beam Exampl e...................................................................................424 IMPLICIT BOUNDARY FINITE ELEMENT METHOD FOR 2D MAGNETOSTA TICS...............................................................................................44Solution Structure fo r 2D Magnetos tatics................................................................44Finite Element Formulation for 2D Magnet ostatics.................................................46Solution Structure for 2D Magnet ostatics with Perma nent M agnet.........................49Solution Structure for Multiple Ma terials.................................................................51Magnetic Force Co mputatio n..................................................................................585 IMPLICIT BOUNDARY FINITE ELEMENT METHOD FOR 3D MAGNETOSTA TICS...............................................................................................62Governing Equation and Weak Fo rm for 3D Magnet ostati cs..................................62 PAGE 6 6 Solution Structure fo r 3D Magnetos tatics................................................................64Current Density Computat ion..................................................................................676 OPEN BOUNDARY TECH IQUES FOR IBFEM......................................................68Overview of Open Boundary Techniques fo r IBFEM..............................................68Asymptotic Boundary Cond itions for IBFEM...........................................................70Derivation of Asymptot ic Boundary C onditions.................................................71Two Dimensional Asymptot ic Boundary Co nditions .........................................72Infinite Element s for IB FEM....................................................................................73Decay Interpolation Function A pproach...........................................................74Decay Function Infinite Element for IBFEM......................................................787 RESULTS AND DI SCUSSION S.............................................................................81Two Dimensional Magnetos tatic Probl ems.............................................................81Example 711: 2D Coaxial C able....................................................................81Example 712: 2D Clapper Solenoid Actuator.................................................85Example 713: 2D Clapper Solenoid Actuator with Artificial Damage.............87Example 714: 2D Swit ched Reluctanc e Moto r...............................................88Three Dimensional Magnet ostatic Pr oblems..........................................................91Example 721: Iron Block in a Homogenous Magnet ic Fiel d...........................92Example 722: 3D Coaxial C able....................................................................96Example 723: 3D Plunger Solenoid Actuator.................................................99Example 724: 3D Clapper Solenoid Ac tuator...............................................104Magnetostatic Problems wit h Permanent M agnets...............................................109Example 731: US haped Permanent Magnet...............................................109Example 732: Three Dimensio nal Cylindrical Magnet ..................................113Magnetostatic Problems with Op en Boundary Tec hniques...................................115Example 741: UShaped Perm anent Magnet with Open Boundary Techniques..................................................................................................115Example 742: 2D Clapper Sole noid Actuator with Open Boundary Techniques..................................................................................................118Example 743: Force Between Two Parallel Wires.......................................120Coupled MagnetoElastost atic Probl ems..............................................................123Example 751: 2D Clapper Solenoid Actuator with Cant ilever Beam............124Example 752: 3D Plunger Solenoid Actuator with Structures......................1268 DESIGN AND ANALYSIS OF MAGNET IC ACTUATORS US ING IBFEM............131Design Crit eria...................................................................................................... 131Solenoid Actuators with Clapper Arma ture...........................................................132Solenoid Actuators with Plunger Arma ture...........................................................137Solenoid Actuators with a Combined Plunger & Clapper Armature......................142Coil Actuat ors.......................................................................................................144The Best Actuator among t he Designed Act uators...............................................146Coupled MagnetoElastostatic Probl ems with a Flapping Wing Mo del.................149 PAGE 7 7 9 CONCLUSIONS AND FU TURE WORK...............................................................155Conclusion s..........................................................................................................155Future Wo rk..........................................................................................................157LIST OF REFE RENCES.............................................................................................159BIOGRAPHICAL SKETCH ..........................................................................................164 PAGE 8 8 LIST OF TABLES Table page 61 Comparison of different open boundary techni ques [40]...................................6881 Comparison for three clapper solenoid actuators ( NI =30)...............................13782 Comparison for three plunger solenoid actuators ( NI =30)...............................14183 Comparison for three combined plunger & clapper solenoid actuators ( NI =30)............................................................................................................. 14384 Comparison for coil actuators ( NI =30).............................................................146 PAGE 9 9 LIST OF FIGURES Figure page 11 Confo rmal me sh................................................................................................1612 Conformal mesh ve rsus.....................................................................................1813 3D structured gird s on multima terial .................................................................1814 Flapping wings operated by magnetic actuator.................................................2121 Maxwells equati ons..........................................................................................2522 Block diagram of a magnetic ac tuator...............................................................2823 Reluctance method (Magnet ic circuit method)...................................................3124 Fringi ng flux .......................................................................................................3231 The solution structure with t he essential boundar y condit ion............................3732 One dimensional field so lution using tw o element s...........................................4033 Displacement surfac e plots in ycom ponent...................................................... 4341 2D magnetosta tic prob lem.................................................................................4442 2D magnetostatic prob lem with perm anent m agnet.......................................... 5043 Two grids for two materi als................................................................................5144 Component plots of interface soluti on structure at the material boundary.........5445 Interface solution struct ure at the mate rial boundar y.........................................5546 Three parts with own gr ids.................................................................................5647 Contac t pairs .....................................................................................................5648 Four element s in Part 1.....................................................................................5749 Material in terface boundary ...............................................................................5951 Analysis doma in and boundar ies.......................................................................6352 Sequential analysis fo r 3D Magnetos tatics........................................................6761 Twodimensional interior region with an exte rior annul us..................................69 PAGE 10 10 62 Structured grid and rectangular boundary.........................................................7263 Coordinate transformation for the decay function infi nite element.....................7564 1D shape functions using exponential decay functi ons.....................................7765 Structured grids with the infinite elements.........................................................7866 Domain of integrat ion........................................................................................7967 Solution plots for 1D in finite element example...................................................8071 2D coaxial cable model with structur ed grid s.....................................................8272 Magnitude of H fi eld...........................................................................................8273 Convergence pl ot for H1 norm...........................................................................8374 Convergence pl ot for L2 norm...........................................................................8475 2D Clapper solenoi d with struct ured gr id...........................................................8576 Magnetic vector potentia l and magnetic fl ux li nes.............................................8677 Magnetic force ve rsus air gap length.................................................................8778 Clapper solenoid wit h artificial damage.............................................................8879 2D Planar model of s witched reluctanc e motor.................................................88710 Magnetic vector potent ial and magnetic flux lines in the aligned position..........90711 Magnetic vector potent ial and magnetic flux lines in the unaligned position......91712 Iron cube in homogeneous magnetic field.........................................................92713 Iron objects with t he same grid density..............................................................93714 Crosssections wit h the line y= z=10mm............................................................94715 Components of B al ong the line y= z=10mm......................................................95716 3D coaxial cable model with the stru ctured gr id................................................96717 Magnitude of H field for 3D coaxia l cabl e..........................................................98718 Magnetic field in the hood direction vers us. radi us.............................................98719 Convergence pl ot for H1 norm...........................................................................99 PAGE 11 11 720 Plunger solenoid actuator of axisymmetric geometry......................................100721 Top view of two plunger actuat ors...................................................................101722 3D solid model of solenoid ac tuators with plunger armature s..........................101723 The magnetic flux density in the ydirection for two pl unger solenoids............102724 The magnetic field in the ydire ction for two plung er solenoi ds.......................103725 Magnetic force versus gap length for plunger solenoids ..................................104726 Clapper solenoid actuator of axisymmetric geometry......................................105727 Top view of two clapper actuat ors...................................................................105728 3D solid model of solenoid ac tuator with clapper armatures...........................106729 The magnetic flux density in the ydirection for two cl apper actuators.............107730 The magnetic field in the ydire ction for two clapper actuators........................107731 Magnetic force versus gap length for clapper solenoids ..................................108732 Ushaped permanent magnet model ...............................................................109733 Surface and contour plots for magnetic potential fr om Comsol.......................110734 Surface and contour plots for magnetic potential fr om IBFEM........................111735 The observation line for the com parison ..........................................................111736 Magnetic field in the xdirection on the line ......................................................112737 Magnetic field no rm on the lin e.........................................................................112738 3D solid model for cylindr ical magnet in the air................................................113739 Magnetic field in the ydirect ion........................................................................114740 Magnetic field along y ax is...............................................................................115741 Ushaped permanent magnet with the reduc ed air domai n..............................116742 Contour plots for magnet ic vector pot ential......................................................117743 Magnetic field in the xdire ction on the obser vation lin e...................................117744 2D Clapper sol enoid actuat or...........................................................................118 PAGE 12 12 745 Contour plots of magnetic vector potentia l.......................................................119746 Magnetic force ve rsus gap l ength.....................................................................120747 Two long para llel wires .....................................................................................121748 Contour plots for magnet ic vector pot ential......................................................122749 The magnetic force versus t he distance between two wires .............................123750 Planar clapper solenoid act uator with a cant ilever bea m..................................124751 Displacement at the tip versus NI with varying attach ment location.................125752 Planar clapper solenoid act uator with a cant ilever beam ..................................125753 Displacement on the tip versus NI with varying beam thickness......................126754 Top views of structures atta ched to the plunger armature................................127755 3D plunger solenoid act uators with A) Solid plate, B) plate with one hole, and C) plate with two holes .....................................................................................128756 Deformation due to magnetic fo rce...................................................................129757 Maximum displacement versus NI.................................................................... 13081 Solenoid actuator wit h clapper arma ture..........................................................13382 Computed results using IBFE M........................................................................13583 Clapper solenoi d actuator s...............................................................................13584 Contour plots of magnetic vector potential for cla pper solenoid actuators........13685 Solenoid actuator wit h plunger arma ture..........................................................13786 Computed results using IBFE M........................................................................13987 Plunger solenoi d actuator s...............................................................................14088 Contour plots of magnetic vector potential for pl unger solenoid actuators.......14189 Three combined plunger & cl apper solenoid ac tuators.....................................142810 Contour plots of magnetic vector potential fo r combined plunger & clapper solenoid act uators............................................................................................143811 Three coil ac tuators..........................................................................................144 PAGE 13 13 812 Magnitude of B fields for three coil ac tuators....................................................145813 The best magnetic actuator among several designed actuators.......................146814 Sold model of the 3D coil act uator.................................................................... 147816 The magnitude of the magnetic flux density of the coil act uator.......................148817 Lorentz force of the coil actuator versus NI...................................................... 149818 The coil actuator with flapping wings................................................................149819 Four structures with flapping wings..................................................................150820 Four surface structur es with flappi ng wings ......................................................151821 Structured mesh and boundary conditions of t he first de sign...........................152822 Displacement in the zdirect ion during the wi ng stroke .....................................153823 Tip displacement versus NI for wing upstr oke..................................................154 PAGE 14 14 Abstract of Dissertation Pr esented to the Graduate School of the University of Florida in Partial Fulf illment of the Requirements for t he Degree of Doctor of Philosophy COUPLED MAGNETOELESTOSTATIC AN ALYSIS USING IMPLICIT BOUNDARY FINITE ELEMENT METHOD By SUNGUK ZHANG May 2010 Chair: Ashok V. Kumar Major: Mechanical Engineering Implicit boundary finite element method (IBFEM) uses solution structures constructed using step functions to enforce boundary and interface conditions so that a structured grid can be used to perform the analysis. A structured grid, which consists of regular shaped elements, is much easie r to generate than conforming mesh thus eliminating the difficulties associated with mesh generation for complex assembles. In this study, IBFEM is extended to solve 2D and 3D magnetostatics, compute magnetic forces and to solve coupled magnet oelastostatic problems that typically involve an assembly of parts made of seve ral different material s. The geometry is accurately modeled using equations from CAD models and a separate structured mesh is used for each part in an assembly. Speciall y constructed solution structures are used to represent test and trial functions such that boundary and interf ace conditions are enforced. Several magnetostatic problems with known solutions are modeled to validate the method. The magnetostatic problems are cl assified as unbound problems so that sometimes a very large analysis domain sh ould be modeled to get more accurate results. In order to r educe the analysis domain, two open boundary techniques are PAGE 15 15 developed for IBFEM: asymptotic boundary conditions and decay function infinite element. In addition, a magnetostatic probl em with permanent magnets is solved using IBFEM. PAGE 16 16 CHAPTER 1 INTRODUCTION Overview Implicit boundary finite element method (I BFEM) is a modified finite element method that avoids the need for a confo rming mesh. IBFEM has been applied to solid mechanics and heat transfer problems in past work [13][16], [65]. In this study, IBFEM is extended to perform magnet ostatic analysis as well as c oupled magnetoelastostatic analysis. Magnetostatic analysis and force computati on for magnetic actuators involves modeling an assembly of components with diffe rent material properties. The traditional finite element method has been used for such analysis but it requires a conforming mesh that approximates the geometry of the assembly. The mesh must contain nodes along the external boundaries and the interf aces between parts. The edges / faces of the elements must approximate these inte rfaces. Figure 11 shows a typical mesh where both the boundary and the in terfaces between material s are approximated poorly if the mesh is not very dense. Furthermore the mesh generator has to ensure sufficient node density along the boundary / interface and pr oper connectivity so that element edges are not connected ac ross the interface. Figure 11. Conformal mesh. A) 33 elements, and B) 82 element PAGE 17 17 Generating such a mesh is difficult and, despite decades of research, 3D mesh generation (especially using hexahedral element s) is still not a fully automated process and in fact requires significant user inpu t. To address mesh generation difficulties, several meshless methods [1] have been proposed that still need a wellplaced distribution of nodes but does not requires these nodes to be connected into elements. Some of these methods have been successfully used for magnetostatic analysis [2][7]. These methods use interpolation and approximation schemes that do not need connectivity between nodes. However, com putationally these methods are more expensive and they still approximate boundaries and interfaces using nodes along them. An alternate approach to avoid mesh generati on difficulties is to use a structured background mesh to represent the solution wh ile using accurate equations of curves and surfaces to present the boundaries. A stru ctured mesh consists of uniform regular shaped elements and is t herefore easy to generate. Ext ended finite element method (XFEM) [810] is one such method which uses a structured mesh and implicit equations for the boundaries and interfaces. In the XFEM approach, the solution is enriched near singularities and discontinuities such as cra cks. An important application of this method has been fracture mechanics, where crack propagation [11][12] is simulated by modifying the equations of the crack rat her than regenerating t he mesh. Boundary and interface conditions have been imposed using Lagrange multiplier and Penalty methods for XFEM. The Implicit Boundary FEM (IBFEM), described in this study, uses solution structures constructed utilizi ng implicit equations of the boundaries to enforce boundary and interface conditions. This method has been applied to 2D and 3D elastostatics and PAGE 18 18 steady state heat transfer problems [13][ 16]. Structured mesh, which has uniform and undistorted elements, can be used for t he analysis because the implicit boundary method does not require nodes on the boundary to impose boundary conditions. Compared to the conformal mesh shown in Fi gure 12 B, the structured mesh, such as the examples shown in Figure 12 C, is ea sy to generate since all elements are regular shaped and the grid does not have to conform to the geometry. Figure 12. Conformal mesh versus Struct ured mesh. A) Solid model, B) Conformal mesh, and C) Structured mesh For modeling multiple materials and assembli es, a separate grid is generated for each material or part as shown in Figure 13. Wi thin overlapping elements at the interface, the piecewise interpolation wit hin each grid in combined into a single solution structure as explained in the later chapter. Figure 13. 3D structured girds on multimate rial. A) Material 1, B) Material 2 and C) Multimaterial PAGE 19 19 IBFEM is extended to solve magnetost atic problems. Some magnetostatic problems can be characterized as 2D problems the geometry has cons tant thickness or depth. If the magnetic fiel d can be assumed to be in a plane then the problem can be modeled as 2D planar magnetostatic problem s. Cylindrical shaped models and models symmetric about an axis are classified as t he 2D axisymmetric problem. Such simple structures can be analyzed under the 2D assumption. Magnetic vector potential becomes a scalar in 2D problems whose gov erning equation is a Poissons equation. The Poissons equation is also the governing equation in steady st ate heat transfer problems [14]. Under the 2D assumption, se veral applications with multiple components can be analyzed such as magnetic actuators, coaxial cables and switched reluctance motors. Using the IBFEM, it is possible to use Bspline shape functions instead of the traditional finite element shape functions if needed. The advantage of using Bsplines is that the computed magnetic flux density and field (or stresses and strains) are continuous between elements. Under the 2D assumption, magnetostatic a nalysis is computationally inexpensive; however, 3D magnetostatic analysis is necessa ry when the shape of the structures is not simple. In the case of 3D problems, the number of equations drastically increases because the magnetic vector potential is a 3D vector and the number of nodes per element increases due to the usage of 3D elements. Although the current density is a scalar value in 2D magnetostatic analysis, in 3D, the current density becomes a vector field. Therefore, the curr ent density distribution must be computed by electrostatic analysis prior to magnetostat ic analysis. Under the above considerations, several 3D magnetostatic problems are solved in this study. PAGE 20 20 Since magnetostatic problems are oft en infinite domain problems, an open boundary technique is needed to obtain more accu rate results in order to use small finite domain. Among several open boundary techniques, asymptotic boundary condition and infinite element with decay f unctions are implemented for IBFEM. Several structures including permanent magnet ar e studied and the results are compared to ones from commercial software or analytical solutions. After validating IBFEM for magnetostatic anal ysis, the method is extended to solve coupled magnetoelastostatic problems. One of the applications for the coupled analysis is for micro air vehicles. Micro air vehi cles (MAV) have recently been developed for military or scientific purpos es. MAV are useful for scouting in dangerous or hazardous area where ground vehicles cannot go. Airp lanelike fixed wings and helicopterlike rotary wings are widely used because of hi gher efficiency in the fixed wings and hovering capability in the rotary wings. One way to miniaturize MAVs further is to use flapping wings. In order to design the flappi ng wings actuated using magnetic forces, a coupled magnetoelastostatic analysis is needed because the magnetic force produced by the actuator deforms the flapping wings In traditional magnet ic actuators, the magnetic force is used to create a rigid body motion of a rotor and any attached mechanism. Coupled magnetoelastic analysis is needed when the magnetic forces produce structural deformation. In this research, using coupled analysis, magnetic actuators for flapping wings are designed as shown in Figure 14. Figure 14 A shows the downward wing stroke when the magnetic actuator turns off. Figure 14 B shows the upward wing stroke when it PAGE 21 21 turns on. Several magnetic actuators are char acterized in terms of magnetic force and iron weight to design a proper actuator for flapping wings. Figure 14. Flapping wings oper ated by magnetic actuator. A) Downward wing stroke, and B) Upward wing stroke Goals and Objectives The goal of this research is to extend t he Implicit Boundary Finite Element Method (IBFEM) to develop capability for multima terial system analysis and to perform magnetostatic analysis as well as c oupled magnetoelasto static analysis. The main objectives of this disse rtation are to extend IBFEM to Perform 2D magnetostatic analysis Perform 3D magnetostatic analysis Model permanent magnets PAGE 22 22 Perform multimaterial analysis Compute magnetic forces, and Solve coupled magnetoelastostatic analysis In addition, IBFEM are tested and validated us ing several examples: coaxial cable, solenoid actuator, Ushaped permanent magnet and so on. Based on the analysis capability of IBFEM, it is used as magnetic actuator design tool in order to design magnetically actuated structures. Outline The rest of the dissertati on is organized as follows: Chapter 2 briefly discusses computational electromagnetics, Maxwells equations, overview of magnetic actuators and ov erview of reluctance method. Among the Maxwells equations, Gausss law and Am peres law are used as the governing equations for elastostatics and magnetostatics. In addition, the chapter describes how to analyze magnetic actuators using the reluctance method. Chapter 3 discusses previous works for Im plicit Boundary Finite Element Method. The chapter briefly describes the motivation for IBFEM, applications for solid mechanics, and Bspline element. Chapter 4 discusses 2D magnet ostatic analysis with IBFEM. A modified solution structure is introduced for t he multimaterial analysis. Application of permanent magnet is described. In Chapter 5, 3D magnetostatic analysis with IBFEM is described. The chapter introduces the difficulties for 3D magnetostatic analysis, and 3D coupled magnetoelastostatic analysis is introduced. PAGE 23 23 Chapter 6 discusses open boundary techni ques. Infinite element with decay function and asymptotic boundary condition are introduced for IBFEM. Chapter 7 provides several results and di scussions for 2D and 3D magnetostatics, coupled analysis and perm anent magnet problems. Chapter 8 explains several small actuat ors for flapping wings under the specific specifications. IBFEM is introduced as a design tool. Finally, chapter 9 provides the summary of the results and conclusions. The further work is suggested in the end of the dissertation. PAGE 24 24 CHAPTER 2 MAGNETOSTATIC ANALSIS Computational Electromagnetics Researchers have studied electromagnetic devices using Maxwells equations. By the Maxwells equations, electromagnetic fiel ds are computed to characterize the behaviors of electromagnetic devices such as transistors, electrical machines, waveguides, and so on. Before the advent of computers, the onl y way to solve the Maxwells equations was using elaborate ma thematics such as series expansions, Legendre polynomials, Bessels functions, and so on. Using these methods, the solving process took days or they needed to make drastic simplifying assumption on device geometry, current and charge dist ribution. For example, the geometry is assumed to be circular or rectangular and t he current distribution is suggested to be uniform in the analysis domain. Under those assumptions, t hey can obtain a closed form solution. After the advent of the com puter, several numerical schemes could make solutions without the severe assumptions. With a r ealistic design, researchers can solve Maxwells equations so that they can pr edict the behavior of the electromagnetic devices. Although those schemes are approx imation methods, solutions by those schemes are accurate enough. Sometimes using the simple geometry such as circular or rectangular, the analytic solution by cla ssical methods is more accurate than the numerical solution. However, the numerical schemes are not comparable to the classical methods in term of area of applicat ions. For computationa l electromagnetics, several numerical schemes have been introduc ed such as finite difference method, variational methods, differentia l variational schemes, finite element method and so on PAGE 25 25 [17]. Among those numerical schemes, the finite element method is known as a general method to solve differential equations. Maxwells Equations Maxwells equations are a set of four parti al differential equations that relate the electric and magnetic fields to their sour ces, charge density and current density. The individual equations are known as Gausss law, Gausss law for magnetism, Faradays law of induction and Amperes law. Figure 21 shows concept ual drawings for Maxwells equations Figure 21. Maxwells equations. A) Gausss law, B) Gausss law for magnetism, C) Faradays law of induction, and D) Amperes circuital law. Figure 21 A shows Gausss law that relates electrical charge within a close surface to the surrounding electric field. The differential form of the Gausss law is expressed as D (21) PAGE 26 26 where, D is the electric flux density and is the electrical charge density. The constitutive equation is written as DE (22) where, E is the electric field, and is the permittivity for diel ectric material. Figure 21 B shows Gausss law for magnetism which means the total magnetic flux through a closed surface is zero. As the dipole magnetic change only exists in real world, the divergence of magnetic fields cancels each other out. The differential form for Gausss law for magnetism is stated as 0 B (23) where, B is the magnetic flux density. Figure 21 C shows Faradays law of induction. The law states that a changing magnetic field produces an induced electric field, which is the operating principle for many elec tric generators. The differential form for Faradays law of induction is stated as t B E (24) Amperes law with Maxwells correction is show n in Figure 21 D. The law indicates that two factors can generate magnetic field; electrical current and changing electric flux density. The differential form for Amperes la w with Maxwells correction is written as t D HJ (25) where, H is the magnetic field. The constitutive equation between H and B can be stated as BH (26) where, is the permeability. PAGE 27 27 For electric field problems, dielectric materials are characterized by the permittivity However, conductive materials are characterized by the conductivity so that a different constitutive equation is used which is stated as JE (27) After taking the divergence of Ampere s law, the equation can be written as 0 tt D D HJJ (28) Equation 28 can be restated as t J (29) The equation is called the electrical continui ty equation. Practically all electromagnetic devices guarantee that the input current is equal to the out put current for the devices. Otherwise, electrical charge accumulates in the device or is produced by the device. Therefore, the electrical continui ty equation is normally written as 0 J (210) Overview of Magnetic Actuator Among many electromagnetic devices, magnetic actuators are the focus of this study. Magnetic actuators have widely been used as components of electrohydraulic valves, fuel injectors in engines of autom obiles, biomedical prosthesis devices for artificial organs, head positioners for comput er disk drives, loudspeakers and relays [18]. The magnetic actuator is an energy conversi on device or a transducer. This transducer transforms magnetic energy to mechanical ener gy. In order to use these actuators as precise components, inputs and outputs need to be controllable. For the input, magnetic circuit is used for electrical signal to be able to control the intensit y or direction of the magnetic field. The electrical signal is characterized as directed current (DC) and PAGE 28 28 alternating current (AC). For the output, me chanical system is used to render magnetic force to be used as controlled mechanical ou tput. Figure 22 shows the block diagram of a magnetic actuator. The electrical input can be direct current or alternating current. The mechanical output can be rotary motion or linear motion. The flex ibility of the input and the output enables broad applic ation of magnetic actuator s. In order to design a proper magnetic actuator for a specific appl ication, sometimes the analysis of the magnetic actuator is difficult because of complexities within th ree blocks: magnetic circuit, force factor, and mechanical system show n in Figure 22. As these blocks often involve the complex geometry, the analysis of magnetic actuator becomes complicated so that numerical schemes such as finite element method are required. Figure 22. Block diagram of a magnetic actuator In this study, only direct current (DC) is used as the input. The common types of magnetic actuators using DC are solenoid actuators, coil actuators, proportional actuators and rotary actuators. Solenoid actuators have an armature (movi ng part) and a stator (stationary part). The armature is made of steel laminates so that eddy current effect is reduced. Eddy currents create heat in a device that leads to energy loss. The stator has a solenoid coil which is wound into shapes like cylinder, cubi c, or parallelepiped. Solenoid actuators can produce linear motion of the armature whic h can be designed in a variety of shapes. PAGE 29 29 Rotary actuators generate rotary motion as the mechanical output instead of the linear motion of solenoid actuators. Rotary actuators have a rotor (moving part) and a stator (stationary part). The mo st common rotary magnetic act uator is a step motor. The step motor is used to work with a microproce ssor or a digital sign al processor which produces digital pulses for controlling the moti on. The processors send digital signals to the step motor in order to generate increment al step motion. One example of the step motors is a switched reluctance motor with t he rotor composed of only steel laminations. The advantages of the switched reluctanc e motor are high speed operation and the lowest construction cost because of the si mple structure and t he absence of permanent magnets [19]. In addition, mi nimum switching devices are required because only unidirectional current is needed. Because of above merits, the switched reluctance motor has several applications includ ing electric propulsion, fan and pump. Coil actuators are linear motion actuator s that can produce reversible force. Although forces on ferromagnetic material are used for the solenoid actuators, Lorentz force on currentcarrying coil is used for coil actuators. The Lorentz force is expressed as FNIBl (211) where, N is the number of turns for the coil, I is the amount of curr ent per one turn coil, l is the average turn length and B is the magnetic flux density which is perpendicular to the coil direction. As the force is proportional to the current, the direction of the force can be controlled by the direction of the curr ent. Usually, the magnetic field is provided by permanent magnets so that high flux den sity can be obtained without any power loss and temperature increments on the device. Coil actuators are widely used in PAGE 30 30 loudspeakers, where they are called voice co il actuators. Another popular application is a computer disk drive head actuation. Reluctance Method In order to analyze a m agnetic actuator, reluctanc e method can be used. The reluctance method is also called magnetic circ uit method, which is a simplified method to solve for magnetic fluxes and magnetic fields. For a magnetic actuator model with simple geometry, the reluctance method is an analytical method to estimate actuators characteristics such as approximate magnet ic fluxes, magnetic fields and magnetic force. In this study, magnetic fluxes and magnetic force by the reluctance method is used to compare with IBFEMs results. The reluctance method is based on the Amper es law. Amperes law in integral form can be expressed as a summation form as follows k kkk kk kB NIHllNI Hdl (212) where, NI is the ampereturns, kl is the line segment along the field intensity kH and k is the permeability of path segm ent k. The magnetic flux in the surface integral form can be redefined in discrete form as follows kkk B SBdS (213) where, and k is the magnetic flux and the magnetic flux on the k th path and k B is the magnetic flux density normal to the crosssection surface area kS. Substituting into the Amperes law in the summation fo rm, Equation 212 is stated as kk k kkl NI S (214) PAGE 31 31 As the divergence of flux density is zero, the fluxes through all segments have a same quantity So Equation 214 is rewritten as k k R NI (215) where, k k kkl R S called the reluctance. When the total reluctance R and the ampereturns, which is called magnetomotive forc e (MMF), are known, the flux can be expressed as NI R (216) The equation is similar to the familiar Ohms law of electric circuits: V I R where I is the electric current, V is the voltage potential, and R is the resistance. The reluctance method, which is also known as the magnetic circuit method, is graphically shown in Figure 23. Figure 23. Reluctance method (Magnetic circuit method) Even though the reluctance method is an easy way to estimate magnetic fluxes, the method has a limitation. The reluctance met hod ignores the fringing flux which means the expansion of flux in air as shown in Figure 24. PAGE 32 32 Figure 24. Fringing flux. A) Fluxes without fr inging flux, and B) Flux es with fringing flux When magnetic fluxes flow from the steel to air, the magnet ic flux tends to increase because of a large difference in permeability. The fringing flux makes airgap reluctance decrease so that the total fl ux increases. The upper and the lower steel cores are cylindrical with the diameter D and the air gap length is g as shown in Figure 24. The fringing flux can be ignored when t he ratio g/D is less than about 0.04, which means the reluctance method can be accurate However, the fri nging flux should be considered for larger g/D ratios [18]. Mo reover, most devices including complicated geometries or nonlinear materials such as permanent magnet ar e difficult to be analyzed using this method. PAGE 33 33 CHAPTER 3 IMPLICIT BOUNDARY FINI TE ELEMENT METHOD Overview of Finite Element Method Partial differential equations are used to express physical phenomenon such as heat transfer, electrostatics, magnetostatics, solid mec hanics, and so on. Exact or analytical solutions of those equations can be obtained under restri cted conditions of simple geometry and loading conditions. Other wise it is difficult to obtain analytical solutions. Since complex geometry and loading conditions are common for real physical or engineering situations, alternative methods are necessary to solve those problems. A popular method is the finite element method (FEM), which is a numerical method for solving partial differential equations. FEM has been widely used for solving engineering and mathematical problems both in academia and industry. FEM is a well established numerical method, but it is highly depended on the quality of the domain discretization using a conforming mesh. Oft en mesh generators are not reliable for 3D domain and if the mesh contains a few highly distorted el ements it can lead to i naccurate results. In order to avoid the need for generating conforming mesh, several meshless techniques has been developed [1], where nodes scattered within the domain are used without connecting them to form elements. Several appr oximation schemes have been developed that do not need elem ents or connection between nodes. Researchers have developed different types of meshless te chniques such as moving least square approximation method by Nayroles et al. and Belytschko et al., reproducing kernel particle method (RKPM) by Liu et al., Hp clouds method by Duarte et al. and so on. The meshless approximation schemes do not have Kroneckers delta property that is the value of the approximation at a node is not equal to the nodal value. Therefore, it is PAGE 34 34 difficult to impose essential boundar y conditions. Lagrange multiplier approaches, modified variational principles, penalty me thods and so on are used as remedies. The Lagrange multiplier method is the more ac curate method; how ever, this method increases the number of variables, causes the element matrix to be non positive definite and nonbanded which increases the computat ional time for solving the equations. Another approach to avoid traditional mesh generation is to use structured grid methods. Structured grids have regular, undist orted elements and are therefore easy to generate automatically. A structured grid does not approximate t he geometry of the analysis domain adequately therefore t he geometry has to be independently represented using equations. It is possibl e to use a variety of interpolation or approximation methods with st ructured grid methods. Belytschko et al. developed extended finite element method (XFEM) [8], wh ich uses implicit equations to define the geometry of the domain. And they us ed Lagrange multipliers to apply essential boundary conditions. Kantorovich and Krylov pr oposed a solution structure constructed using implicit equations of the boundaries, a solution stru cture that ensures that essential boundary conditions are satisfied at these boundaries. Rvachev et al used Rfunctions to construct implicit equations for the solution st ructure, and analyzed numerous partial differential equ ations [20]. The Rfunction method is a way to define implicit equations for domains using Bool ean operations between simple primitives. Shapiro and Tsukanov extended the Rfunction method to so lve nonstationary physical problems with timevarying geometries, and us ed transfinite Lagrangian interpolation to apply essential boundary conditions [21]. Hol lig et al. proposed the Webmethod which uses weighted extended Bspline basis to guar antee higher order continuous solutions PAGE 35 35 [22]. Their method has been used with Rfuncti ons or distance functions as implicit equations in order to construct the solution structures. In addition, Apaydin et al have extended the Webmethod to solve onedimensional elec tromagnetic problems and twodimensional electromagnetic wave equations [23][24]. In this study, the implicit boundary finite element method (IBFEM) is used, which has been developed by Kumar et.al [13][16] [65]. It has been used in the past to analyze 2D and 3D elastostatics and steady state heat transfer. IBFEM uses solution structures constructed using approximate step func tions of the bounda ries. Approximate step functions have a unit value inside the dom ain of analysis and transition sharply to zero at the boundary. An advant age of using step functions to construct the solution structure is that all the in ternal elements have identical element matrix. IBFEM has been used with traditional FEM interpolations also known as Lagrange inte rpolation as well as with uniform Bspline approximations [15]. IBFEM for Elastostatics Kantorovich and Krylov first proposed using solution structures constructed using implicit equations of boundaries to impose e ssential boundary conditions when using nonconforming mesh / structured grids. Since the nodes of the structured grid may not lie on boundary, traditional methods for applying boundary conditions cannot be used. For elastostatics, a solution structure is expressed as g asa uDuuuu (31) where u is the displacement vector g u is the grid variable represented by a piecewise interpolation, au is the boundary value functi on with prescribed values, s u is defined as the homogenous part of the solution, and (,..,)dindiagDD D is the diagonal matrix where PAGE 36 36 iD is an approximate step function that has been re ferred to as Dirichlet function [13] At any given point 3 R x, a Dirichlet function (or Dfunction) is defined as 0()0 () ()110() 1() x x xx xkD (32) where, k is an integer. The Dfunction has 1 kC continuity at () x where ()0 x is the implicit equation of the boundary at whic h the essential boundary conditions are applied. The Dfunction transiti ons between 0 and 1 in a narrow band 0() x. For IBFEM, very small values of is used (510 or smaller value) so that this band is very narrow and the Dfunction is a good approximation of a step function. Hollig et al. have used a weighting function that is similar to the Dfunction. Ho wever, they use relatively larger value of so that the weighting function is not a step function. The advantage of using an approximate step function for constructi ng solution structure is that all internal elements have identical element matrix. Ho wever, special techniques are needed to handle the large gradient of the step function within the band. There are several ways to construct the boundary value function. In this study, the boundary value function is defined as aa ii iNu u (33) where, iN is the shape functions which are identical to the basis functions used for the grid variable and a iu is the nodal values of the boundary elements in which the essential boundary conditions are applied. For example, if the prescribed value at a boundary is constant equal to 3, then the value, of a ll the support nodes of the boundary element, is PAGE 37 37 set equal to 3. The other nodes are set to zero. Figure 31 graphically depicts the solution structure. The element e1 is a boundary element and the others are internal elements. Figure 31. The solution structure with the essential boundary condition If the homogenous part of the solution and the boundary value function for the displacement vector are defined as s ssuv u and ,aaauv u. The strain can be stated as s a (34) For 2D problems in Voigt notation these strains can be defined as PAGE 38 38 ,TT ssssaaaa sauvuvuvuv xyyxxyyx (35) Using the constitutive equation, the stress is s asa C C (36) where, C is the elasticity matrix. The governing equation for linear elasticity is stated as 0in b (37) where, b is the boundary force, is the analysis domain. The essential boundary and natural boundary condit ions are defined as on,onut00uu nT (38) where, 0T is the traction vector, and u and t are the prescribed boundaries for the essential boundary and the natur al boundary conditions. Using the solution structure in the principle of virtual work for elastostatics we obtain tddddTsTTTa 0 uT ub (39) where, and u are small perturbation of t he strain and the displacement. Finite Element Formulation for Elastostatics Considering 2D elastostat ic problems, the boundary value function and the homogeneous part of the solu tion for the displacem ent field is stated as 1 1 1 100 00 00 00a N aaa a N s N s sg s NNN u NN v NN u DD NN v uXNX uXNX (310) where, iN is the shape function, and aX and g X are the predefi ned boundary value vector and the unknown nodal grid value vector respectively. Using these definitions, the strain can be expressed as PAGE 39 39 1 1 11 1 1 1 1100 00, 00 00a N a a aa N aa NN N s g N NNu N N x xx N N v yyy NN NN uv yxyx yx N N DD N xx N N DD yy NN NN DDDD yxyx XBX X 1 11 1200 00N g N NN gggDD N xx DD NN yy DDDD NNNN yxyx X BXBXBX (311) The virtual strain and stress are stated as ,ees BX CBX (312) where, eX and eX are the nodal grid values for the element e. Then, the weak form can be expressed in the following discrete form: 1 111e eetNE T ee e i NENENBE TTT T T eee eee iiid ddd a 0XBCBX XNbXB XNT (313) where, NE is the total number of elements in the domain, NBE is the number of elements on the boundaries t hat have natural boundary conditions specified. N contains the shape functions. Bspline element Using Lagrange elements, the tr aditional FEM only warrants C0 continuity of a solution between elements. When Bspline approximation is used, higher order continuity of the solution can be guaranteed. Figure 32 shows one dimensional solution using Lagrange interpolation and one using Bspline approximat ion scheme. Both PAGE 40 40 solutions can provide higher order continuity of the solution within elements. However, only C0 continuity of a solution between t he element e1 and the element e2 can be guaranteed using Lagrange interpolation as graphically shown in Figure 32 A. Figure 32. One dimensi onal field solution using two elements. A) Lagrange interpolation, and B) Bsp line approximation scheme. Bspline presentation has been used in com puter aided geometric design in order to represent curves and surfaces with higher order of continuity. As it is a parametric representation, the basis functions of Bsp line could be adapted to finite elements. Burla and Kumar [15] have solved elas tostatic problems using Bspline element with IBFEM. They reported that solutions of IBFEM hav e faster convergence rate using Bspline element. In IBFEM, uniform Bsplines can be us ed, which means nodes have been equally spaced in the parametric s pace. The Bspline shape func tions for each element are constructed to have independent parameter s pace such that the parametric range is from 1 to 1. Namely, one di mensional Bspline element with C1 continuity has three support nodes and three basis functions. The bas is functions are quadratic functions expressed as follows PAGE 41 41 0000102 1101112 2 22021221 Naaa Naaar Naaar (314) where, iN are the basis functions, r is the parametric value varying within 1,1 and ija are coefficients to be solved. Using these basis functions, the approximated functions of two adjacent elements can be stated as follows 1 1 01210122 23(),()ii TT ii ii iiuu frNNNufrNNNu uu (315) where, i f and 1 i f are the approximated functions and iu are the support nodal values. As given Bspline elements warrants C0 and C1 continuity, two continuity requirements provides the following restricted conditions : 1(1)(1)iiff, and 1(1)(1)iiff rr (316) where, (1)i f is the approximated value at the end point in the i th element 1(1)if is the approximated value at the start point in the i+1 th element, (1)if r is the tangential value at the end point the i th element, and 1(1)if r is the tangential value at the start point in the i+ 1 th element. The two conditions lead to following eight independent equations. PAGE 42 42 0 0 0 1 10 21 21 2 2(1) (1) (1) (1) (1)(1) 0 and 0 (1)(1) (1)(1) (1) (1) N r N N N NN rr NN NN rr N N r (317) Additionally, the basis functions should form a partition of unity, 2 01i iN. This provides another independent equation. Using the nine independ ent equations, the nine coefficients can be solved. The basis function can be stated as follows 0 1 2 2111 848 1 31 0 44 111 848 N N r N r (318) The Bspline element using these basis func tions is called Quadratic Bspline element (QBS). Using the same approach, t he Bspline basis functions with C2 continuity can also be obtained. The Bspline element with the basis functi ons is defined as Cubic Bspline element (CBS). Two and three dimensional Bspline elem ents are created by taking product of one dimensional Bspline elements. For exam ple, 2D quadratic Bspline element has nine supports nodes and nine basis functions 2D cubic Bspline element has sixteen supports nodes and sixteen basis functions. Timoshenko Beam Example Several elastostatic problems have been solved using IBFEM [13][16]. Here we provide one example for elastostatics from [15]. A cantilever beam with relatively large PAGE 43 43 thickness is created. The lengt h of 1.0m and the thickness of 0.2 m are used. According to Timoshenko beam theory, the tip deflection is given as 231 (45) 38tipP tLL EI where P is the applied load at the tip, E is the Youngs modulus, I is the moment of the inertia of the beam cross section, is the Poissons ratio, L is the length of the beam, t is the thickness of the beam. Assuming that the material is isotropi c, Youngs modulus and Poisson ratio are defined as 210 E GPa and 0.25 The uniform shear load of 0.1 M Pa is applied at the tip. The expected tip deflection is 54.886910 m using the Timoshenko beam theory. A B Figure 33. Displacement surface plots in ycomponent. A) Using 4 node bilinear elements B) Using 9 node Bspline elements Figure 33 shows displacement contours us ing two element types. In the case of Figure 33 B, where quadrat ic Bspline elements are used, the computed tip displacement is much closer to the analytical solution. PAGE 44 44 CHAPTER 4 IMPLICIT BOUNDARY FINI TE ELEMENT METHOD FO R 2D MAGNETOSTATICS Solution Structure for 2D Magnetostatics Under magneto static or quasistatic assump tion, if material is homogenous and isotropic, and only current induces the m agnetic field, the governing equation is 1() AJ (41) where, A is the magnetic vector potential, is the permeability of a material and Jis the current density vector. The co nstitutive equation is defined as 11 HBA (42) where, Bis the magnetic flux density and H is the magnetic field density. Figure 41 shows 2D assumption where, J xy Jk, and A xy Ak. That is, the current only flows in the zdirection and the magnetic potential only exists on the direction normal to the plane of analysis. The magnetic flux density and the magnetic field are vector fields in the plane xy. Figure 41. 2D magnetostatic problem PAGE 45 45 If the material is homogeneous, Equation 41 is expressed as the gradient of the scalar function A as following 111 AA A J xxyy (43) The constitutive equation for the 2D problem can be rewritten as 11 A A yx HBij (44) The solution structure for the magnetic vector potential A is defined as ()()()()()()gasaADAAAA xxxxxx (45) where, g A is a grid variable that is defined by piecewise interpolation or using Bspline approximation over a structured grid. Aa is the boundary value function which has a value equal to the prescribed bou ndary conditions at the boundaries, () Dx is a weighting function defined such that ()0 D x at boundaries where essential boundary conditions are applied so that ()=()aAAxx at these boundarie s. The boundary value functiona A is constructed by interpolating nodal values within elements. The nodal values are selected such that at the boundar y it has a value equal to the specified boundary condition. Using the weighted residual me thod with the solution struct ure in Equation 43, the weak form for 2D magnet ostatics is obtained as 1 1()()()()()()sssssa t A AdAJdAHAAd (46) where, A s is the virtual magnetic potential vector and t H is the tangential component of the magnetic field. PAGE 46 46 Finite Element Formulation for 2D Magnetostatics The grid variable g A is interpolated within each element as T g gA NA where, TNis a row matrix containing the shape functions and g A is a column matrix containing the nodal values of the grid vari able. Similarly, the boundary value function is represented within each element as T aaA NA where, aA is column matrix containing the nodal values assigned such that a A has the prescribed value at the boundary. Note that the same shape f unctions are used to interpolate g A and a A If all the essential boundary conditi ons are homogeneous, so that 0 A is the only prescribed boundary condition, then the boundary value function a A is zero everywhere and can be eliminated from the solution structur e. Otherwise, nodes near the boundary are assigned values of a A equal to the prescribed value. For 2D problems, the gradients of the boundary value function a A is expressed as 12[]T aa aaa i i i jN AA AA xxx BA (47) The gradient of the homogenous part of the solution s A is stated as 12[]T ss s gg i ii i jjN AAD ADNA xxxx BA (48) Using Equations 47 and 48, the weak form is rewritten as 1 1 1 111e e eNE T T gg e i NENBENE TTT T TT ggga ete iiid J dHdd ABBA ANANABBA (49) where, NE is the total number of grid elements and NBE is the number of grid elements on the natural boundary. [] B is decomposed into two matrices 1[] B and 2[] B such that PAGE 47 47 only 2[] B contains derivates of t he Dfunction which can have ve ry large values near the boundary. They are expressed as 1 1[]i ij jN BD x B (410) 2 2[]iji j D BN x B (411) The element matrix to be assembled into the global equations can be obtained as 1 123eT eeee ed KBBKKK (412) 1 111eT e ed KBB (413) 1 222eT e ed KBB (414) 11 31221eTT e ed KBBBB (415) As 2 B has terms containing derivatives of the Dfunction, it is nonzero only within the narrow transition band near the boundary. Ther efore, for all in ternal elements and boundary elements without ess ential boundary conditions 2e K and 3e K are zero. For boundary elements 1e K is evaluated by subdividing thes e elements into triangles and integrating only within tri angles that are inside the geometry. For boundary elements with boundary conditions, the volu me integral for computing 2e K and 3e K can be converted to surface integr als because they contain 2 B which is nonzero only within the narrow transition band near the boundary. The components of 2e K, can be expressed using in dex notation as PAGE 48 48 11 2 111 2 2 1 11 2 111 2 1e eN e e NN N e NNNNNN DD d xy SYMNN NNNN DD dnd xy SYMNN K (416) where, dn is the increment on the normal di rection of the boundar y. Using index notation, the above equat ion is rewritten as 1 2ee ijijeKNNd (417) where, 2 01k kD d x (418) In the preceding equation, the volume integral in 2e K has been converted into a combination of surface integral al ong the boundary and an integration over Similarly, the third term, 3e K, can be stated using the index notation as 11 3ej e i ijjike k kkN N K NNd xx (419) where, 01k kD Dd x (420) All element stiff ness matrices of 1e K, 2e K, and 3e K are evaluated using Gaussian quadrature. For surface int egrals, the boundary within elem ent is approximately by sufficiently small straight line segments to achieve accuracy. PAGE 49 49 Solution Structure for 2D Magneto statics with Permanent Magnet When permanent magnets exist on analysis domain, a different constitutive equation is used for the domain with the permanent magnets. The permanent magnetic domain has the magn etization vector, 0M being a part of constitutive equation. The constitutive equation can be stated as 000 r BHM (421) where, 0M is the magnetization vector, 0 is the vacuum permeability called magnetic constant 7410H/m, r is the relative permeability which is the ratio of the permeability of a medium to 0 The magnetic field can be expressed as 000 011rr HBMBM (422) where, 01r is the reluctivity. If the value of the reluctivit y becomes constant, then the governing equation is restated as 00AMJ (423) The equation becomes a nonlinear Poisson s equation for the vector potential. Multiplying the weight function A on both sides and integrating, Equation 423 becomes 00() ddAMAJA (424) Using Greens first identity for vector fields, the left hand side term becomes 00 0000() ()()d dd AMA AMAAMA (425) Using the divergence theorem, the se cond term on the right hand side becomes PAGE 50 50 0000()() ddAMAAMAn. (426) Using the identities of FGGF and F GTFGT Equation 426 becomes 0000()()ddAMAnAAMn (427) The weak form becomes 0000()d ddd AA AMAAMnJA (428) When a homogenous Neumann boundary condition is applied on the boundary, the tangential component of H 00()t HAMn, is set equal to zero. Thus, the line integral can be ignored. The w eak form can be rewritten as 00dddAAAMJA (429) In case of the 2D magnetostatic problem, the magnetic vector pot ential and the current density have only z components and the magnetizat ion vector is in the xyplane as shown in Figure 42. Figure 42. 2D magnetostatic problem with pe rmanent magnet PAGE 51 51 The 2D weak form becomes 00AAAAdddMJ (430) If 0 x y M M Mi j and A= 00AA x yzyx A ijk i j the weak form with permanent magnet can be rewritten as 0xyAAMMAAA ddd yx J (431) Using the solution structure of Equation 45, the weak form with permanent magnet is obtained as 1 1 0xy(A)(A) MM(A)(A)(A)ss ss ssad AA dJdd yx (432) Solution Structure for Multiple Materials When multiple materials are involved in t he analysis, there can be discontinuity in the magnetic field at t he interface even though A is continuous. To allow discontinuity in the magnetic field, separate gr ids are used for each material as shown in Figure 43. Figure 43. Two grids for two materials PAGE 52 52 At the interface the elements from neighboring grids overlap. A solution structure is needed for these overlapping elements to ensure that A is continuous while flux density B and magnetic field H can be discontinuous. Burla et al [16] suggested a modified solution structure using two different structur ed grids and solved an elastostatic problem for composite microstructures. According to [16], the modifi ed solution structure at the material interface boundary is called interface solution structure as shown below 121()() g gg A DADA xx (433) where, g i A is the field interpolated or approx imated within the el ement from grid i, (1,2) i (()) D x is the approximat e step function and () x is the implicit equation of the interface curve (represented using signed distance function). When the solution structure in Equation 433 is used in the weak form, element matrix for elements on the interface boundary contain terms t hat involve the gradient of (()) D x. The gradient of (()) D x is very large near the interface. The te rms in the volume integral for computing the element matrix are separ ated into those that contai n the gradient of the shape functions and those that contain the gradient of the Dfunction. As explained earlier for elements on the external bounda ries, the gradient of the Dfunction is zero outside a narrow band near the interface boundaries. Ther efore, the volume integral for terms containing the gradient of t he Dfunction can be converted into surface integrals for efficient computation. These techniques are described in deta il in [16] for elastostatics and have been adopted here for magnetostatics. Moreover, these techniques have been extended to solve multimaterial problem s with more than two materials in this study. PAGE 53 53 At the interface between materials with different magnetic permeability, the required interface conditions, expressed in terms of the magnetic vector potential are 1122nAnA (434) 12nAnA (435) which means the tangential component of the magnetic field and the normal component of magnetic flux density are continuous. It is obvi ous that if the magnet ic vector potential is continuous, that is 12 AA, then the interface condition (Equation 435) is automatically satisfied. However, the normal component of magnetic field and the tangential component of flux density can be discont inuous. To allow this discontinuity in the magnetic field and the m agnetic flux density, separat e grids are used for each material. Using two grids, the solution structure in Equation 433 combines the interpolation within the overlapping elements at an interface to represent the solution near the interface. When the Dfunction in unity, the solution is given by 2 g A A, which is the solution from the second grid and when the Dfunction is zero, the solution is given 1 g A A, which is the solution from first grid. In the region where the Dfunction varies from zero to unity, the solution is a blend of t he solutions from the tw o grids. This way of constructing the solution structure ensures t he continuity of the solution throughout the analysis domain. It also allows the derivative (and magnetic field and fl ux density) to be discontinuous at the interface. This pr operty can be seen from the gradients of the solution structure as shown below 12 121ggg gg jjjjjAADAD DADA xxxxx (436) In this expression, the firs t term and the third term ar e continuous at the interface boundary while the second term and the fourth term are discontinuous at the interface PAGE 54 54 boundary because the gradient of the Dfunction is zero for 0 and nonzero for 0 Therefore, these terms provide independent slopes on the two sides of the interface allowing discontinuous flux densit y when necessary and at the same time producing a continuous flux density if 12 g g A A Figure 44 graphically shows plots of components for the interface solution structure at the material interface boundary. Figure 44 A and B are weight functions for the grid 1 and the grid 2. When the solution st ructure behaves like a sinusoidal function, Figure 37 C and D represents A field distributions from the grid 1 and the grid 2. A B C D Figure 44. Component plots of interface soluti on structure at the ma terial boundary. A) 1() D B) () D C) 11() g DA, and D) 2() g DA PAGE 55 55 Figure 45. Interface solution st ructure at the material boundary Figure 45 shows the interface solution st ructure at the mate rial boundary following by Figure 44. The xaxis an d the yaxis represent grid elements and the field value. Part 1 and Part 2 represents t he grid 1 and the grid 2. Elem ents of e1 and e2 belong to Part 1 and elements of e2 and e3 belong to Part2. The element e2 is from each grid overlapping at the interface and the solution structure is used to blend solution from the two grids. In order to apply the interfac e solution structure to a m odel containing more than three materials, preprocessing is required to create a contact list containing contact pairs. Each contact pair associated with tw o materials is predefined and then the interface solution structure can be appli ed at the interface boundary based on each contact pair. Three parts with own grids are shown in Figure 46. For t he first part (Figure 46 A), three boundary elements belong to the inte rface boundary elements. One boundary PAGE 56 56 element belongs to the outer boundary elements. For the thir d part (Figure 46 C), three boundary elements belong to the interface boundary elements. Five boundary elements are outer boundary elements. Remai ned elements are inner elements. A B C Figure 46. Three parts with each own grid. A) Part 1 with the grid 1, B) Part 2 with the grid 2, and C) Part 3 with the grid 3. Using the interface boundary elements of each part, three contact pairs are defined. One is the pair of the 1st part and the 2nd part, which is defined as contact 1. Another is the pair of the 1st part and the 3rd part, which is defined as Contact 2. The last one is the pair of the 2nd part and the 3rd part, which is defined as Contact 3. These pairs are shown in Figure 47. A B Figure 47. Contact pairs. A) Contac t 1 and Contact 2, and B) Contact 3 For the contact pairs, the interface soluti on structure can be applied. Figure 48 shows four elements for the first par t. For the outer boundary element e1, the solution structure PAGE 57 57 is applied. For the other elements from e2 to e4, the interface solution structure can be applied. Figure 48. Four el ements in Part 1 In order to express severa l interface boundaries among mult imaterials, the Dfunction is redefined asijD where the subscripts i and j indicates the ith and the j th grids or parts. So then the generalized interf ace solution structure is restated as ijij1 j i g g g A DADA (437) As the element e2 has the 2nd contact pair, the interface solution structure becomes 3 113131 g g g A DADA (438) where, 13D is the Dfunction with the interface boundary between the 1st grid and the 3rd grid. Similarly as the element e3 is related to the 1st contact pair, the interface solution structure is 1212121H g g g A DAA (439) As the element e4 is linked to the 1st and 2nd contact pairs, the interface solution structure can be stated as 3 1121313121211g g gg g A DADADADA (440) PAGE 58 58 Magnetic Force Computation Several techniques for computing magnetic forces can be found in literature [25][30]. These include the virtual work principl e [25][28], Maxwells stress tensor method, equivalent source method, the force density method to list a few. These approaches have been implemented using FEM and theref ore can be used with IBFEM. Assuming that the exact equation of the surface is av ailable (preferable as a parametric equation), it is easier to implement a method that int egrates surface force densities to compute the nodal force. The weak form for solid mechani cs problems is the principle of virtual work, which can be stated as follows: T mddd C f ut u (441) where, C is the stiffness matrix, is the strain matrix, t is the traction force on specified surface, and mf is the magnetic force density. The force density is defined as follows In a current carrying conductor, the force due to magnetic field is the Lorentz force, given by c fJB (442) where, B is the magnetic flux density. In a linear, isotropic and noncompressible ferromagnetic material, 21 2fH f (443) The virtual work done by the force densities can be evaluated as follows: 21 2mdHddf u uJB u (444) PAGE 59 59 The Lorentz force is evaluated easily as the body force term. The equivalent nodal forces due to this body force can be computed by integrating the virtual work over the volume of each element as follows eT bdFNJB (445) The volume integration of each element is evaluated using Gaussian quadrature. To compute the force on the ferromagnetic ma terials, the volume integral of the force density ff must be changed to surface integr al. For a ferromagnetic object with permeability, 1 surrounded by a medium whose permeability is 2 the permeability can be considered to change from one value to the other over a band along the boundary whose width, measured in the normal direction, is n The limit of 0n represents the discontinuous va riation at the boundary. The gradient of the permeability within this band can then be written as n n (446) where, n is the direction normal to the boundary between the two materials with different permeability as shown in Figure 49. Figure 49. Material interface boundary: 1 (Inner material permeability) and 2 (Outer material permeability) PAGE 60 60 Using Equation 446, the force on the fe rromagnetic material is expressed as 2 1011 () 22 1 () 2nddnd n dd T uHHHHnu HHnu (447) The surface force density term can be evaluated by expressing the square of magnetic field as a function of permeability. If ther e is no surface current then the tangential component of the magnetic field tH does not vary across the boundary and can be treated as a constant. Similarly, the normal component of the magnetic flux density nB is constant (by Gausss law) and does not vary across the boundary even though the permeability is different on the two sides of the boundary. The squar e of the magnetic field can be expressed as t he sum of the squares of th e tangential com ponent of the magnetic field tH and the normal component of the magnetic flux nBas 2 22 2 n t B HH HH (448) Both the tangential component of magnetic field and the normal component of the magnetic flux density can be treated as cons tants for the integration in computing surface traction since these quantities do not va ry in the direction normal to the interface between the two materials. Using the pr eceding equation for the square of magnetic field to compute the surface traction due to magnetic forces, Equation 447 is rewritten as 1 22 222 21 2 2111111 (()) 222n ttnB HddHBd nuun (449) As the material is assumed to be linear, t H can easily be changed to tB using the constitutive equation. The surface force density can be shown to be equal to PAGE 61 61 2222 12 112211 22tntn sBBBB fnn (450) where, i and in are the permeability and the out ward normal vector of the i th material. If 21 it follows that 0sf therefore the direction of the surface traction Sf is opposite to n, which means that it acts in the direct ion of decreasing permeability. According to [30] an expression similar to Equation 450 has been deduced for the surface force densities between two linear media from a more general expression fo r magnetic force. The nodal forces at the boundary elements of each grid can be computed by integrating over the piece of the boundary that passes through the element. In other words, for a boundary element whose material property is i the nodal forces are computed as: 221 2eT tn si iiBB d FNn (451) The boundary passing through each element is approximated by lines (2D) or triangles (3D) for the purpose of int egration and Gaussian quadr ature was used along each line segment or triangle. PAGE 62 62 CHAPTER 5 IMPLICIT BOUNDARY FINI TE ELEMENT METHOD FO R 3D MAGNETOSTATICS Governing Equation and Weak Form for 3D Magnetostatics Several alternate formulations have been proposed in literature for 3D magnetostatic analysis using finite element method [31][39]. A formulation based on magnetic vector potential, A is used in this study. The governing equations for 3D magnetostatics can be expressed in terms of magnetic vect or potential as in AJ (51) where, is the domain of analysis. T he boundary of the analysis domain consists of regions with specified nat ural boundary conditions and regions that are open boundaries, which are used to artificially trunc ate the analysis domain when in reality it extends to infinity. Often homogeneous essential boundary conditions are used on these open boundaries as an appr oximation if the boundary is far away from the sources. Several special techniques for modeling such open boundaries have been developed such as the infinite elements and asymptotic boundary condition [40]. Natural boundary conditions can be appli ed on boundaries (denoted as H ) with known tangential component of the magnetic field or on boundaries (denoted as B) with known normal component of the flux densit y. If these boundarie s are planes of symmetry then 0 nA on B and 0 nA on H To ensure uniqueness of the solution, the followin g essential boundary conditions are used to enforce these conditions [33]. B0 on nA (52) H0 on nA (53) PAGE 63 63 Figure 51 shows an example domain of analysi s which may contain regions of different materials (1 and 2) as shown. 12 is the interface surface between the two subdomains 1 and 2 as shown in Fig. 2. At the in terface, the tangent ial component of the magnetic field and the normal component of the flux density are continuous. Figure 51. Analysis domain and boundaries The weak form for these governing equatio ns and boundary conditions, obtained using the weighted residual method, is HBdddAAAHnJA (53) where, A is the vector weighting function. This weak form is used in the traditional FEM to compute the element matrices by int egrating the left hand side over the volume of each element. When a structured mesh is used for the analysis, the boundaries pass through the elements so that it is necessary to integrate over partial volume of the element that is inside the boundary. Severa l techniques [41][42] have been developed for integrating over partial elements approxim ated as polygons. Alternatively, tetrahedral PAGE 64 64 elements could be generated withi n these partial elements for integration purpose. Even though, this is like mesh generation within the boundary elem ents, the generated tetrahedrons are used only for Gaussian quadrat ure and not to represent the solution. The boundary is approximated duri ng these integration techniques by triangles but this approximation is independent of the structured mesh that is used for interpolating or approximating the field variables. So the geometry can be represented reasonable accurately even if a sparse mesh is used for the analysis. Solution Structure for 3D Magnetostatics For three dimensional magnetostatics, a solution structure for the magnetic vector potential ()Ax could be defined as ()()()()()() AxDxAxAxAxAxgasa (54) where, ()Dx is diagonal matrix whose co mponents are defined such that D()0 xii at the boundaries where essential boundar y conditions are applied on the ith component of A. Substituting the solution structure (Equati on 54) into the weak form (Equation 53), a modified weak form of the 3D m agnetostatic equation can be derived as ssssadddAAJAAA (55) where, s A is the virtual magnetic potential vector. The grid variable vector, g A, is interpolated within each element as T g g ANA where, TN is a matrix containing the shape functions and g A is a column matrix containing the nodal values of the grid variable vector. For brick elem ents with 8 nodes, the size of g A is 24 because the nodal degree of freedom is 3. Sim ilarly, the boundary value function, Aa, is defined by interpolating nodal values within each element as T aa ANA where, aA is column PAGE 65 65 matrix containing the nodal values of aA. These nodal values are assigned such that, at the boundaries, this function will have a valu es prescribed by the boundary condition. Using the solution structure, the magnetic flux for the boundary value function can be derived as ABAaCa (56) where, the curl matrix C B for the boundary value function is defined as 12... BBBBCCCC n, where, 0 0 0 Bii C ii i iiNN zy NN zx NN yx (57) for 1,iZN. The curl of s A can be computed as [] ABA s Cg (58) For convenience, C B is defined as a sum of two matr ices such that the first one only contains derivatives of the shape functi on and the second matrix contains the derivatives of the Dfunction. 12 BBBCCC, where 111121 221222... ... BBBB BBBBCCCC n CCCC n (59) 2233 32 11133 31 1122 210 0 0 Bii C ii i iiNN DD x x NN DD x x NN DD xx (510) PAGE 66 66 33 22 32 33 11 2 31 1122 210 0 0 Bii C iii iiD D NN xx D D NN x x DD NN xx (511) In the preceding equations, 1,2...in where, n is the number of nodes per element. The element matrix that is assembled in to the global equations can be defined as 123eT eCCeee ed KBBKKK (512) 111eT eCC ed KBB (513) 222eT eCC ed KBB (514) 31221eTT eCCCC ed KBBBB (515) Since 2 C B contains the derivates of D()xii, it is nonzero only within the narrow transition band near the boundary. Therefore, for all internal elements and boundary elements without essential boundary conditions 2 e K and 3 e K are zero. Within the transition band, the derivatives of D()x can have large magni tude. For the boundary elements with boundary conditions, the volume integral for computing 2 e K and 3 e K must be converted to surfac e integrals as follows to compute them accurately. 222 01 KBBeT eCC edd (516) 31221 01 KBBBBeTT eCCCC edd (517) To derive the preceding equations, we make use of the fact that 2 BC is zero except in the narrow band 0 Therefore, the volume integral is converted into a surface PAGE 67 67 integral along the boundary e and an integral over the tr ansition band (normal to the surface). Note that if is a signed distance function then 1 If the width of the band is very small, then one can assume that the shape functions are constant within the band, allowing the integral over to be determined analytically. Alternatively, the integration over can also be evaluated numerically. Current Density Computation In threedimensional space, magnetic force computation and multimaterial analysis has the same approach as in two dim ensional problems. Just as the magnetic vector potential becomes a vector for 3D problems, the current density also becomes a vector. The current density is a spatial func tion or a vector fiel d for 3D magnetostatic problems so that the curr ent density distribution shoul d be computed prior to 3D magnetostatic analysis. Assuming that 3D problem s are static, it is possible to compute electric and magnetic fields separately. In the first step, the currents in the conductor are computed. The computed values ar e introduced as sources for the 3D magnetostatics in the second step. Figur e 52 shows governing equations for electrostatics and magnetostatics and t he current density as sources. V is the electrical potential and is the electrical conductivity. The current density is defined as V J. As the governing equation for electrostatics is one of Laplaces equations, the solution structure and finite element im plementation for 3D electrosta tics is identical to ones in the previous 2D magnetostatics. Figure 52. Sequential analys is for 3D Magnetostatics. PAGE 68 68 CHAPTER 6 OPEN BOUNDARY TECHIQUES FOR IBFEM Overview of Open Bounda ry Techniques for IBFEM Researchers have developed several tools for the computation of the solution fields for open boundary problems, whic h are boundary value problems that need infinite solution domain. Sinc e electromagnetic field problems are the infinite domain problems, several tools, called open boundary techniques, are developed for the magnetostatic analysis. The techniques in cludes truncation of outer boundaries, ballooning, infinite elements, infinitesima l scaling, spatial transformations, boundary relaxation approach, database approach, hybrid approaches, asymptotic boundary conditions (ABC), and measured equation of invariance (MEI) and so on [40]. Each open boundary technique has both merits and dem erits as shown in Table 61 [40]. Unfortunately, no allpurpose powerful techni que exists, so a researcher should choose one of these according to their application. Table 61. Comparison of diffe rent open boundary techniques [40] Methods Advantages Disadvantages Truncation Simple Large number of unknowns Iteration Simple Remeshing Relaxation Simple More than 1 solutions Ballooning Accurate/sparse Only 2D Laplace equation Infinite elements Sparse Exterior inaccurate Infinitesimal scaling Accurate N onlinear equation and dense matrix Spatial transformation Spar se No current source Hybrid method Accurate Dense matrix 1storder ABC Sparse/efficient Circular/spherical boundaries MEI Sparse Selection of metrons Among open boundary techniques for FEM, only some techniques are suitable for use with IBFEM. Some techniques ar e not suitable for IBFEM because those techniques require the conformal meshing or have other limitations. The unsuitable techniques include the ball ooning method, the infinitesimal scaling approach, and the PAGE 69 69 database approach. One of the most simple and efficient methods is the ballooning method [43][46]. The method needs a special meshing technique as shown in Figure 61. An outer boundary of a discretized regi on called annulus is moved outward from a center recursively. All nodes lie on the radial lines exte nding from a center point o Using a fixed ratio r called the mapping ratio, new annulus can be added in the radial direction. Figure 61. Twodimensional interi or region with an exterior annulus Based on the mesh, the method can create fast and accurate results. However, there are disadvantages such as having to select the center point and the mapping ratio for wellestablished mesh. Additionally, the me thod is only applicable to 2D Laplaces equations. The special meshing technique is not suitable for IBFEM so that it can not adopt this ballooning method. Secondly, the infinitesimal scaling approach [40] also requires such an annulus between interior and exterior region making it unsuitable for IBFEM. Thirdly, the database appr oach, developed by Sun et al [47] and Chen et al [48], PAGE 70 70 requires preprocessing to assemble and elimi nate matrices of the exterior region. The information is saved in a database so that it is used during the solving procedure. In this research, IBFEM is extended to use the following techniques: the truncation approach, the asym ptotic boundary condition, and the infinite element. The approach involves the truncation of out er boundaries is to impose homogeneous essential boundary conditions at the far aw ay boundaries. The technique assumes that the distance of the outer boundary fr om the center is at least five times the size of the region of interest. The method is easy to be implemented, but it is computationally expensive due to the large analysis domain. Ther efore, there is no effort to adopt this method for IBFEM. Unlike the truncation appr oach, the other two techniques have implementation issues that are de scribed in the following sections. Asymptotic Boundary Conditions for IBFEM Considering computational efficiency, t he asymptotic boundary condition is one of the most attractive approaches because the method guarantees matrix sparsity [40]. It is very similar to absorbing boundary conditi ons in high frequency calculations. In static and quasistatic electromagnetic fields, t he asymptotic boundary condition is used. Many researchers had solved open boundary probl ems using this technique. Brauer et al. [49] introduced the magnetic vector potential based problems using vector asymptotic boundary conditions. As the a symptotic boundary conditions are derived using the spherical coordinate system, circ ular or spherical boundaries are needed for the vector asymptotic boundary conditions. In order avoi d the inconvenience of having circular boundaries, Chen, K onrad, and Baronijan [50] developed techniques to solve axisymmetric electrostatic problems wit h rectangular boundaries. Chen, Konard, and Biringer [51] introduced the vector asymptotic boundary condi tions for three dimentional PAGE 71 71 magnetostatic problems with box boundaries. Gratkowski et al. [52] solved the axisymmetric magnetostatic problems with arbitrary boundaries. Derivation of Asymptotic Boundary Conditions The magnetic vector potent ial is defined as follow BA (61) Using Equation (61), the Amperes law can be restated as 2211 BAAAAJ (62) where, 0 A that is the Coulomb gauge condition. When A=0 at infinity, the general solution of the Poissons equation is (') ()' 4 d Jr Arv rr (63) where, (')Jr is current density at position 'r and rr is distance from the current source. When r is much larger than r, the solution asymptotically becomes 00 2()(')'()' 44dd rr ArJrvJrrrv (64) where, r is the unit radial direction vector. When rr, (')'0d Jrv. Therefore, the magnetic vector potentia l approximately becomes 0 221 ()()'(,) 4d rr ArJrrrva (65) where, (,) a is a function in the spherical coor dinate system. The leads to following equation: 2 0rr A A (66) which is called the firstorder as ymptotic boundary condition (ABC). PAGE 72 72 Two Dimensional Asymptot ic Boundary Conditions The asymptotic boundary conditions were orig inally based on circular or spherical boundaries. In two dimensional domain, Gartko wski et al. developed this method for arbitrary outer boundaries [52]. Cartesian c oordinates and rectangul ar outer boundaries are used so that IBFEM can be easy to adopt this approach. For 2D Magnetostatic problems, the magnetic potential vector exists only on the z component. On the rectangular boundaries for IBFEM, the line segment 12 with the unit normal of (1,0,0) is shown in Figure 62. Figure 62. Structured gr id and rectangular boundary The derivative of A in the radial direction on the line segment becomes cossinzzzzzAAAAA xy rxryrxy (67) Substituting Equation 67 into Equation 66, the equation becomes 02 0zz zx AA y A xryrr (68) PAGE 73 73 Since 0cos x r and sin y r. Therefore, ABC on a line s egment 12 in Figure 62 can be written as 002zz zAA y A x xxy (69) which means the derivative of the magnetic ve ctor potential in the normal direction can be decomposed of two terms: the derivative of the magnetic vect or potential in the tangential direction and the magnetic vector potential. According to Gratkowski et al [52], t he coefficients of the boundary matrix can be calculated from 02ej i ijijijeN N KNNyNNd xyy (610) Similarly, the line segments with the unit no rmal of (0,1,0) can have the following ABC 002zz zAA x A y yyx (611) The coefficients of t he boundary matrix become 02ej i ijijijeN N KNNxNNd yxx (612) Infinite Elements for IBFEM Since Ungless and Bettess [53] first proposed the infinite element technique, many researchers have developed the technique with two approaches: decay interpolation function approach (or displacement descent formulation) and mapped finite element approach (or coordinate ascent formulation) Bettess [54][56] dev eloped the decay interpolation function approach for potential and elasticity problems. The method was extended to solve the Helmho ltz equations by Rahman and Davies [56], McDougall and PAGE 74 74 Webb [58] and Towers, McCowen, and Macnab [59]. This approach uses the finite element shape functions multiplied by a decay function such as the reciprocal decay function [53] or the exponential decay function [55]. The modified shape functions make a field variable decrease to the specified far field value. With in the element, the numerical volume integration is done us ing GuassLagueree quadra ture instead of GaussLegendre quadrature used in the finite element method because the integration range is different. The mapped finite element approach was in itially developed by Beer and Meek[60], and Zienkiewicz et al [61]. AbdelFattah et al [62] generalized the method for static analysis in one, twoand thr eedimensional infinite do main. The approach transforms a regular finite domain into an infinite dom ain using the mapping of the shape function and the numerical integrat ion of GaussLegendre quadratu re. Using this approach, open boundary electromagnetic field problems we re solved by Rahman and Davies [56], McDougall and Webb [58] and Gratko wski and Ziolkowski [63]. These two approaches for infinite el ement implementation have similar performance. For IBFEM, t he first approach was implement ed. When a structured mesh is used, the shapes of the infi nite elements are known so that the infinite elements can be treated as virtual elements. Based on the assumption, the decay function infinite element approach can be easily implemented wit hout creating extra structured elements on the infinite domain. Decay Interpolati on Function Approach The mapping functions used for infinite el ements are the same as shape function. They are denoted as i M 1 to in where n is the total number of nodes in the element. PAGE 75 75 When the element is an isoparametric elemen t, the shape functions for the infinite element are multiplied by the decay functions as follows (,)(,)(,)iiiNstfstMst (617) where, (,)iNst are the shape functions on the infinite domain and(,)i f st are decay functions. Figure 63 shows coordinate trans formation for the decay function infinite element. When the decay direction is the radial direction r in the global coordinate system, the infinite element can be transferr ed from the global coor dinate system to the parametric coordinate system as shown in Figure 63. In the parametric coordinate system, the decay direction becomes the positive s direction and 1,s When an element matrix, K is obtained taking a volume integr al of the element, GaussLaquerre integration is used instead of GaussLagendre scheme that is used for the traditional finite element method because t he integration is done from 1 to the infinity in the parametric space. Figure 63. Coordinate transformation fo r the decay function infinite element PAGE 76 76 The decay function should be unity at its own node,(,)iist, however do not have to be zero at the other nodes. Under these assump tions, the decay function can be one of several decay patterns such as reciproc al decay type, exponential decay type, sinusoidal decay type, and so on. The select ion of the decay function is contingent on the type of the problems being solved. For magnetostatics, an exponential decay pattern is suitable. Assuming that the decay direction is the positive sdirection, the exponential decay function is (,)exp[()/]ii f stssL (618) where, L is defined as decay parameter to det ermine the degree of decay. If the decay directions are the positive s and t directions, the decay function can be stated as (,)exp[()/]iii f stststL (619) Figure 64 shows 1D shape functions us ing exponential decay functions. The shape functions using exponent ial decay functions are 11 exp[(1)/] 2 s NsL and 21 exp[(1)/] 2 s NsL Figure 64 shows plots of t he shape functions for different decay parameter. If a decay pattern for the problem is know n, a value of the decay parameter can be estimated. Many unbounded potential problem s are governed by the reciprocal decay pattern. In this case, the decay behavior c an be expressed by using the exponential decay functions with proper decay paramet er. Supposing there are two functions: A r PAGE 77 77 and exp(/) B sL where A and B are arbitrary values, 1rr at 1s and 2rr at 1s two functions have same values at two points as follows 1exp(1/) A B L r (620) 2exp(1/) A B L r (621) After eliminating A and B, the dec ay parameter is estimated using 1r and 2r as follow 1 22lnr L r (622) A B C D Figure 64. 1D shape functions using exponenti al decay functions. A) L=1, B) L=2, C) L=0.5, and L=0.1 PAGE 78 78 Decay Function Infinite Element for IBFEM The decay function infinite element is easy to implement for a structured mesh. As the infinite elements are rectangular, the Jacobi an matrix of the infinite element is easy to compute. When the field variable is ze ro on far away boundaries, the grid elements for the infinite domain need not be created so t hat the infinite elements can be virtual elements. Figure 65 shows structured grids with the infinite elements. Figure 65. Structured grids with the infinite elements Element e1 belongs to stru ctured grid and the other el ements e2, e3, and e4 are the infinite elements. e2 and e4 are extended to infinity in one direction. e3 is extended in infinity in two directions. Figure 66 illustrates an example with thr ee elements: two finite elements and one infinite element. The boundary conditions are the solution 1 P at x=0 and 0P at infinity. PAGE 79 79 Figure 66. Domain of integration The differential equation is 2 20 dP dx (623) For the infinite element, the mapping func tion and the shape function are defined as 11 1, 2 s Ms (624) 1/ 11 1, 2sLs Nes (625) As 0P at infinity, 11K, coefficient of the infinite el ement matrix, is stated as follow 22/ 11 1111 11()sLNN dx Kdsepsds xxds (626) where, 2 1 1111 () M psM sL and 1 dx ds Substituting 1 2 L ss the coefficient becomes 111111 00()() 2ssdsL Kepsdsepsds ds (627) The integration is done using GuassLaguer ee integration wit h following modified weights and abscissas 2/2L newoldL WWe (628) 1 2 L ss (629) PAGE 80 80 A solution for this example is provided by Chari and Salon [63]. Figure 67 shows the reference solution by Chari and Salon [63] and the computed solution varying with the decay factor L It graphically represents that the solution is very sensitive to the decay factor. Thus, it is necessary to choose proper value of the decay factor according a physical phenomenon. Figure 67. Solution plots for 1D infinite element example PAGE 81 81 CHAPTER 7 RESULTS AND DISCUSSIONS Two Dimensional Magnetostatic Problems IBFEM was implemented by modifying a finite element program. Four magnetostatic examples are examined to valid ate this method and show the benefits of IBFEM. These examples are published in [66 ]. The first example is a two dimensional model of a coaxial cable. The second ex ample is a planar solenoid with a clapper armature [18]. Both problems were selected because they have analytical solutions for comparison. The third example is a planar solenoid with a damaged armature. This example illustrates that even when the geomet ry is changed the same grid still provides reasonable answers while with FEM a new me sh is needed that is harder to generate as the geometry gets complicated. The final example is a 2D model of a switched reluctance motor (SRM). For all these examples, the geometry was created in commercial CAD software (Pro/Engineer) and imported into the analysis software (IBFEM). Example 711: 2D Coaxial Cable A coaxial cable, which consists of an inner conductor, an insulator, and an outer conductor, is modeled. Due to circular symme try of the geometry, only a quadrant of the coaxial cable crosssection is created. Figur e 71 shows the coaxial cable model using three separate structured grids. The radii of the inner conductor, the insulator and the outer conductor are a, b and c. The inner and outer conducto rs carry the same amount of total current in opposite directions. T he total current flowing through each conductor is I The current flows in the axial direction (the zdirection) and the current density is assumed to be uniform. PAGE 82 82 Figure 71. 2D coaxial cabl e model with structured grids The analytical solution of the magnetic fiel d in circumferential direction can be derived as 1 2 1 1 1 22222,0 2, 2, 0, raIra rIarb H crcbrIbrc cr (71) where, r is the radial distance. The following va lues of current and radii were used in the numerical model: 1000 I A, 0.5a 1b and 1.5c mm. Figure 72. Magnitude of H field PAGE 83 83 Figure 72 shows the magnitude of the magnetic field t hat was computed using the quadratic Bspline elements. It shows that the maximum magnet ic field value is at the interface between the inner conductor and the insulator and has a value of 23.18210 A/mm. This is very close to the value obtai ned from the analytical solution which is 23.18310 A/mm. 102 103 106 105 104 log(The number of nodes) log(Magnetic flux error norm) IBFEM Q4 IBFEM B9 Figure 73. Convergence plot for H1 norm PAGE 84 84 102 103 106 105 log(The number of nodes) log(Error norm) IBFEM Q4 IBFEM B9 Figure 74. Convergence plot for L2 norm Figure 73 shows the convergence of H1 error norm which is the root mean square error in flux density field over the domain and can be defined as 1 2 1T ehehHd BBBB (72) where, eB is the exact value of the magnetic fl ux from the analytical solution and hB is the corresponding computed value. Fig. 74 shows convergence of L2 error norm of A(x,y) This error norm is defined as 1 2 2T ehehLAAAAd (73) where, e A is the exact value of the magnetic potential from the analytical solution and h A is the computed value. The plot shows fa ster convergence for Bspline elements but the rate of convergence decreases with increasing number of nodes. PAGE 85 85 Example 712: 2D Clapper Solenoid Actuator The clapper armature soleno id is composed of an armature, a stator, and coils. A twodimensional planar model is shown in Fi gure 75 where symmetry is used to model just half of the system. The dime nsions of the components provided in [18] were used in the model to compare with t he approximate solution. Rela tive permeability of the armature and the stator was assumed to be 2000r The stator winding has 200 turns and a current I = 2A. When current flows th rough the coil, an attractive force is generated between the armature and the stator, which produces linear motion similar to clapping. Figure 75 shows a ty pical grid that is used to model this problem where each part has its own grid and elem ents of these grids overlap only at the interfaces. Figure 75. 2D Clapper sol enoid with structured grid Brauer [18] has provided an approximate soluti on for this problem using the reluctance method and using FEM when the gap width is 2mm. The total force obtained by the reluctance method is 122.62N and the forc e computed by FEM is 135.9N. For the PAGE 86 86 IBFEM model, a grid consisting of quadratic Bspline elements was used for the results shown in Figure 76 B. E ssential boundary conditions are applied along the axis of symmetry to enforce symmetry. Figure 76. Magnetic vector potential and magnetic flux lines. A) Comsol, and B) IBFEM Figure 76 shows flux lines and the plot of magnetic vector potential when the width of the air gap is 2mm. Figure 76 A shows t he result using commercial FEA software (COMSOL). The result from IBFEM is shown in Figure 76 B. The patterns of flux lines are similar for both results and the values of magnetic vector potent ial obtained are very close. Using the magnetic fi eld computed using IBFEM, t he total force was 140.0N, which is quite close to the total forces obtained by reluctance method and FEM in [18]. Figure 77 shows the total force on the clapper armature as a func tion of the gap width. In the Figure, the computed values of tota l force for different gap width are compared with approximate (analytical) forces determi ned using reluctance method and Maxwells stress tensor. PAGE 87 87 0.5 1 1.5 2 2.5 3 3.5 4 x 103 0 200 400 600 800 1000 1200 1400 1600 1800 2000 Air gap [m]Total force per unit depth [N/m] Analytic solution Computed forces Figure 77. Magnetic force versus air gap length Example 713: 2D Clapper Solenoid Actuator with Artificial Damage The clapper armature in above example has been remodeled in this example with an artificial damage added by editing one of the edges of the 2D geometry as shown in Figure 78 A. As in the previous exampl e, air gap width equal to 2mm was used to compute the results shown in Figure 78 B. This example was created to show that significant increase in complexity of t he model geometry does not pose a problem and in fact the same grid that was used in the last example was used for this example too. The computed force on the armature wit h the damage was 83.56 N compared to 140N without the damage. Figure 78 B displays the contour and a plot of the magnetic vector PAGE 88 88 potential. The computed magnet ic potential varies smoothl y between the two materials (iron and air) in the damaged area despite the increase in geometric complexity. Figure 78. Clapper solenoid with artifi cial damage. A) Da maged region and boundary condition, and B) Flux lines and surface plot of A Example 714: 2D Switched Reluctance Motor A 2D planar model of the Switched Reluct ance Motor (SRM) is shown in Figure 79. SRM is a DC motor where the stator has win dings around the poles while the rotor does not have any windings. The motor is called 6/4 Pole SRM because the stator has 6 poles and the rotor has 4 poles. Figure 79. 2D Planar model of switched re luctance motor. A) Aligned position, and B) Unaligned position. PAGE 89 89 According to Arumugam et al [64], the main dimensions of the motor are: Stator core outer diameter = 16.51 cm. Stator bore diameter =9.3cm. Length of iron core =10.8cm. Stator pole arc =1.88cm. Height of stator pole =2.355cm. Length of air gap = 0.0255cm. Rotor pole arc =2.83cm. Height of rotor pole =1.95cm. Diameter of the shaft = 1.858 cm Number of turns per pole =222 Current is applied to the coils around the poles of the st ator sequentially to produce a torque on the rotor as it rotate s. The stator and the rotor are made of iron with relative permeability of 2000. A quarte r of the motor is model ed in both its aligned and unaligned orientation. The num ber of turns in the coil is assumed to be 500 and the current is 2 A. Currents only flows into co ils attached to the top pole of the stator. Essential boundary condition (A=0) is imposed on peripheries of the shaft and the stator and the vertical axis. Bspline elements ar e used for the analysis so that accurate results are obtained even with a sparse grid as shown in Figure 710 and Figure 711. The flux lines computed here are similar thos e computed by FEM [64]. Note that the air gap is very small compared to the average element size in the grid. Despite the geometric complexity and large number of parts and materials involved, this approach PAGE 90 90 for analysis yields results using structured gr ids comparable to results from traditional FEM using conforming mesh. Figure 710. Magnetic vector potential and magnetic flux lines in the aligned position PAGE 91 91 Figure 711. Magnetic vector potential and magnetic flux li nes in the unaligned position Three Dimensional Magnetostatic Problems Several 3D examples are examined here to validate th e implicit boundary method. The examples are published in [67]. The first example is an iron material in homogenous magnetic field. This example al lows us to study the continuity or discontinuity of components of the flux density computed using IBFEM at material interfaces. The second example is a 3D coax ial cable that has an analytical solution for comparison. The third and the fourth exampl es are two solenoid actuators with plunger PAGE 92 92 and clapper armatures. For each solenoid actuator, a cylindric al and a block like shape are analyzed. For the cylindric al actuator, the computed fo rce and magnetic flux density are compared with 2D FEM solutions and the analytical solution from [18]. Example 721: Iron Block in a Homogenous Magnetic Field The example of an iron cube in air subj ect to homogenous magnetic field has been used to verify a variety of formulations [4][5], [7], [51]. Figure 712 shows oneeighth of the system modeled considering it s symmetry. The relative permeability of iron cube is 1000. The modeled region is subjected to a homogenous magnetic flux density 0B in the zdirection. Figure 712. Iron cube in homogeneous magnetic field PAGE 93 93 The halflength of t he iron cube edge is a. The sy mmetry planes are x=0, y=0 and z=0. The planes; x=b, y=b, and z=b r epresent the outer boundaries. A homogeneous magnetic field is applied in the zdirection with the aid of boundary conditions. On the outer boundaries, the Dirichlet boundary conditions are: 02y B b A and 0zA on x=b and 02x B b A and 0zA on y=b. The dimensions used are 20amm 40bmm and the flux density magnitude is 01.0 T B. In addition to these, the following essential boundary conditions are also imposed B0 on (x=0 and y=0) nA (74) H0 on (z=0 and z=b) nA (75) As the magnetic flux density only exists in the zdirection, t he normal components of magnetic flux density must be zero on the symmetry planes and the tangential component of magnetic field must be zero on planes normal to the zdirection. Figure 713. Iron objects with t he same grid density. A) Cube B) Octagonal prism, and C) Cylinder In addition to modeling the iron cube, we have modeled two other shapes: octagonal prism and cylinder for the iron part as show n in Figure 713. Figure 713 shows the PAGE 94 94 cube, the right octagonal prism, and the cylinder with the same grid density. The height of three objects is 20mm. The edge length of the right octagonal prism is 20mm. The cylinder is 20mm in radius. Along the obs ervation line y=z=10mm, the location of boundary varies with change in the shape of t he iron parts as shown in Figure 714. Figure 714. Crosssections with the line y=z=10mm. A) C ube, B) Octagonal prism, and C) Cylinder For the cube, the discontinuity happens at x=20mm, for the ri ght octagonal prism, at x=24.14mm, for the cylinder at x=17.32mm. The total number of elements is 12167. Using the same number of elements and t he same boundary conditions as above, the variation on the three components of B al ong the line y=z=10mm are obtained as shown in Figure 715. Figure 715 A shows that magnetic flux densit y in the xdirection is continuous only when the shape of the iron is cube because only in this case the normal component is parallel to the xdirection. The discontinuity or continuity of magnetic flux density is successfully shown. But the discontinuity between two materials is not as sharp as in FEM. In order to improve IBFEM result, local grid refinement is needed. PAGE 95 95 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x (mm)Bx (T) Cube Octagonal Prism Cylinder A 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x (mm)By (T) Cube Octagonal Prism Cylinder B 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.5 1 1.5 2 2.5 3 x (mm)Bz (T) Cube Octagonal Prism Cylinder C Figure 715. Components of B along the line y=z=10mm. A) Bx, B) By and C) Bz PAGE 96 96 Example 722: 3D Coaxial Cable A coaxial cable consists of an inner conductor, an insulator, and an outer conductor. A coaxial cable problem is a wellknown problem for 2D magnetostatic problem. As the current density and the m agnetic vector potential has only one component in the zdirection, the analytical solution is easily obtained using the amperes law. Using this example, we tri ed to solve an electromagnetostatic problem sequentially. In the first step, the current density in each conductor is calculated through 3D electrostatics. In the next step, using the computed cu rrent density, magnetic field in circumferential direction is obtained thr ough 3D magnetostatic analysis. The calculated magnetic field is compared with the analytical solution. The governing equations for the electro magnetostatic problem are described as follows 0 in in V AJ (76) where, the current densit y for magnetostatics is V J. Due to circular symmetry of the geometry, only a quadrant of the coaxia l cable is created. Figure 716 shows the coaxial cable model using thr ee separate structured grids. Figure 716. 3D coaxial cable model with t he structured grid. A) Inner conductor, B) Insulator, and C) Outer conductor PAGE 97 97 The radii of the inner conductor, the insulator and the outer conductor are a, b and c. The height of the three components is d. The inner and outer conductors carry the same amount of total current in opposite di rections. The total current flowing through each conductor is I The current flows in the axial direction (the zdirection) and the current density is assumed to be uniform. The anal ytical solution of the magnetic field in circumferential direction can be Equation 71 from the 2D coaxial cable example. The following values of current and radii were used in the numerical model: 1000 I A, 0.5a 1 b 1.5c and 0.2dmm. The current density of the i nner conductor is computed to be 1273 A/mm2 and 254.6 A/mm2 in the outer conducto r. The electric conductivity of the conductors is set equal to 310/ Smm, and the insulator is equal to 310/ Smm. In order to obtain 1000 I A, the voltage difference in the top and the bottom surfaces is set to 0.25 V in the inner conductor and 0.05 in the outer conductor. Figure 717 shows the magnitude of the magnetic field that was computed using 8 node hexahedral elements. It shows that the maximum magn etic field value is at the interface between the inner conductor and the insulator and has a value of 23.35110/ A mm This is close to the value obtained from the analytical solution which is 23.18310/ A mm. Figure 718 shows the magnetic field in the hoop direction varying with the radius. H8 means hexahedron element with 8 nodes. H8S stands for H8 with smoothing. After smoothing, the result of IBFEM is very close to the analytical solution. Figure 719 shows the convergence pl ot for H1 norm using 8 node hexahedral elements. H1 norm is defined as Equation 72. The H1 norm decreases as the larger number of elements is used. PAGE 98 98 Figure 717. Magnitude of H field for 3D coaxial cable 0 0.5 1 1.5 0 50 100 150 200 250 300 350 Radius, r [mm]H in hoop direction [A/mm] Analytical Solution IBFEM H8 IBFEM H8S Figure 718. Magnetic field in the circumferential di rection versus. radius PAGE 99 99 103 105 Log(Num of nodes)Log(H1 Norm) Figure 719. Convergence plot for H1 norm Example 723: 3D Plunger Solenoid Actuator Solenoid actuators have armature (movi ng part) and stator (stationary part). The armature is made of steel laminates in order to reduce eddy current effect. The stator has solenoid coil which is wound into shapes such as cylinder, cubic, parallelepiped and so on. Solenoid actuators c an produce linear motion of t he armature and have several types of actuators depending on the shape of the armature. PAGE 100 100 Figure 720. Plunger solenoid act uator of axisymmetric geometry Figure 720 shows a solenoid actuator with plunger armature from [18]. The plunger is cylindrical and the solenoid is axisymmetric. The magnetic force only acts at the end of the plunger. The number of turn s N=400 and the current I=4A. The relative permeability of the stator, the arma ture and the stopper is 2000r and 1r in the coil. The dimensions are provided in t he Figure 720. Using the relu ctance method, the magnetic flux density in the air gap is B=0.1715 T and the magnetic force is F =14.7 N. Brauer [18] also provided the following 2D FEM re sults: B=0.170 T and F =19.34 N. Using the similar dimensions, two different plunger so lenoid actuators with different shapes were created as shown in Figure 721. PAGE 101 101 Figure 721. Top views of two plunger actuators. A) Cylindr ical plunger, and B) Brick plunger Figure 721 A is the top view of the cylinder actuator. Figure 721 B is the top view of the brick shaped actuator. The corners of the brick plunger are rounded. The former can be modeled as 2D axisymmetric magnetosta tic problem. However, the latter one can only be analyzed using 3D magnetostatics. Figure 722. 3D solid model of solenoid act uators with plunger armatures. A) Cylindrical plunger, and B) Brick plunger. PAGE 102 102 Figure 722 shows 3D solenoid actuator model created by Pro/Engineering. In order to reduce the computat ion effort, we created one fourth of the whole model and modeled air only over the armature. The armatu re is moveable in the ydirection. In order to obtain the same current density as in the 2D axisymmetric problem, a voltage difference is applied between the two faces of the coil part. The voltage difference is equal to 0.0546 V when the conduc tivity of the coil is 610/ Sm. The essential boundary conditions of 10.0273V and 20.0273V are imposed on each face of coil parts. The same voltage boundary conditions are applied for the brick actuator. Using the computed current density from 3D electrostatics, we can calculate magnetic flux density and magnetic field. The boundary conditions at the symmetric planes are applied as 0 nA. The computed magnetic flux density is shown in Figure 723. Figure 723. The magnetic flux density in the ydirection for two plunger solenoids. A) Cylindrical plunger, and B) Brick plunger According to [18], the analytical solution is 0.1715 T and the FE M result is about 0.170 T. The computed magnetic flux dens ity is between 0.114 and 1.977 T for the PAGE 103 103 cylindrical plunger. The comput ed magnetic field density in the ydirection is shown in Figure 724. Figure 724. The magnetic field in the ydi rection for two plunger solenoids. A) Cylindrical plunger, and B) Brick plunger Figure 724 A shows the computed magnetic fi eld in the ydirection is approximately equal to 51.36410 A/m which is computed using the re luctance method [18]. Using the magnetic field and flux density, the computed m agnetic force on the cylindrical armature is 18.33 N which is quite close to the fo rce calculated by 2D FEM. For the brick armature, the computed force is 19.56 N. A ccording to this analysis, the force on the brick armature is larger than that on t he cylindrical armature because the crosssectional area of the brick armature (321.59210 m) is larger than the area of the cylindrical armature (321.25710 m).Figure 725 shows magnetic force versus gap length. When the gap is changed from 2mm to 10mm, the magnetic force decreases. Even though the gap varies, the same mesh is used for all the models. For the cylindrical plunger solenoid, the number of nodes is 10 984. For the brick plunger solenoid, the total number of nodes is 10876. When the shape of the actuator is a PAGE 104 104 cylinder, the analytical solution using the reluctance method can be compared with the results of 3D cylindrical plunger model. The computed values are higher than the analytical values because the reluctance met hod ignores fringing effect. Magnetic force of the brick plunger is a little higher t han the force of the cylindrical plunger. 2 3 4 5 6 7 8 9 10 0 20 40 60 80 100 120 140 160 180 Gap length [mm]Magnetic foce [N]Plunger solenoid actuator Analytical solution Cylinder Brick Figure 725. Magnetic force versus gap length for plunger solenoids Example 724: 3D Clapper Solenoid Actuator One of other popular solenoid actuators is a clapper sol enoid actuator. Figure 726 shows a solenoid actuator with clapper armatu re from [18]. Clapper armature can have linear movement by using voltage or current control. Armature a nd stator are made of thin steel laminations with relative permeabi lity of 2000. This clapper solenoid actuator is axisymmetric. The number of turns N= 2000 and the current I=1A. The dimensions are provided in the Figure 726. Us ing the reluctance method, t he magnetic flux density in PAGE 105 105 the right air gap is B=0.56 T and the magnetic fo rce is F =240 N. According to [18], 2D FEM results are B=0.82 T on the inner pol e and 0.45 T on the outer pole and F =279.41 N. Using the similar dimensions, two differ ent shapes of the clapper solenoid actuators are created in shown in Figure 727. Figure 726. Clapper solenoid act uator of axisymmetric geometry Figure 727. Top views of two clapper actuators. A) Cylindr ical clapper, and B) Brick clapper. PAGE 106 106 Figure 727 A is the top view of the cylinder actuator. Figure 727 B is the top view of the brick actuator. The co rners of the brick armature are rounded. The former can be approximated for 2D axisymmetric magentos tatic analysis. However, the latter one should be analyzed using 3D magnetostatics. Figure 728. 3D solid model of solenoid act uator with clapper armatu res. A) Cylindrical clapper, and B) Brick clapper Figure 728 shows 3D solenoid actuator model created by Pro/Engineering. In order to reduce the computat ion effort, we created one fourth of the whole model and the air model that only covers the armature. T he armature is moveable in the ydirection. In order to obtain the same current density as in 2D axi symmetric problem, the voltage difference, equal to 0.4188 V, was applied between two fa ces of the coil part. The conductivity of the coil is 610/ Sm. The essential boundary conditions of 10.2094V and 20.2094V are imposed on each face of coil parts. The same voltage boundary conditions are applied for the brick actuator. PAGE 107 107 After computing current density distribut ion on conductor, we can calculate magnetic flux density and magnetic field. The boundary conditions at the symmetric planes are applied as 0 nA. The computed magnetic flux density is shown in Figure 729. Figure 729. The magnetic flux density in t he ydirection for two clapper actuators. A) Cylindrical clapper, an d B) Brick clapper The computed magnetic flux density is approxim ately 1.3 T for the cylindrical clapper. The computed magnetic field density in t he ydirection is shown in Figure 730. Figure 730. The magnetic field in the ydi rection for two clapper actuators. A) Cylindrical clapper, an d B) Brick clapper PAGE 108 108 Figure 730 A shows that the computed magnetic field in the ydirection is approximately 57.610 A/m Using the magnetic field and flux density, the computed magnetic force on the cylindrical armature with the forced area of 322.82710 m is 284.76 N which is quite close to the force calc ulated by 2D FEM. For the brick armature with the forced area of 323.51410 m, the computed force is 320 N. 1.5 2 2.5 150 200 250 300 350 400 450 500 550 Gap length [mm]Magnetic foce [N]Clapper solenoid actuator Analytical solution Cylinder Brick Figure 731. Magnetic force versus gap length for clapper solenoids Figure 731 shows magnetic force varyi ng with the gap length. The gap length is changed from 1.5mm to 2.5mm. The analytical solution fo r the cylindrical clapper actuator is lower than the computed one. The magnetic force of the brick clapper actuator is higher than that of the cylindric al clapper actuator because the forced area on the armature is larger. PAGE 109 109 Magnetostatic Problems with Permanent Magnets Two examples including the permanent magnet are examined here to validate the implicit boundary method. The first example is borrowed from the m odel library of the commercial software (Comsol). Using a Ushaped permanent magnet, the computed magnetic fluxes are compared to the result s of Comsol. The second example is a 3D cylindrical permanent magnet that has an analytical solution for comparison. Example 731: UShaped Permanent Magnet A Ushaped permanent magnet is analyz ed using IBFEM and compared with a commercial software (Comsol). The drawing of the model is shown in Figure 732 A. The example is a 2D magnetostatic problem t hat includes three different regions; iron, air, magnets. R1 and R2 represent the magnets or magnetization regions. R3 represents iron. R4 indicates air. R1 and R2 are magnetized in the xdirection. R1 has magnetization equal to 7.5e5 A/m, and R2 has magnetization equal to 7.5e5 A/m. Homogenous essential boundary condition s are applied on the boundaries of R4 Figure 732. Ushaped permanent magnet m odel. A) Geometry model, and B) Conformal mesh PAGE 110 110 The finite element model uses 10286 quadrilate ral elements to solve this problem as shown in Figure 732 B. The density of conf ormal mesh is denser near the permanent magnet. Figure 733 and Figure 734 show surfac e plots and contour plots for magnetic vector potential. Figure 733 is the result from Comsol and Figure 734 is from IBFEM. In IBFEM, four node bilinear elements are us ed and the total number of the elements is 6600. According to the figures, both result s are similar not only graphically but also numerically in terms of the magnetic vector potential. Figure 733. Surface and contour plot s for magnetic potential from Comsol PAGE 111 111 Figure 734. Surface and contour plot s for magnetic potential from IBFEM In order compare two results precisely, t he magnetic field in the xdirection and the magnetic field density norm are obtained along a line. The line for comparison is shown in Figure 735. Figure 736 shows the magnet ic field in the xdirection. IBFEM 3220e stands for IBFEM result using 3220 elements. The line plots show that the two IBFEM results are quite close to Comsol result. Figure 735. The observation line for the comparison PAGE 112 112 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 5 4 3 2 1 0 1 2 3 4 x 105 x axisHx Comsol IBFEM 3220e IBFEM 6600e Figure 736. Magnetic field in the xdirection on the line 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 105 x axisH norm Comsol IBFEM 3220e IBFEM 6600e Figure 737. Magnetic field norm on the line PAGE 113 113 Example 732: Three Dimens ional Cylindrical Magnet A 3D cylindrical permanent magn et in free space is modeled as shown in Figure 738. Considering the symmetry, one fourth of the permanent magnet is created in order to reduce computation. Figure 738. 3D solid model for cylindrical magnet in the air Boundary conditions are 0 nA on the symmetry planes (x=0 and z=0). The radius r and the height h of the magnet are defined as 1 m. The magnetization vector (0 M ) is 1 A/m in the ydirection. U nder the conditions, the analytic al solution of the magnetic field along the y axis is given by [2]: 02yM BB H AA (76) where, 221 2 Rz A dd 1 2 z B d and 0 if /2 2 if PAGE 114 114 Figure 739. Magnetic fiel d in the ydirection Figure 739 shows the magnetic field in t he ydirection. Eight node hexahedral elements are used and the total number of the elements is 20000. T he result graphically shows the discontinuity at the material interfac e. Figure 740 shows the line plots for the magnetic field along the y axis. The resu lt of Comsol was computed using a 2D axisymmetric model. As shown in the Figure, IBFEM result is close to the analytical solution and the comsol result. Some discrepan cy of the magnetic field is founded at the location where it is far from the interface boundary. The re sult may be more accurate when the number of elements increases, the quadratic or cubic el ements are used or an open boundary technique is applied. PAGE 115 115 0 0.2 0.4 0.6 0.8 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 y axisHy [A/m] Analytical Comsol IBFEM Figure 740. Magnetic field along y axis Magnetostatic Problems with Open Boundary Techniques Several examples with the open boundary techniques are examined here to validate the implicit boundary method. The first example is Example 731 with smaller analysis domain. The second is the entire model of the previous sol enoid actuator from Example 712. The third example is to compute Lorentz force between two wires. These examples show how the open boundary techniques could improve results. Example 741: UShaped Permanent Ma gnet with Open Boundary Techniques For the previous example (Example 731), the large unbounded domain (6 x 4) was necessary with the truncation approach. The small air domain (1.6 x 1.4) is created as shown in Figure 741. PAGE 116 116 Figure 741. Ushaped permanent magnet with the reduced air domain Two open boundary techniques (asympt otic boundary conditions and decay function infinite method) are applied on the outer boundaries. Figure 742 shows contour plots for magnetic ve ctor potential for each technique. The total number of four node bilinear elements used in the model is 8888. The decay factor1 L for the decay function is used in infinite element met hod. Figure 742 A show s different contour pattern when the homogenous essential bo undary is applied on the outer boundaries. Figure 743 shows magnetic fiel d in the xdirection on the observation line as shown in Figure 635. EBC stands for the homogenou s essential boundary conditions. The results using ABC are quite similar to the Comsol result. The result obtained using EBC has the largest error among them. PAGE 117 117 Figure 742. Contour plots fo r magnetic vector potential. A) Homogenous essential boundary conditions, B) Asymptotic boundary conditions and C) Decay function infinite element with L =1. 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 5 4 3 2 1 0 1 2 3 4 x 105 x axisHx Comsol EBC ABC Infinite Figure 743. Magnetic fi eld in the xdirectio n on the observation line PAGE 118 118 Example 742: 2D Clapper Solenoid Ac tuator with Open Boundary Techniques The second example for open boundary techniqu es is extended from the previous solenoid actuator (Example 712). The whol e actuator including smaller air domain is created as shown in Figure 744. Three open boundary techni ques are applied on the other boundaries. The magnetomotive force NI and material properties of components are the same as in Example 712. Figure 744. 2D Clapper solenoid actuator When the gap length is 2mm, Figure 745 show s the contour plots for magnetic vector potential. Among the three resu lts, the result using the homogenous essential boundary conditions differs from the other two. The quadratic bspline elemen ts are used for the analysis. The total number of elements is 7049. The computed forces for all open boundary techniques are as follows: Homogenous essential boundary conditions (EBC): 269 N/m Asymmetric boundary conditions (ABC): 276 N/m Decay function infinite element me thod with L=1 (Infinite): 282 N/m. The computed force using EBC is the lowest value. These computed forces can be also compared with the following va lues: the approximate solu tion using the reluctance PAGE 119 119 method is 245.2N/m and the force computed by FEM with the ballooning method is 272 N/m from [18]. The computed force using ABC is quite close to the reference force from [18]. Therefore, using ABC can produce the most reliable answer in three techniques. A B C Figure 745. Contour plots of magnetic vector potentia l. A) Homogenous essential boundary conditions, B) Asymptotic boundary conditions and C) Decay function infinite element with L =1. PAGE 120 120 1.5 2 2.5 3 3.5 4 x 103 50 100 150 200 250 300 350 400 450 500 Gap length [m]Total force per meter depth [N/m] Approximate solution EBC ABC Infinite Figure 746. Magnetic fo rce versus gap length Figure 746 show the total force on the clapper armature as a function of the gap width. All computed values of total forc e follow the approximate solution using the reluctance method. The computed force us ing ABC is the best answer among them, a force that exits between the force using EB C and the force using the infinite element method. Example 743: Force Betw een Two Parallel Wires Two long straight parallel wire s are separated by a distance d The wires carry the same amount of current in the same di rection as shown in Figure 747. Then an attractive force is crea ted between the two wires. PAGE 121 121 Figure 747. Two long parallel wires Using the Amperes law, t he magnetic field at the 2nd wire caused by the current 1 I can be obtained as 01 12 I B d According to the Lorentz force equation, the magnetic force per unit length at the 2nd wire is stated as 012 22 I I F ld (77) where, l is the wire length. The same amount of the magnetic force is obtained at the 1st wire. When 1 md and 121 A II the magnetic force per unit length becomes 7210 N/m. The radius of the wire is 0.05 m, and the surrounding air area of 21.91.9 m is modeled. Open boundary condition is appl ied on the outer boundary using the truncation method, the asymptotic boundary c onditions, or infinite element method. Quadratic Bspline el ements were used. The total number of nodes was 10520. The contour plots of the computed magnetic vect or potential are shown in Figure 748. The magnetic force per unit length is computed as follows; 72.5310 N/m using the PAGE 122 122 truncation approach, 71.6310 N/m using the asymptotic boundary condition, and 72.4210 N/m using the decay function infinite element method with L =0.10. Using the asymptotic boundary condition, the computed forc e is closer to the analytical solution when 1 md. The other methods show similar resu lts in terms of the contour plot and the magnetic force. It is because the decay length is so short that the decay function infinite element method is as sa me as the truncation approach. Figure 748. Contour plots fo r magnetic vector potential. A) Homogenous essential boundary conditions, B) asymptotic boundary conditions and C) Decay function infinite element with L =0.1. Figure 749 shows the magnetic force per uni t length versus the distance between the two wires. The distance d was varied from 0.2 m to 1.6m Within a range from 0.2 m to 0.8 m, the decay factor L =0.15 was used for the decay func tion infinite element method and using this value the computed force using this method were close to the analytical solution. The computed forces using all open boundary techniques are similar when the distance is equal to 0.2 m, 0.4m, or 0.6m; howev er, the force comp utation using the asymptotic boundary condition is much closer to the analytic solution as the distance increases, that is when the source approac hes the boundary. As the magnetic source (currant carrying wire) is closer to the out boundary, the error in the magnetic force increases when the truncation approach or the infinite element method is used. PAGE 123 123 Therefore, the force comput ation using ABC is clearly better than the other two open boundary techniques. Figure 749. The magnetic force vers us the distance between two wires Coupled MagnetoElastostatic Problems Two coupled magnetoelastostatic problems are examined to validate the implicit boundary method. The first example is a 2D example containing a clapper solenoid actuator and a cantilever beam. The cantilever beam is attached on the top surface of the clapper armature. The magne tic force from the actuat or deforms the cantilever beam. The second example is a 3D plunger sol enoid actuator workin g with one of three circular shaped plates. Even though the geomet ries of the plates are different and complicated, the density of t he structured mesh is the same for all three plates. PAGE 124 124 Example 751: 2D Clapper Solenoi d Actuator with Cantilever Beam A cantilever beam is attached on the top su rface of the armature. The cantilever beam is a beam fixed at one end. At the other end, this c antilever beam is attached to the armature. The attachment location varies as shown in Figure 750. The distances from the center axis to the attached location are 10mm, 20mm, and 28mm. The thickness of the beam is 5mm and the length of the beam is 80mm. The beam is made of aluminum of which material properties are 69 Gpa E 0.29 and 1r Figure 750. Planar clapper solenoid actuat or with a cantilever beam. A) d=10mm B) d=20mm and C) d=28mm This coupled analysis can involve nonlinear properties, nonlinearity that can occur from the large deformation of the beam, changing force according to the gap length, and the contact between the armature and the stator In this research, the nonlinearity is ignored during the analysis. The computed m agnetic force is used as loads for the cantilever beam. Four node bi linear elements are used for the analysis. The total number of elements is 1550. Fi gure 751 shows tip displacement versus ampereturns. As NI increases, the displacement on the tip increases. When d=20mm or 28mm, the displacements on the tip is much larger than when d=10mm; howev er, there is very small difference in tip deflection between the two cases: d=20mm and d=28mm. PAGE 125 125 100 200 300 400 500 600 700 800 900 1000 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 NI [Ampereturns]Displacement on the tip [mm] d=10mm d=20mm d=28mm Figure 751. Displacement at the tip versus NI with varying attachment location Next, the distance between the symmetric ax is and the tip of the cantilever beam is fixed at d=20mm, and the thickness of the bar is varied as shown in Figure 752. The thickness value, equal to 3m m, 5mm or 7mm, was used. Figure 752. Planar clapper solenoid actuator with a cantilever beam. A) t=3mm, B) t=5mm and C) t=7mm PAGE 126 126 Figure 753 shows the tip displacement versus NI When t=3mm, the larger displacement can be observed. When NI increases, the difference by the thickness also increases. The tip deflection is very sensit ive to the thickness of the cantilever beam. 100 200 300 400 500 600 700 800 900 1000 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 NI [Ampereturns]Displacement on the tip [mm] t=3mm t=5mm t=7mm Figure 753. Displacement on the tip versus NI with varying beam thickness Example 752: 3D Plunger Sole noid Actuator with Structures The 3D model of a plunger solenoid act uator shown in Figure 728 A works with one of three plates, one that is attached to t he top surface of the plunger armature. The top views of each structure and dimensions are included in the Figure 754. The first model is a solid plate, the second is a pl ate with one hole, and the third one is a plate with two holes. The thicknesses of all the plates equal to 2mm. All of the three structures are attached to the top surface of the plunger arma ture as shown in Figure 755. The structures ar e made of aluminum. PAGE 127 127 Figure 754. Top views of structures attached to the plunger armature. A) Solid plate, B) plate with one hole, and C) plate with two holes When NI =800 and the fixed boundary conditions are applied on the su rfaces of the structure at R=70mm, the m agnetic force from the actu ator results in downward deflection of the plunger armature. The struct ures are deformed by the movement of the armature as shown in Figure 756. The maximum displacements for the structures are 77.2410 mm in the solid plate, 61.08610mm in the plate with one hole, and 78.53910mm in the plate including two holes. Figure 757 shows the maximum displacement versus NI (ampereturns) for each plate. Among the structures, t he plate with one hole has the largest deformation. On the other hand, the solid plate has the highest stiffness among them. PAGE 128 128 A B C Figure 755. 3D plunger solenoi d actuators with A) Solid pl ate, B) plate with one hole, and C) plate with two holes PAGE 129 129 A B C Figure 756. Deformation due to magnetic force. A) Solid plate, B) plate with one hole, and C) plate with two holes PAGE 130 130 1000 1500 2000 2500 3000 2 4 6 8 10 12 14 16 x 106 NI (Ampereturns)Max. displacement (mm) Solid One hole Two holes Figure 757. Maximum displacement versus NI PAGE 131 131 CHAPTER 8 DESIGN AND ANALYSIS OF MAGN ETIC ACTUATORS USING IBFEM Design Criteria In this chapter, IBFEM is used as an act uator design tool. So lenoid actuators with clapper armature or plunger armature and coil actuators are designed under given design criteria. Among them, a lighter and more efficient design are needed for flapping wing micro air vehicle. The micro air vehicles are significantly small size aerial vehicle with the maximum dimension of about 15 cm. Fo r a micro air vehicle actuator, there are three given design criteria: coil winding area, iron weight and size. The given coil winding area wS is assumed to be 212 mm. The entire coil winding area cannot be assumed to be only the copper because the coil includes copper and other materials such as insulator. Thus, the relation between the coil winding area and the copper conductor is defined as the packing factor pF [18]. When the coil is wound tightly, the packing factor can be up to 75%. Thus, the packing factor is assumed to be 70%. The area of the coil wire is 325.0110cSmm when the number of the American Wire Gauge (AWG) is 40. Therefor e, the number of coil turns N is given as following equation 1680pw cFS N S (81) When the coil is wound to a cylindrical bobb in with the diameter of 4 mm, the total length of the coil becomes about 21.1 m. As the resistance per meter for the given AWG is 3.44 /m, the total resistance of the coil is approximately 72 In this research, it is assumed that the coil wire can allow current of 100mA to flow. The allowable current of the wire is based on plastic insulation. In case of the current source, LT3092, PAGE 132 132 manufactured by Linear technology, the ma ximum output current is 200mA. The input voltage range is from 1.2 V to 40 V so that LT3092 can be operated by using a small size battery. Thus, the current of 100mA can be obtained using LT3092. If N =1680 and the amount of current I can be controlled by a digital processor, NI (ampereturns) can vary from 0 to 168. NI is also called magnetomotive forc e. In order to compare several designed actuators in the later se ction, the magnetomotive force, NI =30, is used. When the maximum current (100mA) flows in the coils, the coils can have maximum energy dissipation as 20.72 [] WIR W. Another design criterion is the iron weight of the act uator. The iron occupy large portion of the weight of the ac tuator; however, iron has high re lative permeability so that usage of the iron material can intensify a m agnetic force of the ac tuator. Therefore, there is a tradeoff between the iron weight of the actuator and the usage of the iron. As the second design criterion, we assume that t he total iron weight of the iron portion can be up to 10g. In order to compute the iron weight of an actuator design, the mass density of iron is used as 3337874/7.87410/ kgmgmm. The third design criterion is the size of the actuator. In or der to install an actuator insi de the micro air vehicle, the actuator should fit inside of a cube with edge length of 2cm. In the beginning of a solenoid design, t he reluctance method is used so that a simple design is created under the design criter ia. Based on the initial design, geometry of the design is changed in order to determine a better design. Solenoid Actuators with Clapper Armature A solenoid actuator with a clapper armature is shown in Figure 81. The clapper armature can move in the ydirection when the coils carry current. The armature and the PAGE 133 133 stator are made of thin steel laminations with relative pe rmeability of 2000, a thin steel lamination that can reduce heating caused by eddy current. When the current flows through a coil, magnetic flux is created as shown in Figure 81. The magnetic flux follows a closed path. If magnit ude of the magnetic flux is known, a magnetic force can be estimated using the Maxw ell stress tensor method. Figure 81. Solenoid actuat or with clapper armature Supposing that the dimensions are 130.5llmm 253llmm 466.5llmm, 1 gmm and 1231wwwmm, then reluctances in the armature, the stator, and the gap can be obtained as follows 3 0123 1.59210A/Wb 1armature rlll R w (82) PAGE 134 134 3 0456 6.36410A/Wb 1stator rlll R w (83) 6 02 1.59110A/Wb 1gapg R w (84) Using the reluctance method, the magnetic flux can be stated as 76.28510Wbiarmaturestatorgap iNINI NI RRRR (85) As the magnetic flux is constant on the given path, the magnetic flux density and the magnetic field density can be obtained in each airgap as follows 46.28510T 1gapBNI w (86) 01 500A/mgapgapHBNI (87) Using Maxwells stress tensor method, t he normal magnetic pressure can be obtained as 2 2 4 0416.28310N 2gapFHwNI (88) When 30NI, the magnetic force per unit lengt h becomes 0.565 N/m. The magnetic flux and the magnetic field density are 21.88610T and 41.510A/m. Using the same geometry of the actuator Figure 82 A shows contour plot of magnetic vector potential computed using IBFEM. Figure 82 B shows magnetic field density in the ydirection. The computed magnetic field densit y in the gap is approximately equal to the analytical value of 41.510A/m. For this analysis, the asymptotic boundary conditions are applied on the outer boundaries. The com puted force is 0.633 N/m. The computed force and the analytical solution are quite close. The difference between them results from a fringing flux, which the reluctanc e method ignores. As the iron domain is 240mm, the iron weight per unit depth is 0.315/gmm. PAGE 135 135 Figure 82. IBFEM results. A) Contour pl ot of magnetic vect or potential, and B) Magnetic field density Based on the initial design in shown in Figure 81, three different designs are created by changing the geometry of the armature and the st ator. For all the designs, the rectangular area of the 2D model including the air domain, is 21313mm. Figure 83 shows three designs with dimensions. The dimensions shown are all in mm Even though the geometry is changed for each actuator the mesh densities of the structured mesh are same for all designs. Nine node Bsp line elements were used for the analysis. The total number of nodes was 2586. The asymptotic boundary conditions were applied on the outer boundaries. T he magnetomotive force,30NI was applied. Figure 83. Clapper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3 PAGE 136 136 Figure 84 shows contour plots of magnetic vector potentials for all designs. Based on those computations, Table 81 shows com parison for three designs in terms of the magnetic force and the iron weight. Figure 84. Contour plots of magnetic vector potential for cl apper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3 PAGE 137 137 Table 81. Comparison for th ree clapper solenoid actuators ( NI =30) Design 1 Design 2 Design 3 Force per m 0.691 N/m 0.509 N/m 0.681 N/mm Iron weight per mm 0.315/gmm 0.283/gmm 0.320/gmm According to Table 81, the first actuator design produces the largest force among them. The second actuator is the lightest one. Solenoid Actuators with Plunger Armature A solenoid actuator with a plunger armature is shown in Figure 85. The shape of the plunger armature can be a brick, cyli nder, or conics. There are two gaps of g and gs that the flux can enter and leave the plunger armature. Useful magnetic force is produced only at g Figure 85. Solenoid actuat or with plunger armature PAGE 138 138 The armature and the stator are made of steel with the re lative permeability of 2000. When the dimensions are 17lmm, 22.5lmm 30.9lmm 44.5lmm 51.5lmm, 63.5lmm, 0.1gsmm, 1gmm and 12341wwwwmm all reluctances can be calculated as follows 3 034 2.14810A/Wb 1armature rll R w (89) 3 01256 5.76910A/Wb 1stator rllll R w (810) 4 07.95810A/Wb 1gsgs R w (811) 5 07.95810A/Wb 1gg R w (812) Using the reluctance method, the magnetic flux is 61.13210Wbiarmaturestatorggs iNINI NI RRRRR (813) Based on the magnetic flux, t he magnetic flux density and the magnetic field density in each airgap can be computed as follow 31.13210T 1gapBNI w (814) 01 900.8A/mgapgapHBNI (815) Using the magnetic field in the gap, the us eful magnetic force can be calculated as follow 2 2 3 0211.0210N 2gapFHwNI (816) when 30NI, the force becomes 0.918 N. The m agnetic flux and the magnetic field density are 23.410T and 42.70210A/m. Figure 86 shows the computed results: contour plot of magnetic vect or potential and magnetic field density in the ydirection. PAGE 139 139 For the analysis, the same mesh density wa s used in the clapper solenoid actuator model. Nine node Bspline elements were us ed. Asymptotic boundary conditions were applied on the outer b oundaries. The computed total force on the plunger armature is 0.9130 N/m. As the tota l area of the iron is 238.8 mm, the iron weight per unit depth is 0.306/ gmm. Figure 86. IBFEM results. A) Contour pl ot of magnetic vect or potential, and B) Magnetic field in the ydirection Modifying the initial design in shown in Figure 85, three different designs are created as shown in Fi gure 87. The dimensi ons shown are all in mm For all designs, the area of the 2D model is as same as one of the clapper solenoi d actuator. The same mesh density is applied for all the designs. Nine node Bspline elements were used for the analysis. The total number of nodes wa s 2548. The asymptotic boundary conditions were applied on the outer boundarie s. The magnetomotive force, 30NI, was used. PAGE 140 140 Figure 87. Plunger solenoid ac tuators. A) Design 1, B) Design 2, and C) Design 3 Figure 88 shows contour plots of magnetic vector potentials for all designs. Based on those results, Table 82 shows comparis on for three designs in terms of the magnetic force and the iron weight. According to Table 82, the third plunger solenoid actuator model can produce the largest forc e among them. The second actuator model is the lightest one. PAGE 141 141 Figure 88. Contour plot s of magnetic vector potential for plunger solenoid actuators. A) Design 1, B) Design 2, and C) Design 3. Table 82. Comparison for thr ee plunger solenoid actuators ( NI =30) Design 1 Design 2 Design 3 Force per m 1.154 N/m 0.262 N/m 2.801 N/mm Iron weight per mm 0.306/gmm 0.294/gmm 0.310/gmm PAGE 142 142 Solenoid Actuators with a Combined Plunger & Clapper Armature A combined plunger & clapper armature is an armature that includes a clapper & plunger to increase the force. Three designs were studied. The sizes of those actuators are similar to one of the clapper solenoid act uator in the previous section. Magnetic force and iron weight for each ac tuator are characterized. Figure 89 shows three solenoid act uators with a mixture armature. The dimensions shown are all in mm For all designs, the area of the 2D model is as same as one of the clapper solenoid actuator. For the analysis, the same mesh density is applied for all the designs. Nine node Bspline elements were used for the analysis. The total number of nodes was 2674. The asymptot ic boundary conditions were applied on the outer boundaries. T he magnetomotive force NI is equal to 30. Figure 89. Three combined plunger & clapper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3 After the analysis using these models, contour plots of magnetic vector potentials are show in Figure 810. Based on those resu lts, Table 83 shows comparison for three designs in terms of the magnetic force and t he iron weight. According to Table 83, the third one can produce the largest force among them. The second actuat or is the lightest one. PAGE 143 143 Figure 810. Contour plots of magnetic vector potential for combined plunger & clapper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3. Table 83. Comparison for three combi ned plunger & clapper so lenoid actuators ( NI =30) Design 1 Design 2 Design 3 Magnetic force per m 0.961 N/m 0.413 N/m 1.30 N/m Iron weight per mm 0.342/ gmm 0.339/ gmm 0.340/ gmm PAGE 144 144 Coil Actuators In this section, three coil actuators ar e studied under the given design criteria. The coil actuators include permanent magnets so that higher magnet ic flux can be produced without the increase of NI. Moreover, the usage of the permanent magnets can reduce heating on the actuator. The co il actuator has magnetic force in the conductive coil instead of one on the armature of the sol enoid actuator, a magnet ic force that is computed using the Lorentz force equation. Figure 811 shows three coil actuators. The dimensions shown are all in mm For all the designs, the rectangular area of t he 2D model including the air domain is 21818 mm The coil winding area is the same as in previous solenoid actuators. Two permanent magnets have different direction with the same remanent flux of 0.8 T. The first design is based on a typical voice coil actuator in a loudspeaker. The second is a design to remove the iron lami nates from the first design. The third design has a gap at the bottom portion of the iron laminate, whic h allows the coil to move freely in the downward direction. For the analysis, the same mesh density is applied for all designs. Four node bilinear elements were used. The total number of nodes was 12036. The asymptotic boundary conditions were applied on the outer boundaries. NI is equal to 30. Figure 811. Three coil act uators. A) Design 1, B) Design 2, and C) Design 3 PAGE 145 145 After the analysis using these models, the magnitudes of the magnetic flux densities are show in Figure 812. As the Lor entz force is proporti onal to the magnetic flux density, a stronger magnetic flux density near the moving coils can create a larger force. The first design has the strongest m agnetic flux near the moving coils. Based on those results, Table 83 shows comparison fo r three designs in terms of the Lorentz force on the moving coil and the iron weight of the actuator. According to Table 84, the first coil actuator can produce the largest force among them. The second design is the best design if the actuator weight is the most critical criterion. Figure 812. Magnitude of B fields for three coil actuators. A) Design 1, B) Design 2, and C) Design 3. PAGE 146 146 Table 84. Comparison for coil actuators ( NI =30) Design 1 Design 2 Design 3 Force per m 21.0 N/m 2.64 N/m 9.12 N/m Iron weight per mm 0.428/gmm 0 /gmm 0.394/gmm The Best Actuator among the Designed Actuators Using IBFEM, four types of actuators were examined su ch as the clapper solenoid actuator, the plunger soleno id actuator, the combined plunger & clapper solenoid actuator, and the coil actuator. Among them, the coil actuator is the best actuator considering the given design criteria. Among se veral designs for the coil actuators, the first design is used for flapping wings of the micro air vehicle as shown in Figure 813. The first design could generate the largest force under the given magnetomotive force so that the magnetic force can make the largest deformation of flapping wings among them. Figure 813. The best magnet ic actuator among several designed actuators Considering the 2D model, the best coil act uator is designed as 3D solid model as shown in Figure 814. Consi dering the symmetry, one fourth of the model is created using commercial software, Pro/engineering. Figure 815 shows the dimension of the PAGE 147 147 3D coil actuator and the directions of t he permanent magnets. T he dimensions shown are all in mm Figure 814. Sold model of the 3D coil actuator Figure 815. The top view of the 3D coil actuator PAGE 148 148 The moving coils carry current. The coil has a conductivity value of 610 S/m. In order to perform the analysis, eight node brick elements were used. The total number of the nodes was 9651. According to the design criterion, the ma gnetomotive force, NI can vary from 0 to 168. When NI is equal to 30, the magnetic flux density of the coil actuator is shown in Figure 816. Figure 816. The magnitude of the magnetic flux densit y of the coil actuator The magnetic flux density, B, vary from 1.17 to 2.33 T nearby the moving coils, so the computed force is 0.114 N. When NI is equal to 10, 30, 60, 90, 120 or 168, the computed Lorentz force is shown in Figure 817. Figure 817 shows the linear relation between the magnetomotive forc e and the Lorentz force. T he maximum force is 0.683 N. PAGE 149 149 0 20 40 60 80 100 120 140 160 180 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 NI [A.turns]Lorentz force [N] Figure 817. Lorentz force of the coil actuator versus NI Coupled MagnetoElastostatic Pr oblems with a Fla pping Wing Model The best coil actuator is embedded in a mi cro air vehicle with flapping wings as shown in Figure 818. The coil ac tuator can fit inside of the fuselage. The moving coil of the coil actuator is atta ched to the structure. Figure 818. The coil act uator with flapping wings PAGE 150 150 Four structures are designed as shown in Fi gure 819 structures that are created as surface models. The thickness of the st ructure assumed to be equal to 0.45mm. Additionally, it is assumed t hat the structure is made of aluminum. All the structures have a similar top view; however, they have different front and side views. The dimensions shown are all in mm. As the length in chordwis e direction is 26 mm, two coil actuators can be placed along the chordwise direction because the length of the coil actuator is 12.2mm according to Figure 815 Figure 820 also shows the four surface models. Figure 819. Four structures wit h flapping wings. A) Design 1, B) Design 2, C) Design 3, and D) Design 4. PAGE 151 151 Figure 820. Four surface structures with fl apping wings. A) Design 1, B) Design 2, C) Design 3, and D) Design 4. In order to perform analysis of thin shelllike structures the surface model analysis, IBFEM [65] have been extended to use 3D s hell elements that are 3D elements with three degrees of freedom per node. These she ll elements were used for the analysis. The total number of nodes in model is 666. Considering the symmetr y of the geometry, half of the structure was modeled for the analysis. The structured mesh and the boundary conditions are shown in Figure 821. The magnetic force acts downward so that the wing produces upstr oke. When one coil actuator is used, the applied force PAGE 152 152 varies from 0 to 0.342 N at the edge. Usi ng two coil actuators, the magnetic force can be double. Figure 821. Structured mesh and boundar y conditions of the first design As the wing produces upstroke, the wing can also create downstroke. Using the coil actuator, the same amount of Lorentz force can be produced in the opposite direction by changing the direction of the current. When NI is equal to 30, the wing upstrokes and downstrokes of the four designs are s hown in Figure 822. For the upstroke, the maximum displacement on the tip is 32.63310 mm in the first design, 34.23810 mm in the second design, 23.76510mm in the third design, and 11.96810mm in the fourth design. Magnitudes of the tip deflections during the wing upstroke are the same as ones during the wing downstroke because t he magnetic force is proportional to the displacement. Thus, the tip deflection can be double including both strokes. Among the four designs studied here, the last desi gn produces the largest tip deflection. PAGE 153 153 Figure 822. Displacement in the zdirecti on during the wing stroke. A) Design 1, B) Design 2, C) Design 3, and D) Design 4. PAGE 154 154 When NI varies from 1 to 168, the tip deflect ion as a function of the magnetomotive force is shown in Figure 823. The first thr ee designs are too stiff for our application. 0 20 40 60 80 100 120 140 160 180 0 0.2 0.4 0.6 0.8 1 1.2 1.4 NI [A.turns]Tip displacement [mm] Design 1 Design 2 Design 3 Design 4 Figure 823. Tip displacement versus NI for wing upstroke PAGE 155 155 CHAPTER 9 CONCLUSIONS AND FUTURE WORK Conclusions The main research contributions made in this thesis are to extend Implicit Boundary Finite Element Method to Perform magnetostatic analysis Perform coupled magnetoelastostatic analysis Model multimaterial system Implicit Boundary Finite Element Method (IBFEM) for analysis using structured mesh has been demonstrated for 2Dand 3Dmagnetostatic models and coupled magnetoelastostatic models. This approach directly uses the geometry imported from CAD systems without generating a conforming mesh. Structured meshes are easy to generate and the elements ar e regular and not distorted as in traditional finite element mesh. Furthermore, the inter nal elements are identical to each other and have the same stiffness matrix thus reducin g the computation required. In order to perform magnetos tatic analysis, IBFEM require s the capability for the multimaterial analysis. Ea ch material has its own structured mesh, and a modified solution structure is applied at the interfac e elements. The multimaterial analysis in magnetostatic problems allows magnetic force to be computed, a m agnetic force that includes the magnetic surface force dens ity and the Lorentz fo rce. The magnetic surface force density is a force on ferromat erial, and the Lorentz force is a force on a current carrying coil due to the magnetic flux. A fter computing these forces, the forces can then be used in a subsequent structural analysis to perform coupled magnetostaticelastostatic analysis. PAGE 156 156 IBFEM has been evaluated using several 2D and 3D magnetostatic examples. Although those examples were complicated in terms of geometry, IBFEM could solve the problems without any dist orted element. The results of those examples were compared to one using the reluctance method or to analytic solution when available. Based on the results of the magnetic flux density and the magneti c field density, the magnetic force was computed accurately Through 2D problems, IBFEM presented faster convergence using quadratic Bspline elements, the method that created very accurate solutions with relatively fewe r nodes. Moreover, IBFEM could solve large electromagnetic problems. In case of a 3D model with complex geom etry, the structured mesh could allow relatively fewer elem ents to be used for the analysis. Consequently, the number of equations is smaller than the traditional finite element method. Three open boundary techniques were im plemented in the IBFEM software and evaluated using several examples. The open boundary techniques studied here include the truncation method, the asymptotic boundary conditions, and the decay function infinite element method. Several examples were used for comparing these three methods including permanent magnet, sol enoid actuator, and two wire examples. Among the three methods studied, the asym ptotic boundary conditions produced the most robust results no matter where magnetic sources were located. The computed magnetic force could be s ubsequently used in the structural analysis. So, IBFEM could reveal relations hip that exists bet ween the magnetomotive force NI and the displacement of t he structure. This relation is an important factor in designing a magnetically actuated structure. Several examples in 2D or 3D were created using a variety of structur e geometry. Through the coupled magneto PAGE 157 157 elastostatic analysis, the relationships bet ween magnetomotive force and displacement were studied for these examples. IBFEM was used as a magnetic actuator desi gn tool. In order to design an efficient magnetic actuator for a flappi ngwing micro air vehicle, IBFEM was used to analyze several magnetic actuators created under the gi ven specifications. These specifications include the range of NI the size of a magnetic actuat or and the weight of a magnetic actuator. Several magnetic actuators were ex amined such as clapper solenoid actuator, plunger solenoid actuator, combined clap per & plunger solenoid actuator and coil actuator. Among them, the coil actuator created the largest magnetic force for the given specifications. This magnetic force was used in structural analysi s of flappingwings modeled using 3D shell elements. Through t he coupled magnetoelastostatic analysis, the relationship between t he magnetomotive force of th e coil actuator and the tip displacement of the flappingwing structure was determined. Future Work This research focused on the steady static analysis so that the coupled problem is modeled as the weakly coupled magnetoelastostatic analysis. This assumption allows us to perform the magnetostatic and elasto static analysis sequentia lly. If the magnetic field changes due to the elastic deforma tion, the problem becomes nonlinear and a strongly coupled nonlinear analys is is needed. The magnetic field may change because the structural deformation may cause perm anent magnets or circuits attached to the structure to also move re lative to each other. If the structure undergoes large deformation then geometric nonlinearities must be included. In additi on, if the coupled problem is a dynamic problem, then the electr ical circuits, the magnetic circuits and the PAGE 158 158 structural analysis must be solved as a coupl ed problem. Another source of nonlinearity is due to contact between the armature and the stator. The open boundary techniques should be further studied and extended for 3D magnetostatic problems. The open boundary technique s were studied in this thesis only for two dimensional problems. Asymptotic Boundary Condition (ABC) was found to be the most effective; however, th is approach requires the user to specify the origin of the magnetic source. The accuracy of the results depends on careful choice of this origin. Further, research is needed to explore more effective ways of applying open boundary condition. The magnetic force computation should be extended to nonlinear materials and permanent magnets. In this thesis, the m agnetic force computat ion was implemented and studied only for linear ferrromaterial and cu rrent carrying coils. In order to analyze magnetic actuators that involve moving permanent magnets, we need to implement force computation for permanent magnets. In order to do dynamic analysis for flapping wings operated by magnetic actuators, nonlinear characteristics in geometry, forc e and kinematics must be considered. The wing motion creates large deformation so that a moment can be a function of deformation. When an armature of magnetic actuator moves toward stator, the magnetic force can vary according to the location of the armature. In addition, the contact problem occurs between the armature and the stator. PAGE 159 159 LIST OF REFERENCES [167] [1] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, "Meshless methods: An overview and recent developments," Computer Methods in Applied Mechanics and Engineering vol. 139, pp. 347, 1996. [2] K. M. Lee, Q. Li, and H. Sun, "Effects of numerical formulation on magnetic field computation using meshless methods," IEEE Transactions on Magnetics vol. 42, pp. 21642171, 2006. [3] Q. Li and K. M. Lee, "An adaptive meshless method for magnetic field computation," IEEE Transactions on Magnetics vol. 42, pp. 19962003, 2006. [4] A. Bruyere, L. Illoul, F. Chines ta, and S. Clenet, "Comparison between NEM and FEM in 2D magnetostatics using an error estimator," IEEE Transactions on Magnetics vol. 44, pp. 13421345, 2008. [5] X. Liu, Y. M. Deng, Z. W. Zeng, L. Udpa, and S. S. Udpa, "ModelBased Inversion Technique Using ElementFree Galerkin Method and State Space Search," IEEE Transactions on Magnetics vol. 45, pp. 14861489, 2009. [6] O. Bottauscio, M. Chiampi, and A. Manz in, "Computation of higher order spatial derivatives in the multiscale expansi on of electromagneticfield problems," IEEE Transactions on Magnetics vol. 44, pp. 11941197, 2008. [7] K. Muramatsu, Y. Yokoyama, and N. Takahashi, "3D magnetic field analysis using nonconforming mesh with edge elements," IEEE Transactions on Magnetics vol. 38, pp. 433436, 2002. [8] T. Belytschko, C. Parimi, N. Moes N. Sukumar, and S. Usui, "Structured extended finite element me thods for solids defined by implicit surfaces," International Journal for Nume rical Methods in Engineering vol. 56, pp. 609635, 2003. [9] N. Moes, E. Bechet, and M. Tourbier "Imposing Dirichlet boundary conditions in the extended finite element method," International Journal for Numerical Methods in Engineering vol. 67, pp. 16411669, 2006. [10] Y. Abdelaziz and A. Hamouine, "A su rvey of the extended finite element," Computers & Structures vol. 86, pp. 11411151, 2008. [11] P. M. A. Areias and T. Belytschko, "A nalysis of threedimensional crack initiation and propagation using the extended finite element method," International Journal for Numerical Methods in Engineering vol. 63, pp. 760788, 2005. [12] T. Belytschko and T. Black, "Elastic crack growth in finite elements with minimal remeshing," International Journal for Nume rical Methods in Engineering vol. 45, pp. 601620, 1999. [13] A. V. Kumar, R. Bu rla, S. Padmanabhan, and L. X. Gu, "Finite element analysis using nonconforming mesh," Journal of Computing and Information Science in Engineering vol. 8, 2008. [14] A. V. Kumar, S. P admanabhan, and R. Burla, "Implicit boundary method for finite element analysis using nonconforming mesh or grid," International Journal for Numerical Methods in Engineering vol. 74, pp. 14211447, 2008. [15] R. K. Burla and A. V. Kumar, "Implicit boundary method for analysis using uniform Bspline basis and structured grid," International Journal for Numerical Methods in Engineering, vol. 76, pp. 19932028, 2008. PAGE 160 160 [16] R. K. Burla, A. V. Kumar, and B. V. Sankar, "Imp licit boundary method for determination of effective properties of composite microstructures," International Journal of Solids and Structures vol. 46, pp. 25142526, 2009. [17] S. R. H. Hoole, Computeraided analysis and design of electromagnetic devices New York: Elsevier, 1989. [18] J. R. Brauer, Magnetic Actuators and Sensors : John Wiley & Sons, 2006. [19] K. Vijayakumar, R. Karthikeyan, S. Paramasivam, R. Arumugam, and K. N. Srinivas, "Switched Reluctance Motor Modeling, Design, Simulation, and Analysis: A Comprehensive Review," IEEE Transactions on Magnetics vol. 44, pp. 46054617, 2008. [20] V. L. Rvachev and T. I. Sheiko, "RF unctions in Boundary Value Problems in Mechanics," Applied Mechanics Reviews vol. 48, pp. 151188, 1995. [21] V. Shapiro and I. Tsukanov, "Meshfr ee simulation of deforming domains," ComputerAided Design vol. 31, pp. 459471, 1999. [22] K. Hollig, C. Appric h, and A. Streit, "Introducti on to the Webmethod and its applications," Advances in Computational Mathematics vol. 23, pp. 215237, 2005. [23] G. Apaydin, S. Seker, and N. Ar i, "Weighted extend ed bsplines for onedimensional electrom agnetic problems," Applied Mathemati cs and Computation vol. 190, pp. 11251135, 2007. [24] G. Apaydin, S. Seker, and N. Ari, "Application of w ebspline method in electromagnetics," AeuInternational Journal of Electronics and Communications vol. 62, pp. 163173, 2008. [25] J. L. Coulomb and G. Meunier, "Finit eElement Implementati on of Virtual Work Principle for Magnetic or Electr ic Force and Torque Computation," IEEE Transactions on Magnetics vol. 20, pp. 18941896, 1984. [26] Z. Ren and A. Bossavit, "A new approach to eddy current problems in deformable conductors and some numer ical evidence about its validity," International Journal for Applie d Electromagnetics in Materials vol. 3, pp. 3946, 1992. [27] Z. Ren and A. Razek, "Local Force Computation in Defo rmableBodies Using Edge Elements," IEEE Transactions on Magnetics vol. 28, pp. 12121215, 1992. [28] Z. Ren, "Comparison of Different Forc e Calculation Methods in 3d FiniteElement Modeling," IEEE Transactions on Magnetics vol. 30, pp. 34713474, 1994. [29] J. Bastos and N. Sadowski, Electromagnetic modeling by finite element methods : CRC Press, 2003. [30] R. S. Grandia, V. A. Ga lindo, A. U. Galve, and R. V. Fos, "General formulation for magnetic forces in linear ma terials and permanent magnets," IEEE Transactions on Magnetics vol. 44, pp. 21342140, 2008. [31] J. Simkin and C. W. Trowbridge, "Use of the Total Scalar Potential in the NumericalSolution of Field Problems in Electromagnetics," International Journal for Numerical Methods in Engineering vol. 14, pp. 423440, 1979. [32] J. L. Coulomb, "FiniteElement 3 Dimensional MagneticField Computation," IEEE Transactions on Magnetics vol. 17, pp. 32413246, 1981. PAGE 161 161 [33] O. Biro and K. Preis, "O n the Use of the Magnetic Vect or Potential in the FiniteElement Analysis of 3Dim ensional Eddy Currents," IEEE Transactions on Magnetics vol. 25, pp. 31453159, 1989. [34] K. Preis, I. Bardi, O. Biro, C. Magele, W. Renhart, K. R. Richter, and G. Vrisk, "NumericalAnalysis of 3d Magnetostatic Fields," IEEE Transactions on Magnetics vol. 27, pp. 37983803, 1991. [35] R. Albanese and G. Rubinacci, "Magnetostati c Field Computations in Terms of 2Component Vector Potentials," International Journal for Numerical Methods in Engineering vol. 29, pp. 515532, 1990. [36] R. C. Mesquita and J. P. A. Bastos, "An Incomp lete Gauge Formulation for 3d Nodal FiniteElement Magnetostatics," IEEE Transactions on Magnetics vol. 28, pp. 10441047, 1992. [37] N. A. Demerdash and R. Wang, "Theoretical and Numerical Difficulties in 3D Vector Potential Methods in FiniteE lement Magnetostatic Computations," IEEE Transactions on Magnetics vol. 26, pp. 16561658, 1990. [38] A. S. Semenov, H. Kessler, A. Liskow sky, and H. Balke, "On a vector potential formulation for 3D electromechani cal finite element analysis," Communications in Numerical Methods in Engineering vol. 22, pp. 357375, 2006. [39] J. Jin, The Finite Element Method in Electromagnetics 2 ed: WileyIEEE Press, 2002. [40] Q. Chen and A. Konrad, "A review of finite element open boundary techniques for static and quasistatic elec tromagnetic field problems," IEEE Transactions on Magnetics vol. 33, pp. 663676, 1997. [41] G. Ventura, "On the elimination of quadrature s ubcells for discontinuous functions in the eXtended FiniteElement Method," International Journal for Numerical Methods in Engineering vol. 66, pp. 761795, 2006. [42] S. Natarajan, S. Bordas, and D. R. Mahapatra, "Numerical integration over arbitrary polygonal domains based on Sc hwarzChristoffel conformal mapping," International Journal for Nume rical Methods in Engineering vol. 80, pp. 103134, 2009. [43] D. A. Lowther, C. B. Rajanathan, and P. P. Silvester, "FiniteElement Technique for Solving 2D Open Boundary Problems," IEEE Transactions on Magnetics vol. 14, pp. 467469, 1978. [44] J. R. Brauer, "Open Boundary FiniteElements for Axisymmetric Magnetic and Skin Effect Problems," Journal of Applied Physics vol. 53, pp. 83668368, 1982. [45] C. Antunes, "Approxim ate Ballooning Techniques," IEEE Transactions on Magnetics vol. 19, pp. 25552557, 1983. [46] C. Antunes, E. M. Free man, D. Lowther, and P. Silves ter, "A Static Ballooning Technique for 2D Open BoundaryProblems," Journal of Applied Physics vol. 53, pp. 83608362, 1982. [47] Y. S. Sun, G. S. Zhang, and W. Chen, "A DataBase Method for 3dOpen Boundary Field Computation," IEEE Transactions on Magnetics vol. 26, pp. 807810, 1990. [48] W. Chen, G. S. Zhang, Y. S. Sun, M. J. Chen, and Y. Zhao, "The Application of the DataBase Method for 3D Open Boundary Field Computation," IEEE Transactions on Magnetics vol. 28, pp. 16861689, 1992. PAGE 162 162 [49] J. R. Brauer, S. M. Schaefer, J. F. Lee, and R. Mittra "Asymptotic BoundaryCondition for 3Dimensional M agnetostatic FiniteElements," IEEE Transactions on Magnetics vol. 27, pp. 50135015, 1991. [50] Q. S. Chen, A. Konrad, and S. Bar onijan, "Asymptotic BoundaryConditions for Axisymmetrical FiniteElement Electrostatic Analysis," IEEE Transactions on Magnetics vol. 30, pp. 43354337, 1994. [51] Q. Chen, A. Konrad, and P. P. Biringer, "Computati on of 3Dimensional Unbounded EddyCurrent Problems Using Asymptotic Boundar yConditions," IEEE Transactions on Magnetics vol. 31, pp. 13481351, 1995. [52] S. Gratkowski, L. Pi chon, and H. Gajan, "Asymptotic boundary conditions for open boundaries of axisymmetric magnet ostatic finiteelement models," IEEE Transactions on Magnetics vol. 38, pp. 469472, 2002. [53] R. L. Ungless, "An infinite finite el ement," vol. Master of Science. Canada: The University of British Columbia, 1973. [54] P. Bettess, "I nfinite Elements," International Journal for Numerical Methods in Engineering vol. 11, pp. 5364, 1977. [55] P. Bettess, "More on Infinite Elements," International Journal for Numerical Methods in Engineering vol. 15, pp. 16131626, 1980. [56] P. Bettess, Infinite elements Sunderland, UK: Penshaw Press, 1993. [57] B. M. A. Rahman and J. B. Davies "FiniteElement Analysis of Optical and Microwave Waveguide Problems," IEEE Transactions on Microwave Theory and Techniques vol. 32, pp. 2028, 1984. [58] M. J. McDougall and J. P. Webb, "Infinite Elements for the Analysis of Open Dielectric WaveGuides," IEEE Transactions on Microwave Theory and Techniques vol. 37, pp. 17241731, 1989. [59] M. S. Towers, A. McCowen, and J. A. R. Macnab, "Electromagnetic Scattering from an Arbitrary, Inhomogeneo us 2D Object a Finite and Infinite Element Solution," IEEE Transactions on Antennas and Propagation vol. 41, pp. 770777, 1993. [60] G. Beer and J. L. Meek "Infinite Domain Elements," International Journal for Numerical Methods in Engineering vol. 17, pp. 4352, 1981. [61] O. C. Zienkiewicz, C. Emson, and P. Bettess, "A Novel Boundary Infinite Element," International Journal for Numerical Methods in Engineering vol. 19, pp. 393404, 1983. [62] T. T. AbdelFattah, H. A. Hodhod, and A. Y. Akl, "A novel formulation of infinite elements for static analysis," Computers & Structures vol. 77, pp. 371379, 2000. [63] M. V. K. Chari and S. J. Salon, Numerical methods in electromagnetism San Diego: Academic press, 2000. [64] R. Arumugam, D. A. Lo wther, R. Krishnan, and J. F. Lindsay, "MagneticField Analysis of a Switched Reluctance Motor Using a TwoDimensional FiniteElement Model," IEEE Transactions on Magnetics vol. 21, pp. 18831885, 1985. [65] P. D. S. Periyasamy, "Finite element anal ysis of shell like structures using implicit boundary method," in Mechanical and Aerospace engineering vol. Master of Science. Gainesville: Universi ty of Florida, 2009, pp. 86. [66] S. Zhang and A. V. Kumar, "Magnetic Field Computation us ing Implicit boundary Finite Element Method," IEEE Transactions on Magnetics in press, 2010. PAGE 163 163 [67] S. Zhang and A. V. Kumar, "3D Magnet ostatic Analysis using Implicit Boundary Finite Element Method," Computers & Structures submitted 2010. PAGE 164 164 BIOGRAPHICAL SKETCH SungUk Zhang was born in Seoul, Sout h Korea. He graduated with a bachelors degree in electrical engineering from the S ogang University, South Korea in February 2003. In December 2006, he achieved a ma sters degree in bi omedical engineering from the University of Flor ida. Then, he enrolled in the doctoral problem in mechanical and aerospace engineering at the University of Florida in January 2007. His areas of interest include digital signal processing, electrical impedance tomography, inverse problem, finite element method, and computeraided geometric design. 