BELA

 

Finite Element Electrostatics Solver

 

 

 

Version 1.0

User’s Manual

 

 

 

March 29, 2003

 

David Meeker

dmeeker@ieee.org

http://femm.berlios.de/bela.htm

 

©2003

Acknowledgements

 

Many thanks to those people who offered feedback on the beta version of Bela:

 

·        Mark Smith
 

 

 

1. Introduction

BELA is a suite of programs for solving linear electrostatic problems. The programs currently address problems on two-dimensional planar and axisymmetric domains. BELA is divided into three parts: 

 

  1. Preprocessor (beladraw.exe). This is a CAD-like program for laying out the geometry of the problem to be solved and defining material properties and boundary conditions. Autocad DXF files can be imported to facilitate the analysis of existing geometries.
  2. Solver (belasolv.exe). The solver takes a set of data files that describe problem and solves the relevant Maxwell’s equations to obtain values for the electric field through the solution domain.
  3. Postprocessor (belaview.exe). This is a graphical program that displays the resulting fields in the form of contour and density plots. The program also allows the user to inspect the field at arbitrary points, as well as evaluate a number of different integrals and plot various quantities of interest along user-defined contours. 

 

Two additional programs are also called to perform specialized tasks. These are:

 

 

The Lua scripting language (see http://www.lua.org) is also integrated into the pre- and post-processors. Lua allow “batch” runs to be performed without user interaction. In addition, all edit boxes in the user interface are parsed by Lua, allowing equations or mathematical expressions to be entered into any edit box in lieu of a numerical value. In any edit box in BELA, a selected piece of text can be evaluated by Lua via a selection on the right mouse button click menu.

 

The purpose of this document is to give a brief explanation of the kind of problems solved by BELA and to provide a fairly detailed documentation of the programs’ use.


 

2. Overview

 

The goal of this section is to give the user a brief description of the problems that BELA solves. This information is not really crucial if you are not particularly interested in the approach that BELA takes to formulating the problems. You can skip most of Overview, but take a look at Section 2.2. This section contains some important pointers about assigning enough boundary conditions to get a solvable problem.

 

Some familiarity with magnetism and Maxwell’s equations is assumed, since a review of this material is beyond the scope of this manual. However, the author has found several references that have proved useful in understanding the derivation and solution of Maxwell’s equations in various situations. A very good introductory-level text is Plonus’s Applied electromagnetics [1]. A good intermediate-level review of Maxwell’s equations, as well as a useful analogy of electromagnetism to similar problems in other disciplines is contained in Hoole’s Computer-aided analysis and design of electromagnetic devices [2]. For an advanced treatment, the reader has no recourse but to refer to Jackson’s Classical electrodynamics [3].

2.1. Relevant Maxwell’s Equations

 

For the electrostatics problems addressed by BELA, only a subset of Maxwell’s equations are required.  These problems describe the behavior of electric field intensity, E, and electric flux density (alternatively electric displacement), D.  There are two conditions that these quantities must obey.  The first condition is the integral form of Gauss’ Law, which says that the flux out of any closed volume is equal to the charge contained within the volume:

 

                                                         (1)

 

where r represents charge density.  The second is the integral form of Ampere’s loop law:

 

                                                        (2)

 

Displacement and field intensity are also related to one another via the constitutive relationship:

 

                                                           (3)

 

where e is the electrical permittivity.  Although some electrostatics problems might have a nonlinear constitutive relationship between D and E, BELA only considers linear problems.

 

To simplify the computation of fields which satisfy these conditions, BELA employs the electric scalar potential, V, defined by its relation to E as:

 

                                                        (4)

 

Because of the vector identity  for any scalar y, Ampere’s loop law is automatically satisfied.  Substituting into Gauss’ Law and applying the constitutive relationship yields the second-order partial differential equation:

 

                                                     (5)

 

which applies over regions of homogeneous e. This is the equation that BELA solves—the program solves for voltage V over a user-defined domain with user-defined sources and boundary conditions.

2.2. Boundary Conditions

 

Some discussion of boundary conditions is necessary so that the user will be sure to define an adequate number of boundary conditions to guarantee a unique solution. Boundary conditions for BELA come in three flavors: 

 

  1. Dirichlet. In this type of boundary condition, the value of  is explicitly defined on the boundary, e.g. . The most common use of Dirichlet-type boundary conditions is to define  along a boundary to keep flux from crossing the boundary.
  2. Neumann. This boundary condition specifies the normal derivative of  along the boundary. Usually,  is defined along a boundary to force flux to pass the boundary at exactly a  angle to the boundary. This sort of boundary condition is consistent with an interface with a very highly permeable metal.
  3. Robin. The Robin boundary condition is sort of a mix between Dirichlet and Neumann, prescribing a relationship between the value of  and its normal derivative at the boundary. An example of this boundary condition is:

 

                                            (6)

 

This boundary condition is most often used by BELA as a computationally cheap way of emulating unbounded solution domain. 

 

If no boundary conditions are explicitly defined, each boundary defaults to a  Neumann boundary condition. However, a non-derivative boundary condition must be defined somewhere so that the problem has a unique solution.

2.3. Finite Element Analysis

 

Although the differential equations that describe V appear relatively compact, it is very difficult to get closed-form solutions for all but the simplest geometries. That’s where finite element analysis comes in. The idea of finite elements is to break the problem down into a large number regions, each with a simple geometry (e.g. triangles). For example, Figure 1 shows a map of the Massachusetts broken down into triangles.

Figure 1. Triangulation of Massachusetts

Over these simple regions, the “true” solution for V is approximated by a very simple function. If enough small regions are used, the approximate V closely matches the exact V.

 

The advantage of breaking the domain down into a number of small elements is that the problem becomes transformed from a small but difficult to solve problem into a big but relatively easy to solve problem. Specifically, triangulation of the problem results in a linear algebra problem with perhaps tens of thousands of unknowns. However, techniques exist that allow the computer to solve for all the unknowns in only seconds.

 

BELA uses triangular elements. Over each element, the solution is approximated by a linear interpolation of the values of V at the three vertices of the triangle. The linear algebra problem is formed by choosing V on the basis of minimizing the total energy of the problem.


 

3. Preprocessor

 

The preprocessor is used for drawing the problems geometry, defining materials, and defining boundary conditions.  Drawing a valid geometry usually consists of four (though not necessarily sequential) tasks: 

 

·        Drawing the endpoints of the lines and arc segments that make up a drawing.

·        Connecting the endpoints with either line segments or arc segments

·        Adding “Block Label” markers into each section of the model to define material properties and mesh sizing for each section.

·        Specifying boundary conditions on the outer edges of the geometry. 

 

This section will describe exactly how one goes about performing these tasks and creating a problem that can be solved.

3.1. Preprocessor Drawing Modes

 

The key to using the preprocessor is that the preprocessor is always in one of five modes: the Point mode, the Segment mode, Arc Segment mode, the Block mode, or the Group mode. The first four of these modes correspond to the four types of entities that define the problems geometry: nodes that define all corners in the solution geometry, line segments and arc segments that connect the nodes to form boundaries and interfaces, and block labels that denote what material properties and mesh size are associated with each solution region. When the preprocessor is in a one of the first four drawing modes, editing operations take place only upon the selected type of entity. The fifth mode, the group mode, is meant to glue different objects together into parts so that entire parts can be manipulated more easily.

 

One can switch between drawing modes by clicking the appropriate button on the Drawing Mode potion of the toolbar. This section of the toolbar is pictured in Figure 2. 

Figure 2. Drawing Mode toolbar buttons.

The buttons correspond to Point, Line Segment, Arc Segment, Block Label, and Group modes respectively. The default drawing mode when the program begins is the Point mode.

3.2. Keyboard and Mouse Commands

 Although most of the tasks that need to be performed are available via the toolbar, some important functions are invoked only through the use of “hot” keys. A summary of these keys and their associated functions is contained in Table 1.


 

Point Mode Keys

Key

Function

Space

Edit the properties of selected point(s)

Tab

Display dialog for the numerical entry of coordinates for a new point

Escape

Unselect all points

Delete

Delete selected points

 

Line/Arc Segment Mode Keys

Key

Function

Space

Edit the properties of selected segment(s)

Escape

Unselect all segments and line starting points

Delete

Delete selected segment(s)

 

Block Label Mode Keys

Key

Function

Space

Edit the properties of selected block labels(s)

Tab

Display dialog for the numerical entry of coordinates for a new label

Escape

Unselect all block labels

Delete

Delete selected block label(s)

F3

Scale mesh size in all block labels by ½

F4

Scale mesh size in all block by 2

 

Group Mode Keys

Key

Function

Space

Edit the properties of selected objects

Escape

Unselect all

Delete

Delete selected objects(s)

 

View Manipulation Keys

Key

Function

Left Arrow

Pan left

Right Arrow

Pan right

Up Arrow

Pan up

Down Arrow

Pan down

Page Up

Zoom in

Page Down

Zoom out

Home”

Zoom “natural

 

Table 1: BELADRAW  hot keys

 

Likewise, specific functions are associated with mouse button input. The user employs the mouse to create new object, select obects that have already been created, and inquire about object properties. Table 2 is a summary of the mouse button click actions.


 

Point Mode

Action

Function

Left Button Click

Create a new point at the current mouse pointer location

Right Button Click

Select the nearest point

Right Button DblClick

Display coordinates of the nearest point

 

Line/Arc Segment Mode

Action

Function

Left Button Click

Select a start/end point for a new segment

Right Button Click

Select the nearest line/arc segment

Right Button DblClick

Display length of the nearest arc/line segment

 

Block Label Mode

Action

Function

Left Button Click

Create a new block label at the current mouse pointer location

Right Button Click

Select the nearest block label

Right Button DblClick

Display coordinates of the nearest block label

 

Group Mode

Action

Function

Right Button Click

Select the group associated with the nearest object

 

Table 2. BELADRAW Mouse button actions

3.3. View Manipulation

Generally, the user needs to size or move the view of the problem geometry displayed on the screen. Most of the view manipulation commands are available via buttons on the preprocessor toolbar. The functionality of can generally also be accessed via the ‘View Manipulation Keys’ listed in Table 1. The View Manipulation toolbar buttons are pictured in Figure 3. 

 

Figure 3. View Manipulation toolbar buttons.

The meaning of the View Manipulation toolbar buttons are:  The arrows on the toolbar correspond to moving the view in the direction of the arrow approximately  of the current screen width. The “blank page” button scales the screen to the smallest possible view that displays the entire problem geometry. The “+” and “-” buttons zoom the current view in and out, respectively. The “page with magnifying glass” button allows the view to be zoomed in on a user-specified part of the screen. To use this tool, first push the toolbar button. Then, move the mouse pointer to one of the desired corners of the “new” view. Press and hold the left mouse button. Drag the mouse pointer to the opposite diagonal corner of the desired “new” view. Last, release the left mouse button. The view will zoom in to a window that best fits the user’s desired window.

 

Some infrequently used view commands are also available, but only as options off of the Zoom selection of the main menu. This menu contains all of the manipulations available from the toolbar buttons, plus the options Keyboard, Status Bar, and Toolbar.

 

The Keyboard selection allows the user to zoom in to a window in which the window’s corners are explicitly specified by the user via keyboard entry of the corners’ coordinates. When this selection is chosen, a dialog pops up prompting for the locations of the window corners. Enter the desired window coordinates and hit “OK”. The view will then zoom to the smallest possible window that bounds the desired window corners. Typically, this view manipulation is only done as a new drawing is begun, to initially size the view window to convenient boundaries.

 

The Status Bar selection can be used to hide or show the one-line status bar at the bottom of the Beladraw window. Generally, it is desirable for the toolbar to be displayed, since the current location of the mouse pointer is displayed on the status line.

 

The Toolbar selection can be used to hide or show the toolbar buttons. The toolbar is not fundamentally necessary to running Beladraw, because any selection on the toolbar is also available via selections off of the main menu. If more space on the screen is desired, this option can be chosen to hide the toolbar. Selecting it a second time will show the toolbar again. It may be useful to note that the toolbar can be undocked from the main screen and made to “float” at a user-defined location on screen. This is done by pushing the left mouse button down on an area of the toolbar that is not actually a button, and then dragging the toolbar to its desired location. The toolbar can be docked again by moving it back to its original position.

3.4. Grid Manipulation

To aid in drawing your geometry, a useful tool is the Grid. When the grid is on, a grid of light blue pixels will be displayed on the screen. The spacing between grid points can be specified by the user, and the mouse pointer can be made to “snap” to the closest grid point.

 

The easiest way to manipulate the grid is through the used of the Grid Manipulation toolbar buttons. These buttons are pictured in Figure 4. 

 

Figure 4. Grid Manipulation toolbar buttons.

;  The left-most button in Figure 0 shows and hides the grid. The default is that the button is pushed in, showing the current grid. The second button, with an icon of an arrow pointing to a grid point, is the “snap to grid” button. When this button is pushed in, the location of the mouse pointer is rounded to the nearest grid point location. By default, the “snap to grid” button is not pressed. The right-most button brings up the Grid Properties dialog. This dialog is shown in Figure 5. 

 

Figure 5. Grid Properties dialog.

The Grid Properties dialog has an edit box for the user to enter the desired grid sizing. When the box appears, the number in this edit box is the current grid size. The edit box also contains a drop list that allows the user to select between Cartesian and Polar coordinates. If Cartesian is selected, points are specified by their (x,y) coordinates for a planar problem, or by their (r,z) coordinates for an axisymmetric problem. If Polar is selected, points are specified by an angle and a radial distance from the origin. The default is Cartesian coordinates.

3.5. Edit

Several useful tasks can be performed via the Edit menu off of the main menu.

 

Perhaps the most frequently used is the Undo command. Choosing this selection undoes the last addition or deletion that the user has made to the model’s geometry.

 

For selecting many objects quickly, the Select Group command is useful. This command allows the user to select objects of the current type located in an arbitrary rectangular box. When this command is selected, move the mouse pointer to one corner of the region that is to be selected. Press and hold the left mouse button. Then, drag the mouse pointer to the opposite diagonal corner of the region. A red box will appear, outlining the region to be selected. When the desired region has been specified, release the left mouse button. All objects of the current type completely contained within the box will become selected.

 

Any objects that are currently selected can be moved, copied, or pasted. To move or copy selected objects, simply choose the corresponding selection off of the main menu’s Edit menu. A dialog will appear prompting for an amount of displacement or rotation.

3.6. Problem Definition

The definition of problem type is specified by choosing the Problem selection off of the main menu. Selecting this option brings up the Problem Definition dialog, shown in Figure 6. 

Figure 6. Problem Definition dialog.

The first selection is the Problem Type drop list. This drop box allows the user to choose from a 2-D planar problem (the Planar selection), or an axisymmetric problem (the Axisymmetric selection).

 

Next is the Length Units drop list. This box identifies what unit is associated with the dimensions prescribed in the model’s geometry. Currently, the program supports inches, millimeters, centimeters, meters, mils, and meters.

 

The first edit box is the Depth specification. If a Planar problem is selected, this edit box becomes enabled. This value is the length of the geometry in the “into the page” direction. This value is used for scaling integral results in the post processor (e.g. force, inductance, etc.) to the appropriate length. The units of the Depth selection are the same as the selected length units.

 

The second edit box is the Solver Precision edit box. The number in this edit box specifies the stopping criteria for the linear solver. The linear algebra problem could be represented by:

 

Mx=b                                                              (7)

 

where  is a square matrix,  is a vector, and  is a vector of unknowns to be determined. The solver precision value determines the maximum allowable value for . The default value is .

 

Lastly, there is an optional Comment edit box. The user can enter in a few lines of text that give a brief description of the problem that is being solved. This is useful if the user is running several small variations on a given geometry. The comment can then be used to identify the relevant features for a particular geometry.

3.7. Definition of Properties

To make a solvable problem definition, the user must identify boundary conditions, block materials properties, and so on. The different types of properties defined for a given problem are defined via the Properties selection off of the main menu.

When the Properties selection is chosen, a drop menu appears that has selections for Materials, Boundary, Point, and Conductors. When any one of these selections is chosen, the dialog pictured in Figure 7 appears. 

 

Figure 7. Property Definition dialog box

This dialog is the manager for a particular type of properties. All currently defined properties are displayed in the Property Name drop list at the top of the dialog. At the beginning of a new model definition, the box will be blank, since no properties have yet been defined. Pushing the Add Property button allows the user to define a new property type. The Delete Property button removes the definition of the property currently in view in the Property Name box. The Modify Property button allows the user to view and edit the property currently selected in the Property Name box. Specifics for defining Point, Segment, and Block properties are addressed in the following subsections.

3.7.1. Point Properties

If a new point property is added or an existing point property modified, the Nodal Property dialog box appears. This dialog box is pictured in Figure 8.

 

Figure 8. Nodal Property dialog.

The first selection is the Name edit box. The default name is “New Point Property,” but this name should be changed to something that describes the property that you are defining.

 

Next are edit boxes for defining the voltage at a given point, or prescribing a point charge density at a given point. The type of point property is chosen via the radio buttons, and the value is entered in the enabled edit box.

3.7.2. Boundary Properties

 

The Boundary Property dialog box is used to specify the properties of line segments or arc segments that are to be boundaries of the solution domain. When a new boundary property is added or an existing property modified, the Boundary Property dialog pictured in Figure 8 appears. 

 

Figure 9. Boundary Property dialog.

The first selection in the dialog is the Name of the property. The default name is “New Boundary,” but you should change this name to something more descriptive of the boundary that is being defined.

 

The next selection is the BC Type drop list. This specifies the boundary condition type. Currently, BELA supports the following types of boundaries:

 

·        Fixed Voltage With this type of boundary condition, potential V is set to a prescribed along a given boundary

·        Mixed This denotes a boundary condition of the form:

                                                                                                                                                           (8)

The parameters for this class of boundary condition are specified in the Mixed BC parameters box in the dialog. By the choice of coefficients, this boundary condition can either be a Robin or a Neumann boundary condition.  By carefully selecting the  coefficient and specifying , this boundary condition can be applied to the outer boundary of your geometry to approximate an unbounded solution region. For more information on open boundary problems, refer to the Appendix.

·        Surface Charge Density  This selection is used to apply distributions of line charge to segments or arc segments in the problem.  Unlike all of the other boundary conditions, this BC type is often used on internal boundaries between materials or on isolated segments.  Typically, other BCs are only used on exterior boundaries.

·        Periodic  This type of boundary condition is applied to either two segments or two arcs to force the potential to be identical along each boundary. This sort of boundary is useful in exploiting the symmetry inherent in some problems to reduce the size of the domain which must be modeled. The domain merely needs to be periodic, as opposed to obeying more restrictive  or  line of symmetry conditions. Another useful application of periodic boundary conditions is for the modeling of “open boundary” problems, as discussed in Appendix 3. Often, a periodic boundary is made up of several different line or arc segments. A different periodic condition must be defined for each section of the boundary, since each periodic BC can only be applied to a line or arc and a corresponding line or arc on the remote periodic boundary.

·        Antiperiodic  The antiperiodic boundary condition is applied in a similar way as the periodic boundary condition, but its effect is to force two boundaries to be the negative of one another. This type of boundary is also typically used to reduce the domain which must be modeled, e.g. so that an electric machine might be modeled for the purposes of a finite element analysis with just one pole. 

3.7.3 Materials Properties

The Block Property dialog box is used to specify the properties to be associated with block labels. The properties specified in this dialog have to do with the material of which the block is composed. When a new material property is added or an existing property modified, the Block Property dialog pictured in Figure 10 appears.

 

Figure 10. Block Property dialog.

As with Point and Boundary properties, the first step is to choose a descriptive name for the material that is being described. Enter it in the Name edit box in lieu of “New Material.”

 

Next, permittivity for the material needs to be specified. BELA allows you to specify different relative permittivities in the vertical and horizontal directions (ex for the x- or horizontal direction, and ey for the y- or vertical .

 

A volume charge density (r) can also be prescribed by filling in the appropriate box in the material properties dialog.

 

3.7.4. Materials Library

Since one kind of material might be needed in several different models, BELA has a built-in library of Block Property definitions. The user can access and maintain this library through the Properties | Materials Library selection off of the main menu. When this option is selected, the Materials Library dialog pictured in Figure 11 appears. 

 

Figure 11. Materials Library dialog.

The options on this dialog allow the user to exchange Block Property definitions between the current model and the materials library.

 

The materials library should be located in the same directory as the BELA executable files, under the filename statlib.dat. If you move the materials library, Beladraw will not be able to find it.

3.7.5. Conductor Properties

The purpose of the conductor properties is mainly to allow the user to apply constraints on the total amount of charge carried on a conductor.  Alternatively, conductors with a fixed voltage can be defined, and the program will compute the total charge carried on the conductor during the solution process.

 

For fixed voltages, one could alternatively apply a Fixed Voltage boundary condition.  However, applying a fixed voltage as a conductor allows the user to group together several physically disjoint surfaces into one conductor upon which the total net charge is automatically computed.

 

The dialog for entering conductor properties is pictured in Figure 12.

 

Figure 12. Conductor Property dialog

3.8. Spawned Tasks

 

To actually mesh the model, analyze the model, and view the results, the Beladraw editor must spawn external programs. These tasks are most easily performed by the toolbar buttons pictured in Figure 13.

 

Figure 13. Toolbar buttons for spawning external tasks.

The first of these buttons (with the “yellow mesh” icon) runs the mesh generator. The solver actually automatically calls the mesh generator to make sure that the mesh is up to date, so you never have to call the mesher from within Beladraw. However, it is almost always important to get a look at the mesh and see that it “looks right.” When the mesh generation button is pressed, the mesher is called. While the mesher is running, an entry labeled “triangle” will appear on the Windows taskbar. Triangle is actually a console application that runs in a minimized window. After the geometry is triangulated, the finite element mesh is loaded into memory and displayed underneath the defined nodes, segments, and block labels as a set of yellow lines.

If you have a very large model, just keeping all of the mesh information in core can take up a significant amount of memory. If you are about to analyze a very large problem, it might be a good idea to choose the Mesh | Purge Mesh option off of the main menu. When this option is selected, the mesh is removed from memory, and the memory that it occupied is freed for other uses.

 

The second button, with the “hand-crank” icon, executes the solver, Belasolv.exe. Before Belasolve is actually run, the Triangle is called to make sure the mesh is up to date. Then, Belasolve is called. When Belasolve runs, it opens up a console window to display status information to the user. However, Belasolve requires no user interaction while it is running. When Belasolve is finished analyzing your problem, the console window will disappear. The time that Belasolve requires is highly dependent on the problem being solved. Solution times are typically on the order of 1 to 10 seconds, depending upon the size and complexity of the problem and the speed of the machine analyzing the problem.

 

The “big magnifying glass” icon is used to run the postprocessor once the analysis is finished. A detailed description of the postprocessor addressed in Section 4.

3.9. DXF Import/Export

For interfacing with CAD programs and other finite element packages, Beladraw supports the import and export of the AutoCAD dxf file format. Specifically, the dxf interpreter in Beladraw was written to the dxf revision 13 standards. Only 2D dxf files can be imported in a meaningful way.

 

To import a dxf file, select Import DXF off of the File menu. A dialog will appear after the file is seleted asking for a tolerance. This tolerance is the maximum distance between two points at which the program considers two points to be the same. The default value is usually sufficient. For some files, however, the tolerance needs to be increased to import the file correctly. Beladraw does not understand all the possible tags that can be included in a dxf file; instead, it simply strips out the commands involved with drawing lines, circles, and arcs. All other information is simply ignored. Generally, dxf import is a useful feature. It allows the user to draw an initial geometry using their favorite CAD package. Once the geometry is laid out, the geometry can be imported into Beladraw and detailed for materials properties and boundary conditions.

 

Do not despair if Beladraw takes a while to import dxf files (especially large dxf files). The reason that Beladraw can take a long time to import dxf files is that a lot of consistency checking must be performed to turn the dxf file into a valid BELA geometry. For example, large dxf files might take up to a minute or two to import.

The current Beladraw geometry can be exported in dxf format by selecting the Export DXF option off of the File menu. The dxf files generated from Beladraw can then be imported into CAD programs to aid in the mechanical detailing of a finalized design, or imported into other finite element or boundary element programs.


 

4. Postprocessor

The executable Belaview.exe is the postprocessor used to view solutions generated by the Belasolve solver. This program can either be run on its own, from the “Start” menu (to view previously solved problems), or spawned from within Beladraw to view a newly generated solution. Data files for Belaview have the .ans prefix.

4.1. Postprocessor modes

 

Similar to the preprocessor, the postprocessor always operates in one of three modes, depending upon the task to be performed. These modes are:  Point Values Mode In this mode, the user can click on various points in the solution region. Local field values are then listed in the Belaview Output window. Contour Mode This mode allows the user to define arbitrary contours in the solution region. Once a contour is defined, plots of field quantities can be produced along the contour, and various line integrals can be evaluated along the contour. Block Mode The Block Mode lets the user define a sub-domain in the solution region. Once the block has been defined, a variety of area and volume integrals can be taken over the defined sub-domain. Integrals include stored energy (inductance), various kinds of losses, total current in the block, and so on.

 

The current postprocessor mode is controlled via the Analysis Mode toolbar buttons, shown in Figure 14.

Figure 14. Analysis Mode toolbar buttons

The buttons denote, respectively, Point Values mode, Contour Mode, and Block Mode. The depressed button denotes the current mode. The default mode when Belaview starts is Point Values mode.

4.2. View and Grid Manipulation

The aspects of the current view and of the current grid are regulated via the use of toolbar buttons. The view is manipulated by the following toolbar buttons:

 

Figure 15. View Manipulation toolbar buttons.

and the grid settings are manipulated by these grid manipulation toolbar buttons:

 

Figure 16. Grid Manipulation toolbar buttons.

The grid and view manipulation work in exactly the same fashion as these same features in the preprocessor. Refer to Section 3.4 for a detailed description of grid manipulation, and to Section 3.3 for view manipulation.

4.3. Keyboard Commands

Unlike the preprocessor, Belaview is not very dependent on keyboard commands.

In the Point Values mode, there is only one relevant key press. In this mode, the Tab key allows the user to enter the coordinates of a specific point. After the coordinates of a point are entered, the field values at that point are displayed in the Belaview Output window.

 

In the Contours, there are three relevant keys. Pressing the Escape key wipes out any current contour or block definition. Pressing the Delete removes the last point added to the current contour or block edge. Last, pressing the Tab keys allows the user to numerically enter the coordinates of a point to be included in the current contour.

 

In the Block mode, the Escape and Delete keys have the same definition as in Contours mode. In the Block mode, the Tab key does not do anything, since all points on the contour must also be Points defining the model’s geometry.

4.4. Mouse Actions

In contrast, the operation of the postprocessor is very dependent upon input from the mouse.

 

In the Point Values mode, a Left Button Click is used to display field values at the current mouse location. If Snap to Grid is enabled, values are displayed at the closes grid point instead.

 

In the Contours mode, mouse clicks are used to define the contour. A Left Button Click adds the closest Point in the model’s geometry. Via a Right Button Click, the current mouse pointer location is added to the contour. A contour appears as a red line on the screen.

 

Blocks are defined in Block mode in a fashion very similar to the way in which contours are defined. A block is defined by drawing a contour around the region of interest. The contour appears as a green line on the Belaview screen. When the ends of the contour meet, the block is defined. All elements enclosed by the contour (all elements that form the block) then turn green in the Belaview window.

 

A Left Button Click attempts to add the nearest Point in the input geometry to the Block’s contour. However, a block can only be defined along line and arc segments from the input geometry. Each node on the boundary of the block must be selected in order to form the Block boundary. In Block mode, the Right mouse button has no function.

4.5. Contour Plot

One of the most useful ways to get a subjective feel for a solution is by plotting the eqipotentials of voltage. By default, a set of 19 equipotential lines are plotted when a solution is initially loaded into Belaview. The number and type of equipotential lines to be plotted can be altered using the Contours Plot icon in the Graph Mode section of the toolbar (see Figure 17). The Contour Plot icon is the icon with the black contours.

 

Figure 17. Graph Mode toolbar buttons.

When this button is pressed, a dialog pops up, allowing the choice of the number of contours.

 

In the contour plot dialog, a check box is also present titled “Show stress tensor mask”. If this box is checked, the contour lines associated with the last Weighted Stress Tensor integration are also displayed, by default as orange flux lines.

4.6. Density Plot

Density plots are also a useful way to get a quick feel for the flux density in various parts of the model. By default, a flux density plot is not displayed when Belaview first starts. However, the plot can be displayed by pressing the middle button in the Graph Mode section of the toolbar (see Figure 17). A dialog the pops up that allows the user to turn density plotting on.

 

The user can select between density plots of Voltage (V) or the magnitude of  Electric Field Intensity (E) or Electric Flux Density (D). The field at each point is classified into one of 12 contours distributed evenly between either the minimum and maximum flux densities or user-specified bounds.

4.7. Line Plots

When Belaview is in Contours Mode, various field values of interest can be plotted along the defined contour. A plot of a field value defined contour is performed by pressing the “graphed function” icon in the Plot and Integration group of toolbar buttons, shown in Figure 18.

Figure 18. Line Plot and Integration toolbar buttons

When this button is pressed, the X-Y Plot dialog (see Figure 19) appears with a drop list containing the types of line plots available. Choose the desired type of plot and press “OK.”

Figure 19. X-Y Plot dialog.

After “OK” is pressed, the program computes the desired values along the defined contour. These values are then plotted using the femmplot program, which is called automatically to display the plot.

 

By default, the Write data to text file box is not checked. If the user selects this option, the file selection dialog will appear and prompt for a filename to which to write the data. The data is written in two-column text format. If Write data to text file is selected, a femmplot window will not appear.

 

Currently, the type of line plots supported are:  Potential along the contour; Magnitude of the flux density along the contour; Component of flux density normal to the contour; Component of flux density tangential to the contour; Magnitude of the field intensity along the contour; Component of field intensity normal to the contour; Component of field intensity tangential to the contour;

 

In all of these plots, the direction of the normal is understood to be as shown in Figure20. The tangential direction is understood to be the direction in which the contour was defined. 

Figure 20. When in doubt plots and integrals taken on this side of a contour.

In certain cases, the quantity to be plotted can be ambiguous. This can occur, for example, if a plot of the tangential field intensity is requested on a contour running along an interface between two materials of differing permittivity. In this case, there is a discontinuity in the tangential field intensity, and the value of this quantity is different on each side of the interface. Belaview resolves the conflict by always evaluating the plots at a differentially small distance to the “normal” side of the line. Therefore, by defining the same contour but reversing the order in which the points are specified, plots of the quantity of interest on each side of a boundary can be obtained.

4.8. Line Integrals

Once a contour has been specified in Contours mode, Line Integrals can be performed along the specified contour. These integrals are performed by evaluating a large number of points at evenly spaced along the contour and integrating using a simple trapezoidal-type integration scheme.

 

To perform an integration, press the “integral” icon on the toolbar (as shown in Figure 0). A small dialog will appear with a drop list. Choose the desired integral from the drop list and press OK. The amount of time required to perform the integral will be virtually instantaneous for some types of integrals; however, some types may require several seconds to evaluate. When the evaluation of the integral is completed, the answer appears on the screen in a pop-up box.

 

The line integrals currently supported are:

 

·        E.t This integral returns the voltage drop along the defined contour

·        D.n. This integral returns the total electrix flux passing through a volume defined by extruding or sweeping the defined contour.  If this integral is performed over a closed contour, the resulting quantity is equal to the charge contained inside the contour.

·        Contour Length/Area. This integral returns the length of the defined contour in meters, as well as the area of the extruded or swept volume associated with the defined contour.

·        Force from stress tensor. This integral totals the force produced on the contour derived from Maxwell’s stress tensor. Deriving meaningful force results requires some care in the choice of integration path; refer to Section 4.10 for a detailed discussion of force and torque calculation.

·        Torque from stress tensor. This selection integrates the torque about the point (0,0) inferred from Maxwell’s stress tensor. Again, some guidelines must be followed to get accurate torque results (see Section 4.10).

 

4.9. Block Integrals

To select the regions over which a block integral is to be performed, left-click with the mouse in the desired region.  The selected region will appear highlighted in green.  For some block integrals (i.e. weighted stress tensor force and torque), one desires to select a conductor composed of lines and points, rather than a block.  In this case, the desired conductor can be selected by clicking on it with the right mouse button.  The selected conductor will appear in red.

 

To perform an integration, press the “integral” icon on the toolbar (as shown in Figure 18), and a dialog will appear with a drop list. Choose the desired integral from the drop list and press OK. The integral is then performed by analytically integrating the specified kernel over each element in the defined region, and summing the results for all elements.  Volume integrals may take several seconds to evaluate, especially on dense meshes. Be patient. When the evaluation of the integral is completed, the answer appears on the screen in a pop-up box.

 

The block integrals currently supported are: 

 

 

If the user is interested in the contours along which the integral was performed, the "stress tensor mask" box can be checked in the contour plot dialog. A set of orange (by default) lines will be displayed that.

4.10. Force/Torque Calculation

Often, the calculation of electrostatic and torques is a goal of finite element analysis. This section discusses some of the different methods of deducing forces and torques using BELA.

4.10.1. Lorentz Force/Torque

If one is attempting to compute the force on a collection of currents in a region containing only materials with a unit relative permeability, the volume integral of Lorentz torque is always the method to employ. Lorentz force results tend to be very accurate. However, again, they are only applicable for the forces on conductors of with unit permeability (e.g. coils in a voice coil actuator).

4.10.2. Weighted Stress Tensor Volume Integral

This volume integral greatly simplifies the computation of forces and torques. Merely select the blocks or conductors upon which force or torque are to be computed and evaluate the integral. No particular “art” is required in getting good force or torque results (as opposed to the Stress tensor line integral), although results tend to be more accurate with finer meshing around the region upon which the force or torque is to be computed.

 

One limitation of the Weighted Stress Tensor integral is that the regions upon which the force is being computed must be entirely surrounded by air and/or abutting a boundary. In cases in which the desired region abuts a non-air region, force results may be deduced from differentiation of stored electric field energy.

4.10.3. Maxwell Stress Tensor Line Integral

Generally, you should not use the Stress Tensor line integral to compute forces and torques if it can be avoided (i.e. use the volume integral version instead). The indiscriminate use Maxwell’s Stress Tensor can result in bad predictions forces and torques.

 

Maxwell’s stress tensor prescribes a force per unit area produced by the electric field on a surface. The differential force produced is:

(1/2)(D(E.n)+E(D.n)-n(D.E))
 (9)

where n denotes the direction normal to the surface at the point of interest. The net force on an object is obtained by creating a surface totally enclosing the object of interest and integrating the stress over that surface.

 

For best results, never integrate the stress tensor along an interface between materials. Always define the integration contour as a closed path around the object of interest with the contour displaced several elements (at least two elements) away from any interfaces or boundaries. 

 

Always use as fine a mesh as possible in problems where force results are desired.   A good way to proceed in finding a mesh that is “dense enough” is to solve the problem on progressively finer meshes, evaluating the force on each mesh. By comparing the results from different mesh densities, you can get and idea of the level of accuracy (by looking at what digits in the answer that change between various mesh densities). You then pick the smallest mesh density that gives convergence to the desired digit of accuracy.

 

For torque computations, all the same rules apply as for force computations (i.e. define integration contours away from boundaries and interfaces, and use a dense mesh).

4.11. Exporting of Graphics

Ultimately, you may want to export graphics from BELA for inclusion in reports and so on. It is possible to get what you are seeing on the screen onto disk in several different graphics formats.

 

Probably the easiest way to get graphics out of Belaview is to use the Copy as Bitmap or Copy as Metafile selections off of the main menu’s Edit list. These command takes whatever is currently in the Belaview window and copies it to the clipboard as a Device Independent Bitmap (.bmp format) and Extended Metafile (.emf format), respectively. The clipboard data can then be pasted directly into most applications (e.g. Word, MS Paint, etc).

 

LATEX afficionados typically find PostScript to be the most useful type of graphics output. BELA does not support postscript output directly, but it is still relatively easy to create postcript figures with BELA. To obtain a postscript version of the current view, you first must set up a postscript printer driver that outputs to File:. This is done via the following steps:  Choose +Settings/Printers+ off of the Windows Start menu. A window containing the list of currently defined printers will appear. Double click on the Add Printer icon in this list. The Add Printer Wizard will appear on the screen. Choose Local Printer; hit Next; A list of printers will appear. Choose a postscript printer off of this list. The Apple Laserwriter II NT is a good choice. Select FILE: as the port which will be used with this printer. Accept the defaults for all remaining questions.  Now, when you want a postscript picture of the currently displayed screen, just choose +File/Print+ off of Belaview’s main menu. As the printer, choose the postscript printer that you have previously defined. When you print to this printer, you will be prompted for a file name, and graphic will be written as a postscript figure to the specified file name.

4.12. Conductor Results

If conductor properties are used to specify the excitation, a useful byproduct is ready access to the voltage and charge on the conductor. To view the conductor results, select View|Conductor Props off of the Belaview main menu. A dialog, as pictured in Figure 21 will appear. There is a drop list on the dialog, from which the user selects the conductor for which results are desired. When a conductor is selected, the voltage and charge associated with that conductor are displayed.

 

Figure 21. Conductor results dialog.

4.13. Miscellaneous Useful View Commands

There are some additional entries on the Belaview View menu that might be useful to you from time to time. These are:

 

·        Smoothing  By default, a smoothing algorithm is applied to the flux density solution. Because first-order triangles are used as trial functions for potential, the resulting flux density and field intensity distributions are piece-wise constant in each element. The smoothing algorithm uses a nearest neighbor interpolation to obtain linear D and E distributions over each element. The smoothed solution generally looks better on the screen, and somewhat increases the accuracy of D and E near the vertices of each element. However, if you want to toggle smoothing, this can be done by selecting the Smoothing option.

·        Show Points  Especially when making graphics for reports, presentations, etc, it may be desirable to hide the small boxes on the screen that denote input node points. The Show Points option allows the user to toggle whether or not the input points are shown.

·        ToolBar  Use this toggle to hide and show the floating toolbar.

·        Point Props Use this toggle to hide and show the floating dialog box used to display point property information. 


 

5. Lua Scripting Documentation

 

5.1. What Lua Scripting?

The Lua extension language has been used to add scripting/batch processing facilities. The preprocessor and postprocessor can either run Lua scripts through a selection on the Files menu, or Lua commands can be entered in directly to the Lua Console Window in either program.

Lua is a complete, open-source scripting language. Source code for Lua, in addition to detailed documentation about programming in Lua, can be obtained from the Lua homepage at http://www.lua.org. Because the scripting files are text, they can be edited with any text editor (e.g. notepad).

In addition to the standard Lua command set, a number of BELA-specific functions have been added for manipulating files in both the pre- and post-processor. These commands are described in the following sections.

Lua scripts are invoked by selecting the Open Lua Script selection off of the File menu of either the pre- or post-processor. A file selection dialog then appears, and the selected Lua script file is executed.

5.2. Preprocessor Lua Command Set

A number of different commands are available in the preprocessor. Two naming conventions can be used: one which separates words in the command names by underscores, and one that eliminates the underscores. A list of alternate, equivalent preprocessor scripting function names is shown in Table 3.

 

open_bela_file    openbelafile
save_bela_file    savebelafile
create_mesh    createmesh
show_mesh    showmesh
purge_mesh    purgemesh
prob_def    probdef
analyse    analyze
run_post    runpost
add_node    addnode
add_block_label    addblocklabel
add_segment    addsegment
add_arc    addarc
select_node    selectnode
select_label    selectlabel
select_segment    selectsegment
select_arcsegment    selectarcsegment
clear_selected    clearselected
set_node_prop    setnodeprop
set_block_prop    setblockprop
set_segment_prop    setsegmentprop
set_arcsegment_prop    setarcsegmentprop
delete_selected    deleteselected
delete_selected_nodes    deleteselectednodes
delete_selected_labels    deleteselectedlabels
delete_selected_segments    deleteselectedsegments
delete_selected_arcsegments    deleteselectedarcsegments
zoom_natural    zoomnatural

zoom_out    zoomout
zoom_in    zoomin

add_material    addmaterial

add_point_prop    addpointprop
add_conductor_prop    addconductorprop
add_bound_prop    addboundprop
modify_material    modifymaterial
modify_bound_prop    modifyboundprop
modify_point_prop    modifypointprop
modify_conductor_prop    modifyconductorprop
delete_material    deletematerial
delete_bound_prop    deleteboundprop
delete_conductoruit    deleteconductoruit
delete_point_prop    deletepointprop
move_rotate    moverotate
move_translate    movetranslate
copy_rotate    copyrotate
copy_translate    copytranslate
set_edit_mode    seteditmode
select_group    selectgroup
new_document    newdocument
save_bitmap    savebitmap
save_metafile    savemetafile
exit_pre    exitpre
add_bh_point    addbhpoint
clear_bh_points    clearbhpoints
refresh_view    refreshview
message_box    messagebox
grid_snap    gridsnap
show_grid    showgrid
hide_grid    hidegrid
set_grid    setgrid

Table 3. Alternate preprocessor scripting function names

5.2.1. Object Add/Remove Commands

5.2.2. Geometry Selection Commands

5.2.3. Object Labeling Commands

5.2.4. Problem Commands

run_post("c:\\my-lua-script.lua","-lua-var=filename=myfilename")

All variables are passed as strings so lua internal conversion routines must be employed to obtain numbers from the strings. The postprocessor window can also be minimized by including the parameter "-windowhide", e.g.:

run_post("c:\\myluascrip.lua","-windowhide")

 

5.2.5. File Commands

runpost("\"c:\\program files\\femm30\\bin\\testpost.lua\"")

5.2.6. Mesh Commands

5.2.7. Editing Commands

This command will affect all subsequent uses of the other editing commands, if they are used WITHOUT the editaction parameter. 

5.2.8. Zoom Commands

5.2.9. View Commands

 

5.2.10. Object Properties

 

propnum           Symbol Description

0                      BlockName  Name of the material

1                      ex                   x (or r) direction relative permittivity

2                      ey                   y (or z) direction relative permittivity

3                      qs                   Volume charge

 

 

propnum       Symbol             Description
0                      BdryName                Name of boundary property
1                      Vs                               Fixed Voltage

2                      qs                               Prescribed charge density

3                      c0                               Mixed BC parameter

4                      c1                                Mixed BC parameter

5                      BdryFormat           Type of boundary condition:

0 = Prescribed V
1 = Mixed
2 = Surface charge density
3 = Periodic
4 = Antiperiodic

 

 

 

propnum       Symbol             Description
0                      PointName  Name of the point property
1                      Vp                   Prescribed nodal voltage

2                      qp                    Point charge density in C/m

 

 

propnum       Symbol             Description
0                      ConductorName    Name of the conductor property
1                      Vc                                Conductor voltage
2                      qc                                Total conductor charge
3                      ConductorType    0 = Prescribed charge,

1 = Prescribed voltage

5.2.11. Miscellaneous

5.3. Post Processor Command Set

 

There are a number of Lua scripting commands designed to operate in the postprocessor. As with the preprocessor commands, these commands can be used with either the underscore naming or with the no-underscore naming convention. The equivalent function names in the two conventions are shown in Table 4.

 

get_point_values    getpointvalues
exit_post    exitpost
add_contour    addcontour
bend_contour    bendcontour
clear_contour    clearcontour
line_integral    lineintergral
select_block    selectblock
group_select_block    groupselectblock
clear_block    clearblock
block_integral    blockintergral
zoom_natural    zoomnatural
zoom_in    zoomin
zoom_out    zoomout
show_grid    showgrid
hide_grid    hidegrid
show_mesh    showmesh
hide_mesh    hidemesh
set_edit_mode    seteditmode

hide_density_plot    hidedensityplot
show_density_plot    showdensityplot
hide_contour_plot    hidecountourplot
show_contour_plot    showcountourplot
show_points    showpoints
hide_points    hidepoints
grid_snap    gridsnap
set_grid    setgrid
get_problem_info    getprobleminfo
save_bitmap    savebitmap
get_conductor_properties    getconductorproperties
save_metafile    savemetafile
refresh_view    refreshview
select_point    selectpoint
show_point_props    showpointprops
hide_point_props    hidepointprops
message_box    messagebox
make_plot    makeplot

Table 4. Alternate postprocessor scripting function names

5.3.1. Data Extraction Commands

 

type Integral

0         

            1         

            2          Contour length

            3          Force from stress tensor

            4          Torque from stress tensor

 

            This integral returns either 1 or 2 values, depending on the integral type, e.g.:

 

Fx, Fy = lineintegral(3)

 

 

type  Integral
0          Stored Energy

1          Block Volume

2          Average D over the block

3          Average E over the block

4          Weighted Stress Tensor Force

5          Weighted Stress Tensor Torque


Returns one or two floating point values as results, e.g.:

 

Fx, Fy = blockintegral(4)

 

Symbol    Definition

V         Voltage

Dx        x- or r- direction component of displacement

Dy                   y- or z- direction component of displacement

Ex                   x- or r- direction component of electric field intensity

Ey                   y- or z- direction component of electric field intensity

ex                   x- or r- direction component of permittivity

ey                   y- or z- direction component of permittivity

nrg                 electric field energy density

 

Example: To catch all values at (0.01,0) use

 

V,Dx,Dy,Ex,Ey,ex,ey,nrg= getpointvalues(0.01,0) 

 

 

PlotType  Definition

     0         V     (Voltage)

     1         D|   (Magnitude of flux density)

     2         D . n (Normal flux density)

     3         D . t (Tangential flux density)

     4         E|   (Magnitude of field intensity)

     5         E . n (Normal field intensity)

     6         E . t (Tangential field intensity)


Valid file formats are

 

FileFormat    Definition
0             Multi-column text with legend
1             Multi-column text with no legend
2             Mathematica-style formatting


 For example, if one wanted to plot V  to the screen with 200 points evaluated to make the graph, the command would be:
makeplot(0,200)
If this plot were to be written to disk as a metafile, the command would be:
makeplot(0,200,"c:temp.emf")
To write data instead of a plot to disk, the command would be of the form:
makeplot(0,200,"c:temp.txt",0)

 

5.3.2. Selection Commands

5.3.3. Zoom Commands

5.3.4. View Commands

–        legend Set to 0 to hide the plot legend or 1 to show the plot legend.

–        gscale Set to 0 for a colour density plot or 1 for a grey scale density plot.

–        upper_D Sets the upper display limit for the density plot.

–        lower_D Sets the lower display limit for the density plot.

–        type Sets the type of density plot.  A value of 0 plots voltage; 1 plots the magnitude of D, and 2 plots the magnitude of E

–        numcontours Number of equipotential lines to be plotted.

–        upper_V Upper limit for contours.

–        lower_V Lower limit for contours.

If numcontours is -1 all parameters are ignored and default values are used, e.g. show_contour_plot(-1)

5.3.5. Miscellaneous

save_bitmap("\"c:\\temp\\femm30\\bin\\screenshot.bmp\"")

 

 

 

6. Numerical Methods

For those of you interested in what's going on behind the scenes in the Belasolve solver, this section is meant as a brief description of the methods and techniques used by BELA. References are cited as applicable.

6.1. Finite Element Formulation

All elements were derived using variational formulations (based on minimizing energy, as opposed to Galerkin, least squares residual, and so on). Explanations of the variational approach for 2-D planar problems with first-order triangle elements are widely available in the literature ([6] in particular was referred to during the creation of BELA).

6.2. Linear Solvers

For all problems, variations of the iterative Conjugate Gradient solver are used. This technique is appropriate for the sort of problem that BELA solves, because the matrices are symmetric and very sparse. A row-based storage scheme is used in which only the nonzero elements of the diagonal and upper triangular part of the matrix are solved.

 

The preconditioned conjugate gradient (PCG) code is based on the discussion in [7]. Minor modifications are made to this algorithm to avoid computing certain quantities more than once per iteration. Although Silvester promotes the use of the Incomplete Cholesky preconditioner, it is not used in BELA, because it nearly doubles the storage requirements--for each element of the matrix stored, a corresponding element of the preconditioner must also be stored. Instead, the Symmetric Successive Over-Relaxation (SSOR) preconditioner, as described in [8], is used. The advantage of this preconditioner is that it is built on the fly in a simple way using only the matrix elements that are already in storage. In general, the speed of PCG using SSOR is said to be comparable to the speed of PCG with Incomplete Cholesky.

 

Before a problem is formulated, a node renumbering algorithm is employed. Although the conjugate gradient schemes work well without renumbering, the renumbering seems to roughly halve the solution time. There is an overall advantage to using the renumbering because the time required to perform the renumbering is small compared to the time required to run the Conjugate Gradient solver. Although there are many possible approaches to renumbering, BELA uses the Cuthill-McKee method as described in [2]. Although there are newer schemes that yield a tighter profile, Cuthill-McKee does a relatively good job and requires very little time to execute. The renumbering seems to speeds up the CG algorithm by reducing the error between the SSOR approximation of the matrix inverse and the exact matrix inverse. An interesting paper on the effect of the ordering of the unknowns on convergence in conjugate gradient methods is [9].

 

6.3. Field Smoothing

 

Since first-order triangles are used by BELA, the resulting solution for E obtained by differentiating V is constant over each element. If the raw D and E are used by the postprocessor, density plots of D and 2-D plots of field quantities along user-defined contours look “crusty”. Also, the values of D and E are not so accurate at points in an element away from the element's centroid.

 

The use of smoothing to recover the accuracy lost by differentiating V is known as “superconvergence”. There are quite a few researchers actively pursuing this area. Of the greatest interest to BELA are so-called “patch recovery” techniques. The basic idea is the the solutions for D are most accurate at the centroid of the triangular element (known as its Gauss Point). One desires a continuous profile of D that can be interpolated from nodal values, in the same way that voltage V can be represented. However, the “raw” solution of D is multi-valued at any node point, those values being the different constant values of D in each element surrounding the node point. The general approach to estimating the “true” value of D at any node point is to fit a least-squares plane through the values of D at the Gauss points of all elements that surround a node of interest, and to take the value of the plane at the node point's location as its smoothed value of D [10].

 

However, this approach to patch recovery has a lot of shortcomings. For the rather irregular meshes that can arise in finite elements, the least-squares fit problem can be ill-condition, or even singular, at some nodes in the finite element mesh. Furthermore, the superconvergence solution can actually be less accurate than the piece-wise constant solution in the neighborhood of boundaries and interfaces.

 

One can note that the patch recovery method is merely a weighted average of the flux densities in all of the elements surrounding a given node. Instead of a least-squares fit, BELA simply weights the values of flux density in each adjacent element's Gauss point with a value inversely proportional to the distance from the Gauss point to the node point of interest. Away from boundaries, the results seem to be nearly as good as a least-squares fit. At boundaries and interfaces, the smoothed solution is no worse than the unsmoothed solution.


 

[1]  M. Plonus, Applied electromagnetics. McGraw-Hill, 1978.

[2]  S. R. Hoole, Computer-aided analysis and design of electromagnetic devices, Elsevier, 1989.

[3]  J. D. Jackson, Classical electrodynamics, 2nd ed, Wiley, 1975.

[4]  S. McFee, J. P. Webb, and D. A. Lowther, “A tunable volume integration formulation for force calculation in finite-element based computational magnetostatics,” IEEE Transactions on Magnetics, 24(1):439-442, January 1988.

[5]  F. Henrotte, G. Deliege, and K. Hameyer, “The eggshell method for the computation of electromagnetic forces on rigid bodies in 2D and 3D,” CEFC 2002, Perugia, Italy, April 16-18, 2002.

[6]  P. E. Allaire, Basics of the finite element method, 1985.

[7]  P. P. Silvester, Finite elements for electrical engineers, Cambridge University Press, 1990.

[8]  C. A. Fletcher, Computational techniques for fluid dynamics, Springer-Verlag, 1988.

[9]  E. F. D'Azevedo, P. A. Forsyth, and W. Tang, “Ordering methods for preconditioned conjugate gradient methods applied to unstructured grid problems,” SIAM J. Matrix Anal. Appl., 12(4), July 1992.

[10]  O. C. Zienkiewicz and J. Z. Zhu, “ The superconvergent patch recovery and a posteriori estimates, part 1: the recovery technique,” International Journal for Numerical Methods in Engineering, 33:1331-1364, 1992.

[11]  Q. Chen and A. Konrad, “A review of finite element open boundary techniques for static and quasistatic electromagnetic field problems,” IEEE Transactions on Magnetics, 33(1):663-676, January 1997.

[12]  E. M. Freeman and D. A. Lowther, “A novel mapping technique for open boundary finite element solutions to Poissons equation,” IEEE Transactions on Magnetics, 24(6):2934-2936, November 1988.

[13]  D. A. Lowther, E. M. Freeman, and B. Forghani, “A sparse matrix open boundary method for finite element analysis,” IEEE Transactions on Magnetics, 25(4)2810-2812, July 1989.


 

Appendix:

Open Boundary Problems

Typically, finite element methods are best suited to problems with well-defined, closed solution regions. However, a large number of problems that one might like to address have no natural outer boundary. A prime example is a parallel plate capacitor in air. The boundary condition that one would like to apply is V=0 at infinity. However, finite element methods, by nature, imply a finite domain. Fortunately, there are methods that can be applied to get solutions that closely approximate the “open boundary” solution using finite element methods.

A.1. Truncation of Outer Boundaries

 The simplest, but least accurate, way to proceed is to pick an arbitrary boundary “far enough” away from the area of interest and declare V=0  on this boundary. According to [11], a rule of thumb is that the distance from the center of the problem to the outer boundary should be at least five times the distance from the center to the outside of the objects of interest. Truncation is the method employed by most finite element programs, because it requires no additional effort to implement.

 

The down side to truncation is that get an accurate solution in the region of interest, a volume of air much larger than the region of interest must also be modeled. Usually, this large region exterior to the area of interest can be modeled with a relatively coarse mesh to keep solution times to a minimum. However, some extra time and space is still required to solve for a region in which one has little interest.

A.2. Asymptotic Boundary Conditions

 A thorough review of open boundary techniques is contained in [11]. Perhaps the simple way to approximate an “open” boundary (other than truncation) described in [11] is to use asymptotic boundary conditions. The result is that by carefully specifying the parameters for the “mixed” boundary condition, and then applying this boundary condition to a circular outer boundary, the unbounded solution can be closely approximated.

 

Consider a 2-D planar problem in polar coordinates. The domain is a circular shell of radius ro in an unbounded region. As

r®¥, potential V goes to zero. On the surface of the circle, the vector is a prescribed function of q. This problem has an analytical solution, which is:

V(r,q) =

¥
å
m = 1 

 

vm

rm

cos(m q+am)

(10)

where the vm and am parameters are chosen so that the solution matches the prescribed potential on the surface of the circle.

One could think of this solution as describing the solution exterior to a finite element problem with a circular outer boundary. The solution is described inside the circle via a finite element solution. The trick is to knit together the analytical solution outside the circle to the finite element solution inside the circle.

From inspecting (10), one can see that the higher-numbered harmonic, the faster the magnitude of the harmonic decays with respect to increasing r. After only a short distance, the higher-numbered harmonics decay to the extent that almost all of the open-space solution is described by only the leading harmonic. If n is the number of the leading harmonic, the open-field solution for large, but not infinite, r is closely described by:

V(r,q) »

vn

rn

cos(n q+ an)

(11)

Differentiating with respect to r yields:

 

V

r

= -

nvn

rn+1

cos(n q+ an)

(12)

If (12) is solved for an and substituted into (11), the result is:

 

V

r

+

æ
ç
è

 

n

r

 

ö
÷
ø

V = 0

(13)

Now, (13) is a very useful result. This is the same form as the “mixed” boundary condition supported by FEMM. If the outer edge of the solution domain is circular, and the outer finite element boundary is somewhat removed from the area of primary interest, the open domain solution can be closely approximated by applying (13) the circular boundary.

To apply the Asymptotic Boundary Condition, define a new, mixed-type boundary condition. Then, pick the parameters so that:

 

c0

=

 

eo n

 ro

 

(14)

c1

=

0

 

 

where ro is the outer radius of the region in meters (regardless of the working length units), and eo = 8.85418781762e-012.  Note that eo is defined in the Lua implementation in both the pre- and post-processors as the global variable eo, which can be used in any script or edit box in the program.

Although the above derivation was specifically for 2-D problems, it turns out that when a similar derivation is performed for the axisymmetric case, the definition of the mixed boundary condition coefficients are exactly the same as (14).  One subtle difference, however, is that n = 1 in the axisymmetric case corresponds to the case in which there is a net charge with  (i.e. the geometry looks like a point charge when viewed from a distance), whereas the n = 1 case corresponds to a dipole charge distribution in 2D planar problems.  If charge is conserved in the geometry of interest in the axisymmetric case, one needs to use n = 2 in Eq. (14).  It should be noted that this is a departure from the magnetostatic case with the vector potential formulation (e.g. as implemented in FEMM) in which n = 1 corresponds to a dipole arrangement in the axisymmetric case.

Some care must be used in applying this boundary condition. Most of the time, it is sufficient to take n = 1 (i.e the objects in the solution region look like a dipole when viewed from a large distance). However, there are other in which the leading harmonic is something other than n = 1. You need to use your insight into your specific problem to pick the appropriate n for the leading harmonic. You also must put the objects of interest roughly in the center of the circular finite element domain to minimize the magnitude higher-order field components at the outer boundary.

Although the application of this boundary condition requires some thought on the part of the user, the results can be quite good.  Applying the absorbing boundary condition imposes no additional computing cost on the problem. The ABC is computationally no more time-consuming to apply than enforcing V = 0 at the outer boundary. Solution times for the PCG solver are equivalent in either case.

A.3. Kelvin Transformation

A particularly good approach to “open boundary” problems is the Kelvin Transformation, a technique first discussed in the context of computational electromagnetics in [12] and [13]. The strengths of this technique are:

·         the effects of the exterior region are, in theory, exactly modeled by this approach;

·         a sparse matrix representation of the problem is retained (unlike FEM-BEM methods, which give the same “exact solution” but densely couples together the boundary nodes).

·         requires no “special” features in the finite element solver to implement the technique, other than the ability to apply periodic boundary conditions.

The purposes of this note are to explain what the Kelvin transformation is derived and to show how it is implemented in the context of the FEMM finite element program.

In the “far field” region, the material is typically homogeneous (e.g. air and free of sources. In this case, the differential equation that describes potential V is the Laplace equation:

Ñ2 V = 0

(15)

If we write (15) in polar notation, A is described by:

 

1

r

 

r

 

æ
ç
è

r

V

r

 

ö
÷
ø

+

1

r2

 

2 V

¶q2

= 0

(16)

Assume that the “near field” region of the problem can be contained in a circle of radius ro centered at the origin. The far-field region is then everything outside the circle.

One approach to unbounded problems is to attempt to map the unbounded region onto a bounded region, wherein problems can more easily be solved. Specifically, we desire a way to transform the unbounded region outside the circle into a bounded region. One simple way to make such a mapping is to define another variable, R, that is related to r by:

R =

ro2

r

 

(17)

By inspecting (17), it can be seen that this relationship maps the exterior region onto a circle of radius ro.

The next step is to transform (15), the differential equation that the field must satisfy, into the mapped space. That is, (15) must be written in terms of R and q rather than r and q. We can evaluate derivatives in terms of R instead of r by employing the chain rule:

 

r

=

R

 

æ
ç
è

 

d R

d r

 

ö
÷
ø

= -

R

 

æ
ç
è

 

R

ro

 

ö
÷
ø

2

 

 

(18)

 

Now, we can note that at r = R = ro,

 

 

V

r

= -

V

R

 

(19)

 

and we can substitute (18) into (16) to yield, after some algebraic manipulation:

 

 

1

R

 

R

 

æ
ç
è

R

V

R

 

ö
÷
ø

+

1

R2

 

2 V

¶q2

= 0

(20)

Eq. (20), the transformed equation for the outer region, has exactly the same form as inner region, only in terms of R rather than r. The implication is that for the 2-D planar problem, the exterior can be modeled simply by creating a problem domain consisting of two circular regions: on circular region containing the items of interest, and an additional circular region to represent the “far field.” Then, periodic boundary conditions must be applied to corresponding edges of the two circles to enforce the continuity of V  at the edges of the two regions. For a finite element formulation consisting of first-order triangles, (19) is enforced automatically at the boundaries of the two regions as long as the values of V at the boundaries are slaved together. The second circular region exactly models the infinite space solution, but does it on a bounded domain-one could always back out the field for any point in space by applying the inverse of (17).