PAGE 1 1 COUPLED MAGNETO-ELASTOSTATIC AN ALYSIS USING IMPLICIT BOUNDARY FINITE ELEMENT METHOD By SUNG-UK 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 Sung-Uk 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............................................................................24Maxwell’s 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...........................................................38B-spline 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 7-1-1: 2D Coaxial C able....................................................................81Example 7-1-2: 2D Clapper Solenoid Actuator.................................................85Example 7-1-3: 2D Clapper Solenoid Actuator with Artificial Damage.............87Example 7-1-4: 2D Swit ched Reluctanc e Moto r...............................................88Three Dimensional Magnet ostatic Pr oblems..........................................................91Example 7-2-1: Iron Block in a Homogenous Magnet ic Fiel d...........................92Example 7-2-2: 3D Coaxial C able....................................................................96Example 7-2-3: 3D Plunger Solenoid Actuator.................................................99Example 7-2-4: 3D Clapper Solenoid Ac tuator...............................................104Magnetostatic Problems wit h Permanent M agnets...............................................109Example 7-3-1: U-S haped Permanent Magnet...............................................109Example 7-3-2: Three Dimensio nal Cylindrical Magnet ..................................113Magnetostatic Problems with Op en Boundary Tec hniques...................................115Example 7-4-1: U-Shaped Perm anent Magnet with Open Boundary Techniques..................................................................................................115Example 7-4-2: 2D Clapper Sole noid Actuator with Open Boundary Techniques..................................................................................................118Example 7-4-3: Force Between Two Parallel Wires.......................................120Coupled Magneto-Elastost atic Probl ems..............................................................123Example 7-5-1: 2D Clapper Solenoid Actuator with Cant ilever Beam............124Example 7-5-2: 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 Magneto-Elastostatic 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 6-1 Comparison of different open boundary techni ques [40]...................................688-1 Comparison for three clapper solenoid actuators ( NI =30)...............................1378-2 Comparison for three plunger solenoid actuators ( NI =30)...............................1418-3 Comparison for three combined plunger & clapper solenoid actuators ( NI =30)............................................................................................................. 1438-4 Comparison for coil actuators ( NI =30).............................................................146 PAGE 9 9 LIST OF FIGURES Figure page 1-1 Confo rmal me sh................................................................................................161-2 Conformal mesh ve rsus.....................................................................................181-3 3D structured gird s on multi-ma terial .................................................................181-4 Flapping wings operated by magnetic actuator.................................................212-1 Maxwell’s equati ons..........................................................................................252-2 Block diagram of a magnetic ac tuator...............................................................282-3 Reluctance method (Magnet ic circuit method)...................................................312-4 Fringi ng flux .......................................................................................................323-1 The solution structure with t he essential boundar y condit ion............................373-2 One dimensional field so lution using tw o element s...........................................403-3 Displacement surfac e plots in y-com ponent...................................................... 434-1 2D magnetosta tic prob lem.................................................................................444-2 2D magnetostatic prob lem with perm anent m agnet.......................................... 504-3 Two grids for two materi als................................................................................514-4 Component plots of interface soluti on structure at the material boundary.........544-5 Interface solution struct ure at the mate rial boundar y.........................................554-6 Three parts with own gr ids.................................................................................564-7 Contac t pairs .....................................................................................................564-8 Four element s in Part 1.....................................................................................574-9 Material in terface boundary ...............................................................................595-1 Analysis doma in and boundar ies.......................................................................635-2 Sequential analysis fo r 3D Magnetos tatics........................................................676-1 Two-dimensional interior region with an exte rior annul us..................................69 PAGE 10 10 6-2 Structured grid and rectangular boundary.........................................................726-3 Coordinate transformation for the decay function infi nite element.....................756-4 1D shape functions using exponential decay functi ons.....................................776-5 Structured grids with the infinite elements.........................................................786-6 Domain of integrat ion........................................................................................796-7 Solution plots for 1D in finite element example...................................................807-1 2D coaxial cable model with structur ed grid s.....................................................827-2 Magnitude of H fi eld...........................................................................................827-3 Convergence pl ot for H1 norm...........................................................................837-4 Convergence pl ot for L2 norm...........................................................................847-5 2D Clapper solenoi d with struct ured gr id...........................................................857-6 Magnetic vector potentia l and magnetic fl ux li nes.............................................867-7 Magnetic force ve rsus air gap length.................................................................877-8 Clapper solenoid wit h artificial damage.............................................................887-9 2D Planar model of s witched reluctanc e motor.................................................887-10 Magnetic vector potent ial and magnetic flux lines in the aligned position..........907-11 Magnetic vector potent ial and magnetic flux lines in the unaligned position......917-12 Iron cube in homogeneous magnetic field.........................................................927-13 Iron objects with t he same grid density..............................................................937-14 Cross-sections wit h the line y= z=10mm............................................................947-15 Components of B al ong the line y= z=10mm......................................................957-16 3D coaxial cable model with the stru ctured gr id................................................967-17 Magnitude of H field for 3D coaxia l cabl e..........................................................987-18 Magnetic field in the hood direction vers us. radi us.............................................987-19 Convergence pl ot for H1 norm...........................................................................99 PAGE 11 11 7-20 Plunger solenoid actuator of axisymmetric geometry......................................1007-21 Top view of two plunger actuat ors...................................................................1017-22 3D solid model of solenoid ac tuators with plunger armature s..........................1017-23 The magnetic flux density in the ydirection for two pl unger solenoids............1027-24 The magnetic field in the y-dire ction for two plung er solenoi ds.......................1037-25 Magnetic force versus gap length for plunger solenoids ..................................1047-26 Clapper solenoid actuator of axisymmetric geometry......................................1057-27 Top view of two clapper actuat ors...................................................................1057-28 3D solid model of solenoid ac tuator with clapper armatures...........................1067-29 The magnetic flux density in the ydirection for two cl apper actuators.............1077-30 The magnetic field in the y-dire ction for two clapper actuators........................1077-31 Magnetic force versus gap length for clapper solenoids ..................................1087-32 U-shaped permanent magnet model ...............................................................1097-33 Surface and contour plots for magnetic potential fr om Comsol.......................1107-34 Surface and contour plots for magnetic potential fr om IBFEM........................1117-35 The observation line for the com parison ..........................................................1117-36 Magnetic field in the x-direction on the line ......................................................1127-37 Magnetic field no rm on the lin e.........................................................................1127-38 3D solid model for cylindr ical magnet in the air................................................1137-39 Magnetic field in the y-direct ion........................................................................1147-40 Magnetic field along y ax is...............................................................................1157-41 U-shaped permanent magnet with the reduc ed air domai n..............................1167-42 Contour plots for magnet ic vector pot ential......................................................1177-43 Magnetic field in the x-dire ction on the obser vation lin e...................................1177-44 2D Clapper sol enoid actuat or...........................................................................118 PAGE 12 12 7-45 Contour plots of magnetic vector potentia l.......................................................1197-46 Magnetic force ve rsus gap l ength.....................................................................1207-47 Two long para llel wires .....................................................................................1217-48 Contour plots for magnet ic vector pot ential......................................................1227-49 The magnetic force versus t he distance between two wires .............................1237-50 Planar clapper solenoid act uator with a cant ilever bea m..................................1247-51 Displacement at the tip versus NI with varying attach ment location.................1257-52 Planar clapper solenoid act uator with a cant ilever beam ..................................1257-53 Displacement on the tip versus NI with varying beam thickness......................1267-54 Top views of structures atta ched to the plunger armature................................1277-55 3D plunger solenoid act uators with A) Solid plate, B) plate with one hole, and C) plate with two holes .....................................................................................1287-56 Deformation due to magnetic fo rce...................................................................1297-57 Maximum displacement versus NI.................................................................... 1308-1 Solenoid actuator wit h clapper arma ture..........................................................1338-2 Computed results using IBFE M........................................................................1358-3 Clapper solenoi d actuator s...............................................................................1358-4 Contour plots of magnetic vector potential for cla pper solenoid actuators........1368-5 Solenoid actuator wit h plunger arma ture..........................................................1378-6 Computed results using IBFE M........................................................................1398-7 Plunger solenoi d actuator s...............................................................................1408-8 Contour plots of magnetic vector potential for pl unger solenoid actuators.......1418-9 Three combined plunger & cl apper solenoid ac tuators.....................................1428-10 Contour plots of magnetic vector potential fo r combined plunger & clapper solenoid act uators............................................................................................1438-11 Three coil ac tuators..........................................................................................144 PAGE 13 13 8-12 Magnitude of B fields for three coil ac tuators....................................................1458-13 The best magnetic actuator among several designed actuators.......................1468-14 Sold model of the 3D coil act uator.................................................................... 1478-16 The magnitude of the magnetic flux density of the coil act uator.......................1488-17 Lorentz force of the coil actuator versus NI...................................................... 1498-18 The coil actuator with flapping wings................................................................1498-19 Four structures with flapping wings..................................................................1508-20 Four surface structur es with flappi ng wings ......................................................1518-21 Structured mesh and boundary conditions of t he first de sign...........................1528-22 Displacement in the z-direct ion during the wi ng stroke .....................................1538-23 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 MAGNETO-ELESTOSTATIC AN ALYSIS USING IMPLICIT BOUNDARY FINITE ELEMENT METHOD By SUNG-UK 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 o-elastostatic 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 magneto-elastostatic 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 1-1 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 1-1. 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 well-placed 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) [8-10] 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 X-FEM. 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 1-2 B, the structured mesh, such as the examples shown in Figure 1-2 C, is ea sy to generate since all elements are regular shaped and the grid does not have to conform to the geometry. Figure 1-2. 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 1-3. Wi thin overlapping elements at the interface, the piece-wise interpolation wit hin each grid in combined into a single solution structure as explained in the later chapter. Figure 1-3. 3D structured girds on multi-mate rial. A) Material 1, B) Material 2 and C) Multi-material PAGE 19 19 IBFEM is extended to solve magneto-st atic problems. Some magneto-static 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 magneto-static 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 Poisson’s equation. The Poisson’s 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 B-spline shape functions instead of the traditional finite element shape functions if needed. The advantage of using B-splines is that the computed magnetic flux density and field (or stresses and strains) are continuous between elements. Under the 2D assumption, magneto-static a nalysis is computationally inexpensive; however, 3D magneto-static 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 magneto-static analysis, in 3D, the current density becomes a vector field. Therefore, the curr ent density distribution must be computed by electrostatic analysis prior to magneto-stat ic analysis. Under the above considerations, several 3D magneto-static 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 magneto-static anal ysis, the method is extended to solve coupled magneto-elastostatic 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 lane-like fixed wings and helicopter-like 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 magneto-elastostatic 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 magneto-elastic 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 1-4. Figure 1-4 A shows the downward wing stroke when the magnetic actuator turns off. Figure 1-4 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 1-4. 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 multi-ma terial system analysis and to perform magnetostatic analysis as well as c oupled magneto-elasto static analysis. The main objectives of this disse rtation are to extend IBFEM to Perform 2D magneto-static analysis Perform 3D magneto-static analysis Model permanent magnets PAGE 22 22 Perform multi-material analysis Compute magnetic forces, and Solve coupled magneto-elastostatic analysis In addition, IBFEM are tested and validated us ing several examples: coaxial cable, solenoid actuator, U-shaped 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, Maxwell’s equations, overview of magnetic actuators and ov erview of reluctance method. Among the Maxwell’s equations, Gauss’s law and Am pere’s 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 B-spline element. Chapter 4 discusses 2D magnet o-static analysis with IBFEM. A modified solution structure is introduced for t he multi-material 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 Maxwell’s equations. By the Maxwell’s 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 Maxwell’s equations was using elaborate ma thematics such as series expansions, Legendre polynomials, Bessel’s 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 Maxwell’s 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. Maxwell’s Equations Maxwell’s 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 Gauss’s law, Gauss’s law for magnetism, Faraday’s law of induction and Ampere’s law. Figure 2-1 shows concept ual drawings for Maxwell’s equations Figure 2-1. Maxwell’s equations. A) Gauss’s law, B) Gauss’s law for magnetism, C) Faraday’s law of induction, and D) Ampere’s circuital law. Figure 2-1 A shows Gauss’s law that relates electrical charge within a close surface to the surrounding electric field. The differential form of the Gauss’s law is expressed as D (2-1) PAGE 26 26 where, D is the electric flux density and is the electrical charge density. The constitutive equation is written as DE (2-2) where, E is the electric field, and is the permittivity for diel ectric material. Figure 2-1 B shows Gauss’s 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 Gauss’s law for magnetism is stated as 0 B (2-3) where, B is the magnetic flux density. Figure 21 C shows Faraday’s 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 Faraday’s law of induction is stated as t B E (2-4) Ampere’s law with Maxwell’s correction is show n in Figure 2-1 D. The law indicates that two factors can generate magnetic field; electrical current and changing electric flux density. The differential form for Ampere’s la w with Maxwell’s correction is written as t D HJ (2-5) where, H is the magnetic field. The constitutive equation between H and B can be stated as BH (2-6) 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 (2-7) After taking the divergence of Ampere’ s law, the equation can be written as 0 tt D D HJJ (2-8) Equation 2-8 can be restated as t J (2-9) 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 (2-10) 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 electro-hydraulic 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 2-2 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 2-2. 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 2-2. 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 current-carrying coil is used for coil actuators. The Lorentz force is expressed as FNIBl (2-11) 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 actuator’s 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 IBFEM’s results. The reluctance method is based on the Amper e’s law. Ampere’s law in integral form can be expressed as a summation form as follows k kkk kk kB NIHllNI Hdl (2-12) where, NI is the ampere-turns, 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 (2-13) 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 cross-section surface area kS. Substituting into the Ampere’s law in the summation fo rm, Equation 2-12 is stated as kk k kkl NI S (2-14) PAGE 31 31 As the divergence of flux density is zero, the fluxes through all segments have a same quantity So Equation 2-14 is rewritten as k k R NI (2-15) 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 (2-16) The equation is similar to the familiar Ohm’s 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 2-3. Figure 2-3. 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 2-4. PAGE 32 32 Figure 2-4. 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 air-gap 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 2-4. 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 non-linear 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 Kronecker’s 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 non-banded 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 (X-FEM) [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 R-function method is a way to define implicit equations for domains using Bool ean operations between simple primitives. Shapiro and Tsukanov extended the R-function method to so lve non-stationary physical problems with time-varying geometries, and us ed transfinite Lagrangian interpolation to apply essential boundary conditions [21]. Hol lig et al. proposed the Web-method which uses weighted extended B-spline basis to guar antee higher order continuous solutions PAGE 35 35 [22]. Their method has been used with R-functi ons or distance functions as implicit equations in order to construct the solution structures. In addition, Apaydin et al have extended the Web-method to solve one-dimensional 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 B-spline 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 (3-1) where u is the displacement vector g u is the grid variable represented by a piece-wise 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 D-function) is defined as 0()0 () ()110() 1() x x xx xkD (3-2) where, k is an integer. The D-function 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 D-function 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 D-function is a good approximation of a step function. Hollig et al. have used a weighting function that is similar to the D-function. 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 (3-3) 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 3-1 graphically depicts the solution structure. The element e1 is a boundary element and the others are internal elements. Figure 3-1. 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 (3-4) For 2D problems in Voigt notation these strains can be defined as PAGE 38 38 ,TT ssssaaaa sauvuvuvuv xyyxxyyx (3-5) Using the constitutive equation, the stress is s asa C C (3-6) where, C is the elasticity matrix. The governing equation for linear elasticity is stated as 0in b (3-7) where, b is the boundary force, is the analysis domain. The essential boundary and natural boundary condit ions are defined as on,onut00uu nT (3-8) 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 (3-9) 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 (3-10) 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 (3-11) The virtual strain and stress are stated as ,ees BX CBX (3-12) 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 (3-13) 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. B-spline element Using Lagrange elements, the tr aditional FEM only warrants C0 continuity of a solution between elements. When B-spline approximation is used, higher order continuity of the solution can be guaranteed. Figure 3-2 shows one dimensional solution using Lagrange interpolation and one using B-spline 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 3-2 A. Figure 3-2. One dimensi onal field solution using two elements. A) Lagrange interpolation, and B) B-sp line approximation scheme. B-spline 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 B-sp 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 B-spline element. In IBFEM, uniform B-splines can be us ed, which means nodes have been equally spaced in the parametric s pace. The B-spline 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 B-spline 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 (3-14) 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 (3-15) where, i f and 1 i f are the approximated functions and iu are the support nodal values. As given B-spline elements warrants C0 and C1 continuity, two continuity requirements provides the following restricted conditions : 1(1)(1)iiff, and 1(1)(1)iiff rr (3-16) 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 (3-17) 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 (3-18) The B-spline element using these basis func tions is called Quadratic B-spline element (QBS). Using the same approach, t he B-spline basis functions with C2 continuity can also be obtained. The B-spline element with the basis functi ons is defined as Cubic Bspline element (CBS). Two and three dimensional B-spline elem ents are created by taking product of one dimensional B-spline elements. For exam ple, 2D quadratic B-spline element has nine supports nodes and nine basis functions 2D cubic B-spline 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 Young’s modulus, I is the moment of the inertia of the beam cross section, is the Poisson’s ratio, L is the length of the beam, t is the thickness of the beam. Assuming that the material is isotropi c, Young’s 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 3-3. Displacement surface plots in y-component. A) Using 4 node bilinear elements B) Using 9 node B-spline elements Figure 3-3 shows displacement contours us ing two element types. In the case of Figure 3-3 B, where quadrat ic B-spline 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 (4-1) 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 (4-2) where, Bis the magnetic flux density and H is the magnetic field density. Figure 4-1 shows 2D assumption where, ˆ J xy Jk, and ˆ A xy Ak. That is, the current only flows in the z-direction 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 x-y. Figure 4-1. 2D magnetostatic problem PAGE 45 45 If the material is homogeneous, Equation 4-1 is expressed as the gradient of the scalar function A as following 111 AA A J xxyy (4-3) The constitutive equation for the 2D problem can be rewritten as 11 ˆˆ A A yx HBij (4-4) The solution structure for the magnetic vector potential A is defined as ()()()()()()gasaADAAAA xxxxxx (4-5) where, g A is a grid variable that is defined by piece-wise interpolation or using B-spline 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 4-3, the weak form for 2D magnet ostatics is obtained as 1 1()()()()()()sssssa t A AdAJdAHAAd (4-6) 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 (4-7) 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 (4-8) Using Equations 4-7 and 4-8, 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 (4-9) 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 D-function which can have ve ry large values near the boundary. They are expressed as 1 1[]i ij jN BD x B (4-10) 2 2[]iji j D BN x B (4-11) The element matrix to be assembled into the global equations can be obtained as 1 123eT eeee ed KBBKKK (4-12) 1 111eT e ed KBB (4-13) 1 222eT e ed KBB (4-14) 11 31221eTT e ed KBBBB (4-15) As 2 B has terms containing derivatives of the D-function, it is non-zero 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 non-zero 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 (4-16) 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 (4-17) where, 2 01k kD d x (4-18) 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 (4-19) where, 01k kD Dd x (4-20) 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 (4-21) 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 (4-22) where, 01r is the reluctivity. If the value of the reluctivit y becomes constant, then the governing equation is restated as 00AMJ (4-23) The equation becomes a nonlinear Poisson’ s equation for the vector potential. Multiplying the weight function A on both sides and integrating, Equation 4-23 becomes 00() ddAMAJA (4-24) Using Green’s first identity for vector fields, the left hand side term becomes 00 0000() ()()d dd AMA AMAAMA (4-25) Using the divergence theorem, the se cond term on the right hand side becomes PAGE 50 50 0000()() ddAMAAMAn. (4-26) Using the identities of FGGF and F GTFGT Equation 4-26 becomes 0000()()ddAMAnAAMn (4-27) The weak form becomes 0000()d ddd AA AMAAMnJA (4-28) 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 (4-29) 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 4-2. Figure 4-2. 2D magnetostatic problem with pe rmanent magnet PAGE 51 51 The 2D weak form becomes 00AAAAdddMJ (4-30) 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 (4-31) 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 (4-32) 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 4-3. Figure 4-3. 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 (4-33) 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 4-33 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 D-function. As explained earlier for elements on the external bounda ries, the gradient of the D-function is zero outside a narrow band near the interface boundaries. Ther efore, the volume integral for terms containing the gradient of t he D-function 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 multi-material 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 (4-34) 12nAnA (4-35) 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 4-35) 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 4-33 combines the interpolation within the overlapping elements at an interface to represent the solution near the interface. When the D-function in unity, the solution is given by 2 g A A, which is the solution from the second grid and when the D-function is zero, the solution is given 1 g A A, which is the solution from first grid. In the region where the D-function 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 (4-36) 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 D-function is zero for 0 and non-zero 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 4-4 graphically shows plots of components for the interface solution structure at the material interface boundary. Figure 4-4 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 3-7 C and D represents A field distributions from the grid 1 and the grid 2. A B C D Figure 4-4. 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 4-5. Interface solution st ructure at the material boundary Figure 4-5 shows the interface solution st ructure at the mate rial boundary following by Figure 4-4. The x-axis an d the y-axis 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 4-6. For t he first part (Figure 4-6 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 4-6 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 4-6. 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 4-7. A B Figure 4-7. 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 4-8 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 4-8. Four el ements in Part 1 In order to express severa l interface boundaries among mult i-materials, the D-function 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 (4-37) As the element e2 has the 2nd contact pair, the interface solution structure becomes 3 113131 g g g A DADA (4-38) where, 13D is the D-function 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 (4-39) 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 (4-40) 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], Maxwell’s 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 (4-41) 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 (4-42) where, B is the magnetic flux density. In a linear, isotropic and non-compressible ferromagnetic material, 21 2fH f (4-43) The virtual work done by the force densities can be evaluated as follows: 21 2mdHddf u uJB u (4-44) 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 (4-45) 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 (4-46) where, ˆn is the direction normal to the boundary between the two materials with different permeability as shown in Figure 4-9. Figure 4-9. Material interface boundary: 1 (Inner material permeability) and 2 (Outer material permeability) PAGE 60 60 Using Equation 4-46, the force on the fe rromagnetic material is expressed as 2 1011 ˆ () 22 1 ˆ () 2nddnd n dd T uHHHHnu HHnu (4-47) 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 Gauss’s 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 (4-48) 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 4-47 is rewritten as 1 22 222 21 2 2111111 ˆˆ (()) 222n ttnB HddHBd nuun (4-49) 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 (4-50) 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 4-50 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 (4-51) 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 (5-1) 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 (5-2) H0 on nA (5-3) PAGE 63 63 Figure 5-1 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 5-1. Analysis domain and boundaries The weak form for these governing equatio ns and boundary conditions, obtained using the weighted residual method, is HBdddAAAHnJA (5-3) 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 (5-4) 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 5-4) into the weak form (Equation 5-3), a modified weak form of the 3D m agnetostatic equation can be derived as ssssadddAAJAAA (5-5) 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 (5-6) 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 (5-7) for 1,iZN. The curl of s A can be computed as [] ABA s Cg (5-8) 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 D-function. 12 BBBCCC, where 111121 221222... ... BBBB BBBBCCCC n CCCC n (5-9) 2233 32 11133 31 1122 210 0 0 Bii C ii i iiNN DD x x NN DD x x NN DD xx (5-10) 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 (5-11) 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 (5-12) 111eT eCC ed KBB (5-13) 222eT eCC ed KBB (5-14) 31221eTT eCCCC ed KBBBB (5-15) Since 2 C B contains the derivates of D()xii, it is non-zero 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 (5-16) 31221 01 KBBBBeTT eCCCC edd (5-17) 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 three-dimensional space, magnetic force computation and multi-material 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 5-2 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 Laplace’s equations, the solution structure and finite element im plementation for 3D electrosta tics is identical to ones in the previous 2D magnetostatics. Figure 5-2. 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 6-1 [40]. Unfortunately, no all-purpose powerful techni que exists, so a researcher should choose one of these according to their application. Table 6-1. 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 1st-order 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 6-1. 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 6-1. Two-dimensional 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 well-established mesh. Additionally, the me thod is only applicable to 2D Laplace’s 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 quasi-static 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 (6-1) Using Equation (6-1), the Ampere’s law can be restated as 2211 BAAAAJ (6-2) where, 0 A that is the Coulomb gauge condition. When A=0 at infinity, the general solution of the Poisson’s equation is (') ()' 4 d Jr Arv rr (6-3) 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 (6-4) 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 (6-5) where, (,) a is a function in the spherical coor dinate system. The leads to following equation: 2 0rr A A (6-6) which is called the first-order 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 1-2 with the unit normal of (1,0,0) is shown in Figure 6-2. Figure 6-2. Structured gr id and rectangular boundary The derivative of A in the radial direction on the line segment becomes cossinzzzzzAAAAA xy rxryrxy (6-7) Substituting Equation 6-7 into Equation 6-6, the equation becomes 02 0zz zx AA y A xryrr (6-8) PAGE 73 73 Since 0cos x r and sin y r. Therefore, ABC on a line s egment 1-2 in Figure 6-2 can be written as 002zz zAA y A x xxy (6-9) 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 (6-10) Similarly, the line segments with the unit no rmal of (0,1,0) can have the following ABC 002zz zAA x A y yyx (6-11) The coefficients of t he boundary matrix become 02ej i ijijijeN N KNNxNNd yxx (6-12) 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 Guass-Lagueree quadra ture instead of Gauss-Legendre 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]. Abdel-Fattah 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 Gauss-Legendre 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 (6-17) where, (,)iNst are the shape functions on the infinite domain and(,)i f st are decay functions. Figure 6-3 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 6-3. 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, Gauss-Laquerre integration is used instead of Gauss-Lagendre 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 6-3. 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 (6-18) 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 (6-19) Figure 6-4 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 6-4 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 (6-20) 2exp(1/) A B L r (6-21) After eliminating A and B, the dec ay parameter is estimated using 1r and 2r as follow 1 22lnr L r (6-22) A B C D Figure 6-4. 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 6-5 shows structured grids with the infinite elements. Figure 6-5. 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 6-6 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 6-6. Domain of integration The differential equation is 2 20 dP dx (6-23) For the infinite element, the mapping func tion and the shape function are defined as 11 1, 2 s Ms (6-24) 1/ 11 1, 2sLs Nes (6-25) As 0P at infinity, 11K, coefficient of the infinite el ement matrix, is stated as follow 22/ 11 1111 11()sLNN dx Kdsepsds xxds (6-26) where, 2 1 1111 () M psM sL and 1 dx ds Substituting 1 2 L ss the coefficient becomes 111111 00()() 2ssdsL Kepsdsepsds ds (6-27) The integration is done using Guass-Laguer ee integration wit h following modified weights and abscissas 2/2L newoldL WWe (6-28) 1 2 L ss (6-29) PAGE 80 80 A solution for this example is provided by Chari and Salon [63]. Figure 6-7 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 6-7. 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 7-1-1: 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 cross-section is created. Figur e 7-1 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 z-direction) and the current density is assumed to be uniform. PAGE 82 82 Figure 7-1. 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 (7-1) 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 7-2. Magnitude of H field PAGE 83 83 Figure 7-2 shows the magnitude of the magnetic field t hat was computed using the quadratic B-spline 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 7-3. 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 7-4. Convergence plot for L2 norm Figure 7-3 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 (7-2) where, eB is the exact value of the magnetic fl ux from the analytical solution and hB is the corresponding computed value. Fig. 7-4 shows convergence of L2 error norm of A(x,y) This error norm is defined as 1 2 2T ehehLAAAAd (7-3) 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 B-spline elements but the rate of convergence decreases with increasing number of nodes. PAGE 85 85 Example 7-1-2: 2D Clapper Solenoid Actuator The clapper armature soleno id is composed of an armature, a stator, and coils. A two-dimensional planar model is shown in Fi gure 7-5 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 7-5 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 7-5. 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 B-spline elements was used for the results shown in Figure 7-6 B. E ssential boundary conditions are applied along the axis of symmetry to enforce symmetry. Figure 7-6. Magnetic vector potential and magnetic flux lines. A) Comsol, and B) IBFEM Figure 7-6 shows flux lines and the plot of magnetic vector potential when the width of the air gap is 2mm. Figure 7-6 A shows t he result using commercial FEA software (COMSOL). The result from IBFEM is shown in Figure 7-6 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 7-7 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 Maxwell’s 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 7-7. Magnetic force versus air gap length Example 7-1-3: 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 7-8 A. As in the previous exampl e, air gap width equal to 2mm was used to compute the results shown in Figure 7-8 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 7-8 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 7-8. Clapper solenoid with artifi cial damage. A) Da maged region and boundary condition, and B) Flux lines and surface plot of A Example 7-1-4: 2D Switched Reluctance Motor A 2D planar model of the Switched Reluct ance Motor (SRM) is shown in Figure 7-9. 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 7-9. 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. B-spline elements ar e used for the analysis so that accurate results are obtained even with a sparse grid as shown in Figure 7-10 and Figure 7-11. 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 7-10. Magnetic vector potential and magnetic flux lines in the aligned position PAGE 91 91 Figure 7-11. 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 7-2-1: 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 7-12 shows one-eighth 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 z-direction. Figure 7-12. Iron cube in homogeneous magnetic field PAGE 93 93 The half-length 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 z-direction 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 (7-4) H0 on (z=0 and z=b) nA (7-5) As the magnetic flux density only exists in the z-direction, 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 z-direction. Figure 7-13. 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 7-13. Figure 7-13 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 7-14. Figure 7-14. Cross-sections 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 7-15. Figure 7-15 A shows that magnetic flux densit y in the x-direction is continuous only when the shape of the iron is cube because only in this case the normal component is parallel to the x-direction. 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 7-15. Components of B along the line y=z=10mm. A) Bx, B) By and C) Bz PAGE 96 96 Example 7-2-2: 3D Coaxial Cable A coaxial cable consists of an inner conductor, an insulator, and an outer conductor. A coaxial cable problem is a well-known problem for 2D magnetostatic problem. As the current density and the m agnetic vector potential has only one component in the z-direction, the analytical solution is easily obtained using the ampere’s law. Using this example, we tri ed to solve an electro-magnetostatic 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 (7-6) 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 7-16 shows the coaxial cable model using thr ee separate structured grids. Figure 7-16. 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 z-direction) and the current density is assumed to be uniform. The anal ytical solution of the magnetic field in circumferential direction can be Equation 7-1 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 7-17 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 7-18 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 7-19 shows the convergence pl ot for H1 norm using 8 node hexahedral elements. H1 norm is defined as Equation 7-2. The H1 norm decreases as the larger number of elements is used. PAGE 98 98 Figure 7-17. 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 7-18. Magnetic field in the circumferential di rection versus. radius PAGE 99 99 103 105 Log(Num of nodes)Log(H1 Norm) Figure 7-19. Convergence plot for H1 norm Example 7-2-3: 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 7-20. Plunger solenoid act uator of axisymmetric geometry Figure 7-20 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 7-20. 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 7-21. PAGE 101 101 Figure 7-21. Top views of two plunger actuators. A) Cylindr ical plunger, and B) Brick plunger Figure 7-21 A is the top view of the cylinder actuator. Figure 7-21 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 7-22. 3D solid model of solenoid act uators with plunger armatures. A) Cylindrical plunger, and B) Brick plunger. PAGE 102 102 Figure 7-22 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 y-direction. 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 7-23. Figure 7-23. The magnetic flux density in the y-direction 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 y-direction is shown in Figure 7-24. Figure 7-24. The magnetic field in the y-di rection for two plunger solenoids. A) Cylindrical plunger, and B) Brick plunger Figure 7-24 A shows the computed magnetic fi eld in the y-direction 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 7-25 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 7-25. Magnetic force versus gap length for plunger solenoids Example 7-2-4: 3D Clapper Solenoid Actuator One of other popular solenoid actuators is a clapper sol enoid actuator. Figure 7-26 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 7-26. 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 7-27. Figure 7-26. Clapper solenoid act uator of axisymmetric geometry Figure 7-27. Top views of two clapper actuators. A) Cylindr ical clapper, and B) Brick clapper. PAGE 106 106 Figure 7-27 A is the top view of the cylinder actuator. Figure 7-27 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 7-28. 3D solid model of solenoid act uator with clapper armatu res. A) Cylindrical clapper, and B) Brick clapper Figure 7-28 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 y-direction. 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 7-29. Figure 7-29. The magnetic flux density in t he y-direction 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 y-direction is shown in Figure 7-30. Figure 7-30. The magnetic field in the y-di rection for two clapper actuators. A) Cylindrical clapper, an d B) Brick clapper PAGE 108 108 Figure 7-30 A shows that the computed magnetic field in the y-direction 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 7-31. Magnetic force versus gap length for clapper solenoids Figure 7-31 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 U-shaped 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 7-3-1: U-Shaped Permanent Magnet A U-shaped permanent magnet is analyz ed using IBFEM and compared with a commercial software (Comsol). The drawing of the model is shown in Figure 7-32 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 x-direction. 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 7-32. U-shaped 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 7-32 B. The density of conf ormal mesh is denser near the permanent magnet. Figure 7-33 and Figure 7-34 show surfac e plots and contour plots for magnetic vector potential. Figure 7-33 is the result from Comsol and Figure 7-34 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 7-33. Surface and contour plot s for magnetic potential from Comsol PAGE 111 111 Figure 7-34. Surface and contour plot s for magnetic potential from IBFEM In order compare two results precisely, t he magnetic field in the x-direction and the magnetic field density norm are obtained along a line. The line for comparison is shown in Figure 7-35. Figure 7-36 shows the magnet ic field in the x-direction. 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 7-35. 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 7-36. Magnetic field in the x-direction 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 7-37. Magnetic field norm on the line PAGE 113 113 Example 7-3-2: Three Dimens ional Cylindrical Magnet A 3D cylindrical permanent magn et in free space is modeled as shown in Figure 7-38. Considering the symmetry, one fourth of the permanent magnet is created in order to reduce computation. Figure 7-38. 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 y-direction. U nder the conditions, the analytic al solution of the magnetic field along the y axis is given by [2]: 02yM BB H AA (7-6) where, 221 2 Rz A dd 1 2 z B d and 0 if /2 2 if