Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0041475/00001
## Material Information- Title:
- Coupled Magneto-Elastostatic Analysis Using Implicit Boundary Finite Element Method
- Creator:
- Zhang, Sung-Uk
- Place of Publication:
- [Gainesville, Fla.]
- Publisher:
- University of Florida
- Publication Date:
- 2010
- Language:
- english
- Physical Description:
- 1 online resource (164 p.)
## Thesis/Dissertation Information- Degree:
- Doctorate ( Ph.D.)
- Degree Grantor:
- University of Florida
- Degree Disciplines:
- Mechanical Engineering
Mechanical and Aerospace Engineering - Committee Chair:
- Kumar, Ashok V.
- Committee Members:
- Schueller, John K.
Ifju, Peter Sankar, Bhavani V. Nishida, Toshikazu - Graduation Date:
- 4/29/2010
## Subjects- Subjects / Keywords:
- Armatures ( jstor )
Boundary conditions ( jstor ) Clappers ( jstor ) Flux density ( jstor ) Magnetic fields ( jstor ) Magnetic flux ( jstor ) Magnetism ( jstor ) Magnets ( jstor ) Reluctance ( jstor ) Solenoids ( jstor ) Mechanical and Aerospace Engineering -- Dissertations, Academic -- UF fem, ibfem, magnetostatics - Genre:
- Electronic Thesis or Dissertation
born-digital ( sobekcm ) Mechanical Engineering thesis, Ph.D.
## Notes- Abstract:
- COUPLED MAGNETO-ELESTOSTATIC ANALYSIS USING IMPLICIT BOUNDARY FINITE ELEMENT METHOD 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 easier 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 magneto-elastostatic problems that typically involve an assembly of parts made of several different materials. The geometry is accurately modeled using equations from CAD models and a separate structured mesh is used for each part in an assembly. Specially constructed solution structures are used to represent test and trial functions such that boundary and interface conditions are enforced. Several magnetostatic problems with known solutions are modeled to validate the method. The magnetostatic problems are classified as unbound problems so that sometimes a very large analysis domain should be modeled to get more accurate results. In order to reduce the analysis domain, two open boundary techniques are developed for IBFEM: asymptotic boundary conditions and decay function infinite element. In addition, a magnetostatic problem with permanent magnets is solved using IBFEM. ( en )
- General Note:
- In the series University of Florida Digital Collections.
- General Note:
- Includes vita.
- Bibliography:
- Includes bibliographical references.
- Source of Description:
- Description based on online resource; title from PDF title page.
- Source of Description:
- This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
- Thesis:
- Thesis (Ph.D.)--University of Florida, 2010.
- Local:
- Adviser: Kumar, Ashok V.
- Statement of Responsibility:
- by Sung-Uk Zhang.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Zhang, Sung-Uk. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 697291052 ( OCLC )
- Classification:
- LD1780 2010 ( lcc )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

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 PAGE 114 114 Figure 7-39. Magnetic fiel d in the y-direction Figure 7-39 shows the magnetic field in t he y-direction. Eight node hexahedral elements are used and the total number of the elements is 20000. T he result graphically shows the discontinuity at the material interfac e. Figure 7-40 shows the line plots for the magnetic field along the y axis. The resu lt of Comsol was computed using a 2D axisymmetric model. As shown in the Figure, IBFEM result is close to the analytical solution and the comsol result. Some discrepan cy of the magnetic field is founded at the location where it is far from the interface boundary. The re sult may be more accurate when the number of elements increases, the quadratic or cubic el ements are used or an open boundary technique is applied. PAGE 115 115 0 0.2 0.4 0.6 0.8 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 y axisHy [A/m] Analytical Comsol IBFEM Figure 7-40. Magnetic field along y axis Magnetostatic Problems with Open Boundary Techniques Several examples with the open boundary techniques are examined here to validate the implicit boundary method. The first example is Example 7-3-1 with smaller analysis domain. The second is the entire model of the previous sol enoid actuator from Example 7-1-2. The third example is to compute Lorentz force between two wires. These examples show how the open boundary techniques could improve results. Example 7-4-1: U-Shaped Permanent Ma gnet with Open Boundary Techniques For the previous example (Example 73-1), the large unbounded domain (6 x 4) was necessary with the truncation approach. The small air domain (1.6 x 1.4) is created as shown in Figure 7-41. PAGE 116 116 Figure 7-41. U-shaped permanent magnet with the reduced air domain Two open boundary techniques (asympt otic boundary conditions and decay function infinite method) are applied on the outer boundaries. Figure 7-42 shows contour plots for magnetic ve ctor potential for each technique. The total number of four node bilinear elements used in the model is 8888. The decay factor1 L for the decay function is used in infinite element met hod. Figure 7-42 A show s different contour pattern when the homogenous essential bo undary is applied on the outer boundaries. Figure 7-43 shows magnetic fiel d in the x-direction on the observation line as shown in Figure 6-35. EBC stands for the homogenou s essential boundary conditions. The results using ABC are quite similar to the Comsol result. The result obtained using EBC has the largest error among them. PAGE 117 117 Figure 7-42. Contour plots fo r magnetic vector potential. A) Homogenous essential boundary conditions, B) Asymptotic boundary conditions and C) Decay function infinite element with L =1. 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 5 4 3 2 1 0 1 2 3 4 x 105 x axisHx Comsol EBC ABC Infinite Figure 7-43. Magnetic fi eld in the x-directio n on the observation line PAGE 118 118 Example 7-4-2: 2D Clapper Solenoid Ac tuator with Open Boundary Techniques The second example for open boundary techniqu es is extended from the previous solenoid actuator (Example 7-1-2). The whol e actuator including smaller air domain is created as shown in Figure 7-44. Three open boundary techni ques are applied on the other boundaries. The magnetomotive force NI and material properties of components are the same as in Example 7-1-2. Figure 7-44. 2D Clapper solenoid actuator When the gap length is 2mm, Figure 7-45 show s the contour plots for magnetic vector potential. Among the three resu lts, the result using the homogenous essential boundary conditions differs from the other two. The quadratic b-spline elemen ts are used for the analysis. The total number of elements is 7049. The computed forces for all open boundary techniques are as follows: Homogenous essential boundary conditions (EBC): 269 N/m Asymmetric boundary conditions (ABC): 276 N/m Decay function infinite element me thod with L=1 (Infinite): 282 N/m. The computed force using EBC is the lowest value. These computed forces can be also compared with the following va lues: the approximate solu tion using the reluctance PAGE 119 119 method is 245.2N/m and the force computed by FEM with the ballooning method is 272 N/m from [18]. The computed force using ABC is quite close to the reference force from [18]. Therefore, using ABC can produce the most reliable answer in three techniques. A B C Figure 7-45. Contour plots of magnetic vector potentia l. A) Homogenous essential boundary conditions, B) Asymptotic boundary conditions and C) Decay function infinite element with L =1. PAGE 120 120 1.5 2 2.5 3 3.5 4 x 103 50 100 150 200 250 300 350 400 450 500 Gap length [m]Total force per meter depth [N/m] Approximate solution EBC ABC Infinite Figure 7-46. Magnetic fo rce versus gap length Figure 7-46 show the total force on the clapper armature as a function of the gap width. All computed values of total forc e follow the approximate solution using the reluctance method. The computed force us ing ABC is the best answer among them, a force that exits between the force using EB C and the force using the infinite element method. Example 7-4-3: Force Betw een Two Parallel Wires Two long straight parallel wire s are separated by a distance d The wires carry the same amount of current in the same di rection as shown in Figure 7-47. Then an attractive force is crea ted between the two wires. PAGE 121 121 Figure 7-47. Two long parallel wires Using the AmpereÂ’s law, t he magnetic field at the 2nd wire caused by the current 1 I can be obtained as 01 12 I B d According to the Lorentz force equation, the magnetic force per unit length at the 2nd wire is stated as 012 22 I I F ld (7-7) where, l is the wire length. The same amount of the magnetic force is obtained at the 1st wire. When 1 md and 121 A II the magnetic force per unit length becomes 7210 N/m. The radius of the wire is 0.05 m, and the surrounding air area of 21.91.9 m is modeled. Open boundary condition is appl ied on the outer boundary using the truncation method, the asymptotic boundary c onditions, or infinite element method. Quadratic B-spline el ements were used. The total number of nodes was 10520. The contour plots of the computed magnetic vect or potential are shown in Figure 7-48. The magnetic force per unit length is computed as follows; 72.5310 N/m using the PAGE 122 122 truncation approach, 71.6310 N/m using the asymptotic boundary condition, and 72.4210 N/m using the decay function infinite element method with L =0.10. Using the asymptotic boundary condition, the computed forc e is closer to the analytical solution when 1 md. The other methods show similar resu lts in terms of the contour plot and the magnetic force. It is because the decay length is so short that the decay function infinite element method is as sa me as the truncation approach. Figure 7-48. Contour plots fo r magnetic vector potential. A) Homogenous essential boundary conditions, B) asymptotic boundary conditions and C) Decay function infinite element with L =0.1. Figure 7-49 shows the magnetic force per uni t length versus the distance between the two wires. The distance d was varied from 0.2 m to 1.6m Within a range from 0.2 m to 0.8 m, the decay factor L =0.15 was used for the decay func tion infinite element method and using this value the computed force using this method were close to the analytical solution. The computed forces using all open boundary techniques are similar when the distance is equal to 0.2 m, 0.4m, or 0.6m; howev er, the force comp utation using the asymptotic boundary condition is much closer to the analytic solution as the distance increases, that is when the source approac hes the boundary. As the magnetic source (currant carrying wire) is closer to the out boundary, the error in the magnetic force increases when the truncation approach or the infinite element method is used. PAGE 123 123 Therefore, the force comput ation using ABC is clearly better than the other two open boundary techniques. Figure 7-49. The magnetic force vers us the distance between two wires Coupled Magneto-Elastostatic Problems Two coupled magneto-elastostatic problems are examined to validate the implicit boundary method. The first example is a 2D example containing a clapper solenoid actuator and a cantilever beam. The cantilever beam is attached on the top surface of the clapper armature. The magne tic force from the actuat or deforms the cantilever beam. The second example is a 3D plunger sol enoid actuator workin g with one of three circular shaped plates. Even though the geomet ries of the plates are different and complicated, the density of t he structured mesh is the same for all three plates. PAGE 124 124 Example 7-5-1: 2D Clapper Solenoi d Actuator with Cantilever Beam A cantilever beam is attached on the top su rface of the armature. The cantilever beam is a beam fixed at one end. At the other end, this c antilever beam is attached to the armature. The attachment location varies as shown in Figure 7-50. The distances from the center axis to the attached location are 10mm, 20mm, and 28mm. The thickness of the beam is 5mm and the length of the beam is 80mm. The beam is made of aluminum of which material properties are 69 Gpa E 0.29 and 1r Figure 7-50. Planar clapper solenoid actuat or with a cantilever beam. A) d=10mm B) d=20mm and C) d=28mm This coupled analysis can involve nonlinear properties, nonlinearity that can occur from the large deformation of the beam, changing force according to the gap length, and the contact between the armature and the stator In this research, the nonlinearity is ignored during the analysis. The computed m agnetic force is used as loads for the cantilever beam. Four node bi linear elements are used for the analysis. The total number of elements is 1550. Fi gure 7-51 shows tip displacement versus ampere-turns. As NI increases, the displacement on the tip increases. When d=20mm or 28mm, the displacements on the tip is much larger than when d=10mm; howev er, there is very small difference in tip deflection between the two cases: d=20mm and d=28mm. PAGE 125 125 100 200 300 400 500 600 700 800 900 1000 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 NI [Ampereturns]Displacement on the tip [mm] d=10mm d=20mm d=28mm Figure 7-51. Displacement at the tip versus NI with varying attachment location Next, the distance between the symmetric ax is and the tip of the cantilever beam is fixed at d=20mm, and the thickness of the bar is varied as shown in Figure 7-52. The thickness value, equal to 3m m, 5mm or 7mm, was used. Figure 7-52. Planar clapper solenoid actuator with a cantilever beam. A) t=3mm, B) t=5mm and C) t=7mm PAGE 126 126 Figure 7-53 shows the tip displacement versus NI When t=3mm, the larger displacement can be observed. When NI increases, the difference by the thickness also increases. The tip deflection is very sensit ive to the thickness of the cantilever beam. 100 200 300 400 500 600 700 800 900 1000 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 NI [Ampereturns]Displacement on the tip [mm] t=3mm t=5mm t=7mm Figure 7-53. Displacement on the tip versus NI with varying beam thickness Example 7-5-2: 3D Plunger Sole noid Actuator with Structures The 3D model of a plunger solenoid act uator shown in Figure 7-28 A works with one of three plates, one that is attached to t he top surface of the plunger armature. The top views of each structure and dimensions are included in the Figure 7-54. The first model is a solid plate, the second is a pl ate with one hole, and the third one is a plate with two holes. The thicknesses of all the plates equal to 2mm. All of the three structures are attached to the top surface of the plunger arma ture as shown in Figure 755. The structures ar e made of aluminum. PAGE 127 127 Figure 7-54. Top views of structures attached to the plunger armature. A) Solid plate, B) plate with one hole, and C) plate with two holes When NI =800 and the fixed boundary conditions are applied on the su rfaces of the structure at R=70mm, the m agnetic force from the actu ator results in downward deflection of the plunger armature. The struct ures are deformed by the movement of the armature as shown in Figure 7-56. The maximum displacements for the structures are 77.2410 mm in the solid plate, 61.08610mm in the plate with one hole, and 78.53910mm in the plate including two holes. Figure 7-57 shows the maximum displacement versus NI (ampere-turns) for each plate. Among the structures, t he plate with one hole has the largest deformation. On the other hand, the solid plate has the highest stiffness among them. PAGE 128 128 A B C Figure 7-55. 3D plunger solenoi d actuators with A) Solid pl ate, B) plate with one hole, and C) plate with two holes PAGE 129 129 A B C Figure 7-56. Deformation due to magnetic force. A) Solid plate, B) plate with one hole, and C) plate with two holes PAGE 130 130 1000 1500 2000 2500 3000 2 4 6 8 10 12 14 16 x 106 NI (Ampereturns)Max. displacement (mm) Solid One hole Two holes Figure 7-57. Maximum displacement versus NI PAGE 131 131 CHAPTER 8 DESIGN AND ANALYSIS OF MAGN ETIC ACTUATORS USING IBFEM Design Criteria In this chapter, IBFEM is used as an act uator design tool. So lenoid actuators with clapper armature or plunger armature and coil actuators are designed under given design criteria. Among them, a lighter and more efficient design are needed for flapping wing micro air vehicle. The micro air vehicles are significantly small size aerial vehicle with the maximum dimension of about 15 cm. Fo r a micro air vehicle actuator, there are three given design criteria: coil winding area, iron weight and size. The given coil winding area wS is assumed to be 212 mm. The entire coil winding area cannot be assumed to be only the copper because the coil includes copper and other materials such as insulator. Thus, the relation between the coil winding area and the copper conductor is defined as the packing factor pF [18]. When the coil is wound tightly, the packing factor can be up to 75%. Thus, the packing factor is assumed to be 70%. The area of the coil wire is 325.0110cSmm when the number of the American Wire Gauge (AWG) is 40. Therefor e, the number of coil turns N is given as following equation 1680pw cFS N S (8-1) When the coil is wound to a cylindrical bobb in with the diameter of 4 mm, the total length of the coil becomes about 21.1 m. As the resistance per meter for the given AWG is 3.44 /m, the total resistance of the coil is approximately 72 In this research, it is assumed that the coil wire can allow current of 100mA to flow. The allowable current of the wire is based on plastic insulation. In case of the current source, LT3092, PAGE 132 132 manufactured by Linear technology, the ma ximum output current is 200mA. The input voltage range is from 1.2 V to 40 V so that LT3092 can be operated by using a small size battery. Thus, the current of 100mA can be obtained using LT3092. If N =1680 and the amount of current I can be controlled by a digital processor, NI (ampere-turns) can vary from 0 to 168. NI is also called magnetomotive forc e. In order to compare several designed actuators in the later se ction, the magnetomotive force, NI =30, is used. When the maximum current (100mA) flows in the coils, the coils can have maximum energy dissipation as 20.72 [] WIR W. Another design criterion is the iron weight of the act uator. The iron occupy large portion of the weight of the ac tuator; however, iron has high re lative permeability so that usage of the iron material can intensify a m agnetic force of the ac tuator. Therefore, there is a trade-off between the iron weight of the actuator and the usage of the iron. As the second design criterion, we assume that t he total iron weight of the iron portion can be up to 10g. In order to compute the iron weight of an actuator design, the mass density of iron is used as 3337874/7.87410/ kgmgmm. The third design criterion is the size of the actuator. In or der to install an actuator insi de the micro air vehicle, the actuator should fit inside of a cube with edge length of 2cm. In the beginning of a solenoid design, t he reluctance method is used so that a simple design is created under the design criter ia. Based on the initial design, geometry of the design is changed in order to determine a better design. Solenoid Actuators with Clapper Armature A solenoid actuator with a clapper armature is shown in Figure 8-1. The clapper armature can move in the y-direction when the coils carry current. The armature and the PAGE 133 133 stator are made of thin steel laminations with relative pe rmeability of 2000, a thin steel lamination that can reduce heating caused by eddy current. When the current flows through a coil, magnetic flux is created as shown in Figure 8-1. The magnetic flux follows a closed path. If magnit ude of the magnetic flux is known, a magnetic force can be estimated using the Maxw ell stress tensor method. Figure 8-1. Solenoid actuat or with clapper armature Supposing that the dimensions are 130.5llmm 253llmm 466.5llmm, 1 gmm and 1231wwwmm, then reluctances in the armature, the stator, and the gap can be obtained as follows 3 0123 1.59210A/Wb 1armature rlll R w (8-2) PAGE 134 134 3 0456 6.36410A/Wb 1stator rlll R w (8-3) 6 02 1.59110A/Wb 1gapg R w (8-4) Using the reluctance method, the magnetic flux can be stated as 76.28510Wbiarmaturestatorgap iNINI NI RRRR (8-5) As the magnetic flux is constant on the given path, the magnetic flux density and the magnetic field density can be obtained in each airgap as follows 46.28510T 1gapBNI w (8-6) 01 500A/mgapgapHBNI (8-7) Using MaxwellÂ’s stress tensor method, t he normal magnetic pressure can be obtained as 2 2 4 0416.28310N 2gapFHwNI (8-8) When 30NI, the magnetic force per unit lengt h becomes 0.565 N/m. The magnetic flux and the magnetic field density are 21.88610T and 41.510A/m. Using the same geometry of the actuator Figure 8-2 A shows contour plot of magnetic vector potential computed using IBFEM. Figure 8-2 B shows magnetic field density in the ydirection. The computed magnetic field densit y in the gap is approximately equal to the analytical value of 41.510A/m. For this analysis, the asymptotic boundary conditions are applied on the outer boundaries. The com puted force is 0.633 N/m. The computed force and the analytical solution are quite close. The difference between them results from a fringing flux, which the reluctanc e method ignores. As the iron domain is 240mm, the iron weight per unit depth is 0.315/gmm. PAGE 135 135 Figure 8-2. IBFEM results. A) Contour pl ot of magnetic vect or potential, and B) Magnetic field density Based on the initial design in shown in Figure 8-1, three different designs are created by changing the geometry of the armature and the st ator. For all the designs, the rectangular area of the 2D model including the air domain, is 21313mm. Figure 8-3 shows three designs with dimensions. The dimensions shown are all in mm Even though the geometry is changed for each actuator the mesh densities of the structured mesh are same for all designs. Nine node B-sp line elements were used for the analysis. The total number of nodes was 2586. The asymptotic boundary conditions were applied on the outer boundaries. T he magnetomotive force,30NI was applied. Figure 8-3. Clapper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3 PAGE 136 136 Figure 8-4 shows contour plots of magnetic vector potentials for all designs. Based on those computations, Table 8-1 shows com parison for three designs in terms of the magnetic force and the iron weight. Figure 8-4. Contour plots of magnetic vector potential for cl apper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3 PAGE 137 137 Table 8-1. Comparison for th ree clapper solenoid actuators ( NI =30) Design 1 Design 2 Design 3 Force per m 0.691 N/m 0.509 N/m 0.681 N/mm Iron weight per mm 0.315/gmm 0.283/gmm 0.320/gmm According to Table 8-1, the first actuator design produces the largest force among them. The second actuator is the lightest one. Solenoid Actuators with Plunger Armature A solenoid actuator with a plunger armature is shown in Figure 8-5. The shape of the plunger armature can be a brick, cyli nder, or conics. There are two gaps of g and gs that the flux can enter and leave the plunger armature. Useful magnetic force is produced only at g Figure 8-5. Solenoid actuat or with plunger armature PAGE 138 138 The armature and the stator are made of steel with the re lative permeability of 2000. When the dimensions are 17lmm, 22.5lmm 30.9lmm 44.5lmm 51.5lmm, 63.5lmm, 0.1gsmm, 1gmm and 12341wwwwmm all reluctances can be calculated as follows 3 034 2.14810A/Wb 1armature rll R w (8-9) 3 01256 5.76910A/Wb 1stator rllll R w (8-10) 4 07.95810A/Wb 1gsgs R w (8-11) 5 07.95810A/Wb 1gg R w (8-12) Using the reluctance method, the magnetic flux is 61.13210Wbiarmaturestatorggs iNINI NI RRRRR (8-13) Based on the magnetic flux, t he magnetic flux density and the magnetic field density in each airgap can be computed as follow 31.13210T 1gapBNI w (8-14) 01 900.8A/mgapgapHBNI (8-15) Using the magnetic field in the gap, the us eful magnetic force can be calculated as follow 2 2 -3 0211.0210N 2gapFHwNI (8-16) when 30NI, the force becomes 0.918 N. The m agnetic flux and the magnetic field density are 23.410T and 42.70210A/m. Figure 8-6 shows the computed results: contour plot of magnetic vect or potential and magnetic field density in the y-direction. PAGE 139 139 For the analysis, the same mesh density wa s used in the clapper solenoid actuator model. Nine node B-spline elements were us ed. Asymptotic boundary conditions were applied on the outer b oundaries. The computed total force on the plunger armature is 0.9130 N/m. As the tota l area of the iron is 238.8 mm, the iron weight per unit depth is 0.306/ gmm. Figure 8-6. IBFEM results. A) Contour pl ot of magnetic vect or potential, and B) Magnetic field in the y-direction Modifying the initial design in shown in Figure 8-5, three different designs are created as shown in Fi gure 8-7. The dimensi ons shown are all in mm For all designs, the area of the 2D model is as same as one of the clapper solenoi d actuator. The same mesh density is applied for all the designs. Nine node B-spline elements were used for the analysis. The total number of nodes wa s 2548. The asymptotic boundary conditions were applied on the outer boundarie s. The magnetomotive force, 30NI, was used. PAGE 140 140 Figure 8-7. Plunger solenoid ac tuators. A) Design 1, B) Design 2, and C) Design 3 Figure 8-8 shows contour plots of magnetic vector potentials for all designs. Based on those results, Table 8-2 shows comparis on for three designs in terms of the magnetic force and the iron weight. According to Table 8-2, the third plunger solenoid actuator model can produce the largest forc e among them. The second actuator model is the lightest one. PAGE 141 141 Figure 8-8. Contour plot s of magnetic vector potential for plunger solenoid actuators. A) Design 1, B) Design 2, and C) Design 3. Table 8-2. Comparison for thr ee plunger solenoid actuators ( NI =30) Design 1 Design 2 Design 3 Force per m 1.154 N/m 0.262 N/m 2.801 N/mm Iron weight per mm 0.306/gmm 0.294/gmm 0.310/gmm PAGE 142 142 Solenoid Actuators with a Combined Plunger & Clapper Armature A combined plunger & clapper armature is an armature that includes a clapper & plunger to increase the force. Three designs were studied. The sizes of those actuators are similar to one of the clapper solenoid act uator in the previous section. Magnetic force and iron weight for each ac tuator are characterized. Figure 8-9 shows three solenoid act uators with a mixture armature. The dimensions shown are all in mm For all designs, the area of the 2D model is as same as one of the clapper solenoid actuator. For the analysis, the same mesh density is applied for all the designs. Nine node B-spline elements were used for the analysis. The total number of nodes was 2674. The asymptot ic boundary conditions were applied on the outer boundaries. T he magnetomotive force NI is equal to 30. Figure 8-9. Three combined plunger & clapper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3 After the analysis using these models, contour plots of magnetic vector potentials are show in Figure 8-10. Based on those resu lts, Table 8-3 shows comparison for three designs in terms of the magnetic force and t he iron weight. According to Table 8-3, the third one can produce the largest force among them. The second actuat or is the lightest one. PAGE 143 143 Figure 8-10. Contour plots of magnetic vector potential for combined plunger & clapper solenoid actuators. A) Design 1, B) Design 2, and C) Design 3. Table 8-3. Comparison for three combi ned plunger & clapper so lenoid actuators ( NI =30) Design 1 Design 2 Design 3 Magnetic force per m 0.961 N/m 0.413 N/m 1.30 N/m Iron weight per mm 0.342/ gmm 0.339/ gmm 0.340/ gmm PAGE 144 144 Coil Actuators In this section, three coil actuators ar e studied under the given design criteria. The coil actuators include permanent magnets so that higher magnet ic flux can be produced without the increase of NI. Moreover, the usage of the permanent magnets can reduce heating on the actuator. The co il actuator has magnetic force in the conductive coil instead of one on the armature of the sol enoid actuator, a magnet ic force that is computed using the Lorentz force equation. Figure 8-11 shows three coil actuators. The dimensions shown are all in mm For all the designs, the rectangular area of t he 2D model including the air domain is 21818 mm The coil winding area is the same as in previous solenoid actuators. Two permanent magnets have different direction with the same remanent flux of 0.8 T. The first design is based on a typical voice coil actuator in a loudspeaker. The second is a design to remove the iron lami nates from the first design. The third design has a gap at the bottom portion of the iron laminate, whic h allows the coil to move freely in the downward direction. For the analysis, the same mesh density is applied for all designs. Four node bilinear elements were used. The total number of nodes was 12036. The asymptotic boundary conditions were applied on the outer boundaries. NI is equal to 30. Figure 8-11. Three coil act uators. A) Design 1, B) Design 2, and C) Design 3 PAGE 145 145 After the analysis using these models, the magnitudes of the magnetic flux densities are show in Figure 8-12. As the Lor entz force is proporti onal to the magnetic flux density, a stronger magnetic flux density near the moving coils can create a larger force. The first design has the strongest m agnetic flux near the moving coils. Based on those results, Table 8-3 shows comparison fo r three designs in terms of the Lorentz force on the moving coil and the iron weight of the actuator. According to Table 8-4, the first coil actuator can produce the largest force among them. The second design is the best design if the actuator weight is the most critical criterion. Figure 8-12. Magnitude of B fields for three coil actuators. A) Design 1, B) Design 2, and C) Design 3. PAGE 146 146 Table 8-4. Comparison for coil actuators ( NI =30) Design 1 Design 2 Design 3 Force per m 21.0 N/m 2.64 N/m 9.12 N/m Iron weight per mm 0.428/gmm 0 /gmm 0.394/gmm The Best Actuator among the Designed Actuators Using IBFEM, four types of actuators were examined su ch as the clapper solenoid actuator, the plunger soleno id actuator, the combined plunger & clapper solenoid actuator, and the coil actuator. Among them, the coil actuator is the best actuator considering the given design criteria. Among se veral designs for the coil actuators, the first design is used for flapping wings of the micro air vehicle as shown in Figure 8-13. The first design could generate the largest force under the given magnetomotive force so that the magnetic force can make the largest deformation of flapping wings among them. Figure 8-13. The best magnet ic actuator among several designed actuators Considering the 2D model, the best coil act uator is designed as 3D solid model as shown in Figure 8-14. Consi dering the symmetry, one fourth of the model is created using commercial software, Pro/engineering. Figure 8-15 shows the dimension of the PAGE 147 147 3D coil actuator and the directions of t he permanent magnets. T he dimensions shown are all in mm Figure 8-14. Sold model of the 3D coil actuator Figure 8-15. The top view of the 3D coil actuator PAGE 148 148 The moving coils carry current. The coil has a conductivity value of 610 S/m. In order to perform the analysis, eight node brick elements were used. The total number of the nodes was 9651. According to the design criterion, the ma gnetomotive force, NI can vary from 0 to 168. When NI is equal to 30, the magnetic flux density of the coil actuator is shown in Figure 8-16. Figure 8-16. The magnitude of the magnetic flux densit y of the coil actuator The magnetic flux density, B, vary from 1.17 to 2.33 T nearby the moving coils, so the computed force is -0.114 N. When NI is equal to 10, 30, 60, 90, 120 or 168, the computed Lorentz force is shown in Figure 8-17. Figure 8-17 shows the linear relation between the magnetomotive forc e and the Lorentz force. T he maximum force is 0.683 N. PAGE 149 149 0 20 40 60 80 100 120 140 160 180 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 NI [A.turns]Lorentz force [N] Figure 8-17. Lorentz force of the coil actuator versus NI Coupled Magneto-Elastostatic Pr oblems with a Fla pping Wing Model The best coil actuator is embedded in a mi cro air vehicle with flapping wings as shown in Figure 8-18. The coil ac tuator can fit inside of the fuselage. The moving coil of the coil actuator is atta ched to the structure. Figure 8-18. The coil act uator with flapping wings PAGE 150 150 Four structures are designed as shown in Fi gure 8-19 structures that are created as surface models. The thickness of the st ructure assumed to be equal to 0.45mm. Additionally, it is assumed t hat the structure is made of aluminum. All the structures have a similar top view; however, they have different front and side views. The dimensions shown are all in mm. As the length in chordwis e direction is 26 mm, two coil actuators can be placed along the chordwise direction because the length of the coil actuator is 12.2mm according to Figure 8-15 Figure 8-20 also shows the four surface models. Figure 8-19. Four structures wit h flapping wings. A) Design 1, B) Design 2, C) Design 3, and D) Design 4. PAGE 151 151 Figure 8-20. Four surface structures with fl apping wings. A) Design 1, B) Design 2, C) Design 3, and D) Design 4. In order to perform analysis of thin shell-like structures the surface model analysis, IBFEM [65] have been extended to use 3D s hell elements that are 3D elements with three degrees of freedom per node. These she ll elements were used for the analysis. The total number of nodes in model is 666. Considering the symmetr y of the geometry, half of the structure was modeled for the analysis. The structured mesh and the boundary conditions are shown in Figure 8-21. The magnetic force acts downward so that the wing produces upstr oke. When one coil actuator is used, the applied force PAGE 152 152 varies from 0 to 0.342 N at the edge. Usi ng two coil actuators, the magnetic force can be double. Figure 8-21. Structured mesh and boundar y conditions of the first design As the wing produces upstroke, the wing can also create down-stroke. Using the coil actuator, the same amount of Lorentz force can be produced in the opposite direction by changing the direction of the current. When NI is equal to 30, the wing up-strokes and down-strokes of the four designs are s hown in Figure 8-22. For the up-stroke, the maximum displacement on the tip is 32.63310 mm in the first design, 34.23810 mm in the second design, 23.76510mm in the third design, and 11.96810mm in the fourth design. Magnitudes of the tip deflections during the wing up-stroke are the same as ones during the wing down-stroke because t he magnetic force is proportional to the displacement. Thus, the tip deflection can be double including both strokes. Among the four designs studied here, the last desi gn produces the largest tip deflection. PAGE 153 153 Figure 8-22. Displacement in the z-directi on during the wing stroke. A) Design 1, B) Design 2, C) Design 3, and D) Design 4. PAGE 154 154 When NI varies from 1 to 168, the tip deflect ion as a function of the magnetomotive force is shown in Figure 8-23. The first thr ee designs are too stiff for our application. 0 20 40 60 80 100 120 140 160 180 0 0.2 0.4 0.6 0.8 1 1.2 1.4 NI [A.turns]Tip displacement [mm] Design 1 Design 2 Design 3 Design 4 Figure 8-23. Tip displacement versus NI for wing upstroke PAGE 155 155 CHAPTER 9 CONCLUSIONS AND FUTURE WORK Conclusions The main research contributions made in this thesis are to extend Implicit Boundary Finite Element Method to Perform magnetostatic analysis Perform coupled magneto-elastostatic analysis Model multi-material system Implicit Boundary Finite Element Method (IBFEM) for analysis using structured mesh has been demonstrated for 2Dand 3D-magnetostatic models and coupled magneto-elastostatic models. This approach directly uses the geometry imported from CAD systems without generating a conforming mesh. Structured meshes are easy to generate and the elements ar e regular and not distorted as in traditional finite element mesh. Furthermore, the inter nal elements are identical to each other and have the same stiffness matrix thus reducin g the computation required. In order to perform magnetos tatic analysis, IBFEM require s the capability for the multi-material analysis. Ea ch material has its own structured mesh, and a modified solution structure is applied at the interfac e elements. The multi-material analysis in magnetostatic problems allows magnetic force to be computed, a m agnetic force that includes the magnetic surface force dens ity and the Lorentz fo rce. The magnetic surface force density is a force on ferro-mat erial, and the Lorentz force is a force on a current carrying coil due to the magnetic flux. A fter computing these forces, the forces can then be used in a subsequent structural analysis to perform coupled magnetostaticelastostatic analysis. PAGE 156 156 IBFEM has been evaluated using several 2D and 3D magnetostatic examples. Although those examples were complicated in terms of geometry, IBFEM could solve the problems without any dist orted element. The results of those examples were compared to one using the reluctance method or to analytic solution when available. Based on the results of the magnetic flux density and the magneti c field density, the magnetic force was computed accurately Through 2D problems, IBFEM presented faster convergence using quadratic B-spline elements, the method that created very accurate solutions with relatively fewe r nodes. Moreover, IBFEM could solve large electromagnetic problems. In case of a 3D model with complex geom etry, the structured mesh could allow relatively fewer elem ents to be used for the analysis. Consequently, the number of equations is smaller than the traditional finite element method. Three open boundary techniques were im plemented in the IBFEM software and evaluated using several examples. The open boundary techniques studied here include the truncation method, the asymptotic boundary conditions, and the decay function infinite element method. Several examples were used for comparing these three methods including permanent magnet, sol enoid actuator, and two wire examples. Among the three methods studied, the asym ptotic boundary conditions produced the most robust results no matter where magnetic sources were located. The computed magnetic force could be s ubsequently used in the structural analysis. So, IBFEM could reveal relations hip that exists bet ween the magnetomotive force NI and the displacement of t he structure. This relation is an important factor in designing a magnetically actuated structure. Several examples in 2D or 3D were created using a variety of structur e geometry. Through the coupled magneto- PAGE 157 157 elastostatic analysis, the relationships bet ween magnetomotive force and displacement were studied for these examples. IBFEM was used as a magnetic actuator desi gn tool. In order to design an efficient magnetic actuator for a flappi ng-wing micro air vehicle, IBFEM was used to analyze several magnetic actuators created under the gi ven specifications. These specifications include the range of NI the size of a magnetic actuat or and the weight of a magnetic actuator. Several magnetic actuators were ex amined such as clapper solenoid actuator, plunger solenoid actuator, combined clap per & plunger solenoid actuator and coil actuator. Among them, the coil actuator created the largest magnetic force for the given specifications. This magnetic force was used in structural analysi s of flapping-wings modeled using 3D shell elements. Through t he coupled magneto-elastostatic analysis, the relationship between t he magnetomotive force of th e coil actuator and the tip displacement of the flappingwing structure was determined. Future Work This research focused on the steady static analysis so that the coupled problem is modeled as the weakly coupled magneto-elastostatic analysis. This assumption allows us to perform the magnetostatic and elasto static analysis sequentia lly. If the magnetic field changes due to the elastic deforma tion, the problem becomes nonlinear and a strongly coupled nonlinear analys is is needed. The magnetic field may change because the structural deformation may cause perm anent magnets or circuits attached to the structure to also move re lative to each other. If the structure undergoes large deformation then geometric nonlinearities must be included. In additi on, if the coupled problem is a dynamic problem, then the electr ical circuits, the magnetic circuits and the PAGE 158 158 structural analysis must be solved as a coupl ed problem. Another source of nonlinearity is due to contact between the armature and the stator. The open boundary techniques should be further studied and extended for 3D magnetostatic problems. The open boundary technique s were studied in this thesis only for two dimensional problems. Asymptotic Boundary Condition (ABC) was found to be the most effective; however, th is approach requires the user to specify the origin of the magnetic source. The accuracy of the results depends on careful choice of this origin. Further, research is needed to explore more effective ways of applying open boundary condition. The magnetic force computation should be extended to nonlinear materials and permanent magnets. In this thesis, the m agnetic force computat ion was implemented and studied only for linear ferrro-material and cu rrent carrying coils. In order to analyze magnetic actuators that involve moving permanent magnets, we need to implement force computation for permanent magnets. In order to do dynamic analysis for flapping wings operated by magnetic actuators, nonlinear characteristics in geometry, forc e and kinematics must be considered. The wing motion creates large deformation so that a moment can be a function of deformation. When an armature of magnetic actuator moves toward stator, the magnetic force can vary according to the location of the armature. In addition, the contact problem occurs between the armature and the stator. PAGE 159 159 LIST OF REFERENCES [1-67] [1] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, "Meshless methods: An overview and recent developments," Computer Methods in Applied Mechanics and Engineering vol. 139, pp. 3-47, 1996. [2] K. M. Lee, Q. Li, and H. Sun, "Effects of numerical formulation on magnetic field computation using meshless methods," IEEE Transactions on Magnetics vol. 42, pp. 2164-2171, 2006. [3] Q. Li and K. M. Lee, "An adaptive meshless method for magnetic field computation," IEEE Transactions on Magnetics vol. 42, pp. 1996-2003, 2006. [4] A. Bruyere, L. Illoul, F. Chines ta, and S. Clenet, "Comparison between NEM and FEM in 2-D magnetostatics using an error estimator," IEEE Transactions on Magnetics vol. 44, pp. 1342-1345, 2008. [5] X. Liu, Y. M. Deng, Z. W. Zeng, L. Udpa, and S. S. Udpa, "Model-Based Inversion Technique Using ElementFree Galerkin Method and State Space Search," IEEE Transactions on Magnetics vol. 45, pp. 1486-1489, 2009. [6] O. Bottauscio, M. Chiampi, and A. Manz in, "Computation of higher order spatial derivatives in the multiscale expansi on of electromagnetic-field problems," IEEE Transactions on Magnetics vol. 44, pp. 1194-1197, 2008. [7] K. Muramatsu, Y. Yokoyama, and N. Takahashi, "3-D magnetic field analysis using nonconforming mesh with edge elements," IEEE Transactions on Magnetics vol. 38, pp. 433-436, 2002. [8] T. Belytschko, C. Parimi, N. Moes N. Sukumar, and S. Usui, "Structured extended finite element me thods for solids defined by implicit surfaces," International Journal for Nume rical Methods in Engineering vol. 56, pp. 609-635, 2003. [9] N. Moes, E. Bechet, and M. Tourbier "Imposing Dirichlet boundary conditions in the extended finite element method," International Journal for Numerical Methods in Engineering vol. 67, pp. 1641-1669, 2006. [10] Y. Abdelaziz and A. Hamouine, "A su rvey of the extended finite element," Computers & Structures vol. 86, pp. 1141-1151, 2008. [11] P. M. A. Areias and T. Belytschko, "A nalysis of three-dimensional crack initiation and propagation using the extended finite element method," International Journal for Numerical Methods in Engineering vol. 63, pp. 760-788, 2005. [12] T. Belytschko and T. Black, "Elastic crack growth in finite elements with minimal remeshing," International Journal for Nume rical Methods in Engineering vol. 45, pp. 601-620, 1999. [13] A. V. Kumar, R. Bu rla, S. Padmanabhan, and L. X. Gu, "Finite element analysis using nonconforming mesh," Journal of Computing and Information Science in Engineering vol. 8, 2008. [14] A. V. Kumar, S. P admanabhan, and R. Burla, "Implicit boundary method for finite element analysis using non-conforming mesh or grid," International Journal for Numerical Methods in Engineering vol. 74, pp. 1421-1447, 2008. [15] R. K. Burla and A. V. Kumar, "Implicit boundary method for analysis using uniform B-spline basis and structured grid," International Journal for Numerical Methods in Engineering, vol. 76, pp. 1993-2028, 2008. PAGE 160 160 [16] R. K. Burla, A. V. Kumar, and B. V. Sankar, "Imp licit boundary method for determination of effective properties of composite microstructures," International Journal of Solids and Structures vol. 46, pp. 2514-2526, 2009. [17] S. R. H. Hoole, Computer-aided analysis and design of electromagnetic devices New York: Elsevier, 1989. [18] J. R. Brauer, Magnetic Actuators and Sensors : John Wiley & Sons, 2006. [19] K. Vijayakumar, R. Karthikeyan, S. Paramasivam, R. Arumugam, and K. N. Srinivas, "Switched Reluctance Motor Modeling, Design, Simulation, and Analysis: A Comprehensive Review," IEEE Transactions on Magnetics vol. 44, pp. 4605-4617, 2008. [20] V. L. Rvachev and T. I. Sheiko, "R-F unctions in Boundary Value Problems in Mechanics," Applied Mechanics Reviews vol. 48, pp. 151-188, 1995. [21] V. Shapiro and I. Tsukanov, "Meshfr ee simulation of deforming domains," Computer-Aided Design vol. 31, pp. 459-471, 1999. [22] K. Hollig, C. Appric h, and A. Streit, "Introducti on to the Web-method and its applications," Advances in Computational Mathematics vol. 23, pp. 215-237, 2005. [23] G. Apaydin, S. Seker, and N. Ar i, "Weighted extend ed b-splines for onedimensional electrom agnetic problems," Applied Mathemati cs and Computation vol. 190, pp. 1125-1135, 2007. [24] G. Apaydin, S. Seker, and N. Ari, "Application of w eb-spline method in electromagnetics," Aeu-International Journal of Electronics and Communications vol. 62, pp. 163-173, 2008. [25] J. L. Coulomb and G. Meunier, "Finit e-Element Implementati on of Virtual Work Principle for Magnetic or Electr ic Force and Torque Computation," IEEE Transactions on Magnetics vol. 20, pp. 1894-1896, 1984. [26] Z. Ren and A. Bossavit, "A new approach to eddy current problems in deformable conductors and some numer ical evidence about its validity," International Journal for Applie d Electromagnetics in Materials vol. 3, pp. 39-46, 1992. [27] Z. Ren and A. Razek, "Local Force Computation in Defo rmable-Bodies Using Edge Elements," IEEE Transactions on Magnetics vol. 28, pp. 1212-1215, 1992. [28] Z. Ren, "Comparison of Different Forc e Calculation Methods in 3d Finite-Element Modeling," IEEE Transactions on Magnetics vol. 30, pp. 3471-3474, 1994. [29] J. Bastos and N. Sadowski, Electromagnetic modeling by finite element methods : CRC Press, 2003. [30] R. S. Grandia, V. A. Ga lindo, A. U. Galve, and R. V. Fos, "General formulation for magnetic forces in linear ma terials and permanent magnets," IEEE Transactions on Magnetics vol. 44, pp. 2134-2140, 2008. [31] J. Simkin and C. W. Trowbridge, "Use of the Total Scalar Potential in the Numerical-Solution of Field Problems in Electromagnetics," International Journal for Numerical Methods in Engineering vol. 14, pp. 423-440, 1979. [32] J. L. Coulomb, "Finite-Element 3 Dimensional Magnetic-Field Computation," IEEE Transactions on Magnetics vol. 17, pp. 3241-3246, 1981. PAGE 161 161 [33] O. Biro and K. Preis, "O n the Use of the Magnetic Vect or Potential in the FiniteElement Analysis of 3-Dim ensional Eddy Currents," IEEE Transactions on Magnetics vol. 25, pp. 3145-3159, 1989. [34] K. Preis, I. Bardi, O. Biro, C. Magele, W. Renhart, K. R. Richter, and G. Vrisk, "Numerical-Analysis of 3d Magnetostatic Fields," IEEE Transactions on Magnetics vol. 27, pp. 3798-3803, 1991. [35] R. Albanese and G. Rubinacci, "Magnetostati c Field Computations in Terms of 2Component Vector Potentials," International Journal for Numerical Methods in Engineering vol. 29, pp. 515-532, 1990. [36] R. C. Mesquita and J. P. A. Bastos, "An Incomp lete Gauge Formulation for 3d Nodal Finite-Element Magnetostatics," IEEE Transactions on Magnetics vol. 28, pp. 1044-1047, 1992. [37] N. A. Demerdash and R. Wang, "Theoretical and Numerical Difficulties in 3-D Vector Potential Methods in Finite-E lement Magnetostatic Computations," IEEE Transactions on Magnetics vol. 26, pp. 1656-1658, 1990. [38] A. S. Semenov, H. Kessler, A. Liskow sky, and H. Balke, "On a vector potential formulation for 3D electromechani cal finite element analysis," Communications in Numerical Methods in Engineering vol. 22, pp. 357-375, 2006. [39] J. Jin, The Finite Element Method in Electromagnetics 2 ed: Wiley-IEEE Press, 2002. [40] Q. Chen and A. Konrad, "A review of finite element open boundary techniques for static and quasi-static elec tromagnetic field problems," IEEE Transactions on Magnetics vol. 33, pp. 663-676, 1997. [41] G. Ventura, "On the elimination of quadrature s ubcells for discontinuous functions in the eXtended Finite-Element Method," International Journal for Numerical Methods in Engineering vol. 66, pp. 761-795, 2006. [42] S. Natarajan, S. Bordas, and D. R. Mahapatra, "Numerical integration over arbitrary polygonal domains based on Sc hwarz-Christoffel conformal mapping," International Journal for Nume rical Methods in Engineering vol. 80, pp. 103-134, 2009. [43] D. A. Lowther, C. B. Rajanathan, and P. P. Silvester, "Finite-Element Technique for Solving 2-D Open Boundary Problems," IEEE Transactions on Magnetics vol. 14, pp. 467-469, 1978. [44] J. R. Brauer, "Open Boundary Finite-Elements for Axisymmetric Magnetic and Skin Effect Problems," Journal of Applied Physics vol. 53, pp. 8366-8368, 1982. [45] C. Antunes, "Approxim ate Ballooning Techniques," IEEE Transactions on Magnetics vol. 19, pp. 2555-2557, 1983. [46] C. Antunes, E. M. Free man, D. Lowther, and P. Silves ter, "A Static Ballooning Technique for 2-D Open Boundary-Problems," Journal of Applied Physics vol. 53, pp. 8360-8362, 1982. [47] Y. S. Sun, G. S. Zhang, and W. Chen, "A Data-Base Method for 3d-Open Boundary Field Computation," IEEE Transactions on Magnetics vol. 26, pp. 807810, 1990. [48] W. Chen, G. S. Zhang, Y. S. Sun, M. J. Chen, and Y. Zhao, "The Application of the Data-Base Method for 3-D Open Boundary Field Computation," IEEE Transactions on Magnetics vol. 28, pp. 1686-1689, 1992. PAGE 162 162 [49] J. R. Brauer, S. M. Schaefer, J. F. Lee, and R. Mittra "Asymptotic BoundaryCondition for 3-Dimensional M agnetostatic Finite-Elements," IEEE Transactions on Magnetics vol. 27, pp. 5013-5015, 1991. [50] Q. S. Chen, A. Konrad, and S. Bar onijan, "Asymptotic Boundary-Conditions for Axisymmetrical Finite-Element Electrostatic Analysis," IEEE Transactions on Magnetics vol. 30, pp. 4335-4337, 1994. [51] Q. Chen, A. Konrad, and P. P. Biringer, "Computati on of 3-Dimensional Unbounded Eddy-Current Problems Using Asymptotic Boundar y-Conditions," IEEE Transactions on Magnetics vol. 31, pp. 1348-1351, 1995. [52] S. Gratkowski, L. Pi chon, and H. Gajan, "Asymptotic boundary conditions for open boundaries of axisymmetric magnet ostatic finite-element models," IEEE Transactions on Magnetics vol. 38, pp. 469-472, 2002. [53] R. L. Ungless, "An infinite finite el ement," vol. Master of Science. Canada: The University of British Columbia, 1973. [54] P. Bettess, "I nfinite Elements," International Journal for Numerical Methods in Engineering vol. 11, pp. 53-64, 1977. [55] P. Bettess, "More on Infinite Elements," International Journal for Numerical Methods in Engineering vol. 15, pp. 1613-1626, 1980. [56] P. Bettess, Infinite elements Sunderland, UK: Penshaw Press, 1993. [57] B. M. A. Rahman and J. B. Davies "Finite-Element Analysis of Optical and Microwave Waveguide Problems," IEEE Transactions on Microwave Theory and Techniques vol. 32, pp. 20-28, 1984. [58] M. J. McDougall and J. P. Webb, "Infinite Elements for the Analysis of Open Dielectric Wave-Guides," IEEE Transactions on Microwave Theory and Techniques vol. 37, pp. 1724-1731, 1989. [59] M. S. Towers, A. McCowen, and J. A. R. Macnab, "Electromagnetic Scattering from an Arbitrary, Inhomogeneo us 2-D Object a Finite and Infinite Element Solution," IEEE Transactions on Antennas and Propagation vol. 41, pp. 770-777, 1993. [60] G. Beer and J. L. Meek "Infinite Domain Elements," International Journal for Numerical Methods in Engineering vol. 17, pp. 43-52, 1981. [61] O. C. Zienkiewicz, C. Emson, and P. Bettess, "A Novel Boundary Infinite Element," International Journal for Numerical Methods in Engineering vol. 19, pp. 393-404, 1983. [62] T. T. Abdel-Fattah, H. A. Hodhod, and A. Y. Akl, "A novel formulation of infinite elements for static analysis," Computers & Structures vol. 77, pp. 371-379, 2000. [63] M. V. K. Chari and S. J. Salon, Numerical methods in electromagnetism San Diego: Academic press, 2000. [64] R. Arumugam, D. A. Lo wther, R. Krishnan, and J. F. Lindsay, "Magnetic-Field Analysis of a Switched Reluctance Motor Using a Two-Dimensional FiniteElement Model," IEEE Transactions on Magnetics vol. 21, pp. 1883-1885, 1985. [65] P. D. S. Periyasamy, "Finite element anal ysis of shell like structures using implicit boundary method," in Mechanical and Aerospace engineering vol. Master of Science. Gainesville: Universi ty of Florida, 2009, pp. 86. [66] S. Zhang and A. V. Kumar, "Magnetic Field Computation us ing Implicit boundary Finite Element Method," IEEE Transactions on Magnetics in press, 2010. PAGE 163 163 [67] S. Zhang and A. V. Kumar, "3D Magnet ostatic Analysis using Implicit Boundary Finite Element Method," Computers & Structures submitted 2010. PAGE 164 164 BIOGRAPHICAL SKETCH Sung-Uk Zhang was born in Seoul, Sout h Korea. He graduated with a bachelorÂ’s degree in electrical engineering from the S ogang University, South Korea in February 2003. In December 2006, he achieved a ma sterÂ’s degree in bi omedical engineering from the University of Flor ida. Then, he enrolled in the doctoral problem in mechanical and aerospace engineering at the University of Florida in January 2007. His areas of interest include digital signal processing, electrical impedance tomography, inverse problem, finite element method, and computer-aided geometric design. |