Openfoam Simulation For Electromagnetic Problems: Zhe Huang
Openfoam Simulation For Electromagnetic Problems: Zhe Huang
Openfoam Simulation For Electromagnetic Problems: Zhe Huang
Electromagnetic Problems
Zhe Huang
September, 2010
Cover:
[Simulation results from ANSYS software]
[Chalmers Reproservice]
Gothenburg, Sweden 2010
SUMMARY
Simulation is considered as one of the most important and cost efficient approach of
industry research and development. OpenFOAM is one simulation tool with manual
solver compilation ability and 3D calculation capability, used for instance for
computational fluid dynamics (CFD) [1].
This thesis work is based on the OpenFOAM rodFoamcase and the rodFoam solver
used for plasma arc welding simulations which calculates the magnetic field in air. The
thesis work extends the case and the solver to solve electromagnetic field problems for
more materials including copper, linear steel and permanent magnets with different
geometries. In the future, this work can be applied into the design procedures of
electromagnetic devices, like electrical machines.
Based on the Maxwell equations, two different formulations (the A-V formulation and
the A-J formulation) are derived to solve magnetrostatic field problems. Formulations are
compiled manually into OpenFOAM solvers according to mathematic models by specific
program codes. Furthermore, force calculation equations are derived with the Lorenz
Force Method and the Maxwell Stress Tensor method. Then, simple geometries with
specific initial values and boundary conditions are calculated to test the new solvers.
The results of the developed simulation procedures in OpenFOAM are compared to the
results of another simulation software; COMSOL Multiphysics. It is found that the
simulation results show a very good agreement between OpenFOAM simulations and
COMSOL simulations.
Prediction of the future utilization of this thesis work and directions of development
are also discussed.
Contents
Symbols and Abbreviations ..................................................................... 6
1 Introduction .......................................................................................... 7
2 Theoretical Background ........................................................................ 8
2.1 Maxwell Equations .......................................................................................... 8
2.2 Derived Formulations for 3D Electromagnetic Field Problems ........................ 8
2.2.1 A-V formulations ......................................................................................... 9
2.2.2 A-J formulations ......................................................................................... 10
2.2.3 Equations for Magnetostatic Field Problems with Permanent Magnets........ 10
2.2.4 Transient State Formulations for Electromagnetic Field Problem ................ 10
2.3 Formulations for Force Calculations of 3D Problems..................................... 11
2.3.1 Lorenz Force Method ................................................................................. 11
2.3.2 Maxwell Stress Tensor ............................................................................... 11
2.4 Brief Overview of OpenFOAM ..................................................................... 12
2.4.1 Competitive Strength of OpenFOAM ......................................................... 12
2.4.2 File Structure of OpenFOAM Cases ........................................................... 13
2.4.3 Solvers Compilation ................................................................................... 14
4.2 Simulation Results of Single Bar with Copper Ring Geometry ...................... 37
4.3 Simulation Results of the Single Bar with a Copper Ring and Permanent
Magnets .............................................................................................................. 37
4.4 Simulation Results Comparison of Force Calculation .................................... 39
4.4.1 Field Distribution Comparison between OpenFOAM and COMSOL .......... 39
4.4.2 Numerical Results of Lorenz Force Method and Maxwell Stress Tensor ..... 40
Method................................................................................................................ 40
t
A
B
D
E
H
HC
I
J
L
R
S
V
r
0
divergence operator
curl operator
gradient operator
partial derivative
magnetic vector potential
magnetic flux density
electric displacement field
electric field intensity
magnetic field intensity
coercive magnetic field of permanent magnets
current of conductor
current density
total length of conductor
resistance
total area of conductor
electrical scalar potential
electric charge density
permittivity
magnetic permeability
relative magnetic permeability
permeability in the free space
magnetic reluctivity
relative magnetic reluctivity
electrical conductivity
V s/m
T
C / m2
V/m
A/m
A/m
A
A / m2
m
m2
V
C / m3
F/m
Wb / A m
S/m
1 Introduction
Since Energy Shortage and Greenhouse Effect are two of the most severe
problems nowadays, energy efficiency and clean energy are becoming more and
more important. Further, simulation is considered as one of the most important and
cost efficient approach of industry research and development. OpenFOAM is one
simulation tool with manual solver compilation ability and 3D calculation capability,
used for instance for computational fluid dynamics (CFD) [1].
This thesis work aims at expanding the calculation range of OpenFOAM, by using
C++ syntax in OpenFOAM, in order to solve electromagnetic field problems, which
can be applied into the design procedures of electromagnetic devices, like electrical
machines. The work is one part of electromagnetic devices design performed by the
vehicle company Volvo, and will be one step towards reaching the goal of the
highest electrical efficiency in electrical machines.
In order to design the optimized geometry and predict the operation behaviour, an
accurate knowledge of the field quantities inside the magnetic circuit is necessary
[2]. OpenFOAM is a powerful calculation tools for obtaining complex PDE
solutions. It has its own superiority compared to other commercial software, which
is not only on its free-of-charge and open source licence, but also on its numerous
solvers, utilities and libraries [1]. Compared to the black-box solvers in other
commercial softwares, new physics-mathematics models can be compiled into new
solvers based on these features of OpenFOAM.
With implemented methods of the finite element method (FEM), the finite volume
method (FVM) and the finite area method (FAM), OpenFOAM is the majority
investigated simulation software of this thesis work aimed to enlarge the calculation
range to calculate electromagnetic fields. Meanwhile, some of the OpenFOAM cases
are simulated in another commercial software COMSOL with the same initial field
values and boundary settings. This is done in order to check the validity of the edited
solvers with different mathematic models in the OpenFOAM simulations.
These new solvers can be applied into any electromagnetic devices such as an
electrical machine, in order to obtain the lowest cost and optimal design.
Chapter 2 describes different ways of formulating static electromagnetic field
problems according to Maxwell Equations. Based on the obtained electromagnetic
field, two electromagnetic force calculation methods are discussed. Besides, the
structure of OpenFOAM files is briefly illustrated.
Chapter 3 gives four practical examples of electromagnetic field and force
calculations by OpenFOAM. Each calculation contains the solver compilation and
pre-processing, which includes mesh generation, initial values and boundary
condition application, and post-processing.
Chapter 4 shows the results of all applications in chapter 3 and also discusses pros
and cons of OpenFOAM calculations.
Chapter 5 summarises the thesis work and predicts its future utilizations.
Meanwhile other interesting research trends about this topic are also discussed.
2 Theoretical Background
2.1 Maxwell Equations
The complete set of Maxwell equations is shown as below:
H = J +
D
t
B = 0
E =
(2.1)
(2.2)
B
t
D =
(2.3)
(2.4)
where H is the magnetic field intensity, J is the current density, D is the electric flux
density, B is the magnetic flux density, t is the time, E is the electric field intensity,
and is the electric charge density. There is interdependency between all variables
in the equations and, therefore, a unique solution [3].
For stationary or quasi-stationary electromagnetic field distribution cases, the
D
displacement current density term
is neglected, which gives the equation as
t
below
H = J
(2.5)
Since the divergence of the curl of any vector field in three dimensions is equal to
zero, the following equation is obtained
J = ( H ) = 0
(2.6)
The macroscopic material properties are defined by the constitutive relations [3]:
D = E
(2.7)
B = H = r 0 H
(2.8)
J = E
(2.9)
(2.10)
E=
A
V
t
(2.11)
According to (2.10), the magnetic flux density B can be calculated if the magnetic
vector potential A is known. A is defined in the whole problem region. While for the
electric scalar potential V, at the two ends of conductor region, a high value of V and
a low value of V are defined separately. Also at the non-conducting region, the value
of V is zero. Further, the permeability is assumed constant in each material, thus
non-linearity of for instance electrical machine core parts are not considered.
For magnetostatic cases, the displacement magnetic vector potential term
A
is
t
(2.12)
Combining (2.5) and (2.8) and (2.9) and (2.10) and (2.11), we will get
B
1
H = ( ) = ( A)
(2.13)
= J = E = (V )
( A) = V
(2.14)
(2.15)
the A-V formulation in the static case will get the unique solution of A obtained as
below
1
A A + V = 0
(2.16)
where the second term is a penalty term with the Coulomb gage. Due to the
addition of the Coulomb gauge, the solenoidality of the flux density must be
satisfied separately [4]. This is done by combining equations (2.6), (2.9) and (2.12)
into
(V ) = 0
(2.17)
Equations (2.16) and (2.17) are compiled into an OpenFOAM solver and is applied
for the current-carrying single bar case, which will be further discussed in chapter 3.
Another type of formula used for static magnetic field calculations is the A-J
formulation which uses the electric current density J as driving force.
Combining (2.5) and (2.8) and (2.10), we will get
B
1
H = ( ) = ( A) = J
(2.18)
Using the Coulomb Gauge similar to (2.16), the A-J equation is descripted as
1
A A = J
(2.19)
Compared to the A-V formulation, the A-J formulation does not have V terms,
which means it is more time-saving when applying FEM or FVM in the calculation
softwares.
2.2.3 Equations for Magnetostatic Field Problems with Permanent Magnets
When different types of magnetic materials are utilized as driving source, the
magnetostatic equations below can be applied to calculate the electromagnetic fields,
where H C is the coercive magnetic field of permanent magnets.
1
A A H C = 0
(2.20)[5]
1
1
A A H C = J
(2.21)
A
A A H C +
+ V = 0
t
(2.22)[5]
As it is discussed before, due to the Coulomb Gauge [4] eqation (2.15), the
solenoidality of the flux density must be satisfied separately. If formulas (2.6), (2.9)
and (2.11) are combined, equation (2.23) will be obtained.
A
V = 0
(2.23)
The Lorenz Force method can calculate the total force on a current-carrying
conductor by volume integration of the force density acting on each differential
current-carrying element. The Lorenz force can be calculated by
F = V ( J B)dV
(2.24)
The Maxwell Stress Tensor method can calculate force acting on a body by
integrating force density over a surface enclosing the body. The expression of the
Maxwell Stress Tensor is [6]:
Tij =
( Bi B j
1 2
B ij )
2
(2.25)
where i and j can be replaced by (x, y, z) , ij = 1 when i=j and ij = 0 when i j.
Therefore the Maxwell Stress Tensor can be expressed in another way as (2.26)
[6]:
Txx
T = T yx
Tzx
Txy
T yy
Tzy
1
1 2
2
( B x ) 2 B
Txz
1
T yz =
B B
y x
Tzz
1
[B z Bx ]
[B B ]
1
1
(B y ) 2 B 2
2
1
B z By
B B
y z
1
1 2
2
(
B
)
B
z
2
1
[B x Bz ]
(2.26)
The force calculation formula in terms of the force density vector is given by
F = S (
B ( B n)
1 2
B n)dS = S TdS
2
(2.27)[6]
Basic CFD
Incompressible flows
Compressible flows
Multiphase flows
Direct numerical simulations (DNS) and large eddy simulations (LES)
Combustion
Heat transfer and buoyancy-driven flows
Particle-tracking flows
Molecular dynamics
Direct simulation Monte Carlo
Electromagnetics
Stress analysis of solids
Finance
This thesis work is focused on extending the Electromagnetics solver by using
C++ syntax in OpenFOAM.
In OpenFOAM, there is a top-level code which is applied in the .C file and can
represent PDEs directly and flexibly.
For example, if the governing equation is given as equation (2.11):
E=
A
V
t
(2.11)
where fvm stands for Finite Volume Method and returns an fvMatrix, fvc stands
for Finite Volume Calculus and returns a geometricField.
Table 2.1 shows some of the correspondence between the mathematic expression
and the OpenFOAM expression [8]:
Table 2.1 Correspondence between mathematic expression and OpenFOAM expression [8]
Term
description
Mathematical expression
fvm::/fvc:: functions
Laplacian
laplacian(Gamma,phi)
laplacian(phi)
Time derivative
/ t
ddt(phi)
Gradient
grad(phi)
Divergence
( )
div(psi)
Curl
curl(phi)
3. Implementation in OpenFOAM
3.1 Modelling Methods
With the implemented methods of the finite element method (FEM), the finite
volume method (FVM) and the finite area method (FAM), OpenFOAM is the
majority investigated simulation software of this thesis work used to calculate
electromagnetic fields. Meanwhile, some of the OpenFOAM cases are simulated in
another commercial software; COMSOL with the same initial field values and
boundary settings, in order to check the validity of the edited OpenFOAM solvers.
This thesis work is based on the OpenFOAM rodFoamcase and the rodFoam
solver used for plasma arc welding simulations which calculates the magnetic field
in air. The thesis work extends the case and the solver to solve the electromagnetic
field problems for more materials including copper, linear steel and permanent
magnets with different geometries.
Conductor
Air box
x(m)
0.4
4
y(m)
0,4
4
z(m)
6
6
left
down
up
front
back
right
y
z
x
Figure 3.1 Geometry of the single bar case
Solver name:
SingleBarFoam1 (with A-V solver)
12
13
10
1
0
11
y
7
15
14
Figure 3.2 Patches and points defined on the x-y plane (z=0)
Generally, there are four validity constraints used in OpenFOAM: points, faces,
cells and boundary faces. These are basic constraints that a mesh must satisfy. As it
is shown in the figure 3.2, 16 points and 9 faces are defined for the front surface on
the x-y plane (z=0). For the whole 3 dimensional geometry, the length along the zdirection is 6 m. Furthermore, the whole length along the z direction is divided into
10 units.
The OpenFOAM BlockMesh utility is used to generate meshes in the file
constant /polyMesh/ blockMeshDict of the casedirectory. Inside this file, cell shape
and mesh density can be defined. Hexahedral, wedge, pyramid, tetrahedron, and
tetrahedral wedge cell shapes can be defined by ordering of point labels in
accordance with the numbering scheme contained in the shape model. Hexahedral
cell shape is used in this thesis work. One example of geometry definition is given in
the codes below [8]:
8
(
0
0
1
1
0
0
1
1
0)
0)
0)
0)
0.5)
0.5)
0.5)
0.5)
//point 0
//point 1
//point 2
//point 3
//point 4
//point 5
//point 6
//point 7
Name
Description
Patch
generic patch
symmetryPlane
plane of symmetry
Empty
Wedge
Cyclic
cyclic plane
Wall
For the 3D single bar problem,, patch is used for every boundary. The boundary
type wedge and empty are used for the 2D axis symmetric cases later on. One
example of boundary file is given in Appendix 2.
Field initial values of the boundaries are defined in the so called 0 directory. For
the cases using the A-V formulation, the magnetic vector potential A, the electrical
scalar potential V and the electrical conductivity are pre-defined in the 0
directory. Consequently, for the cases using the A-J formulation, the magnetic
vector potential A and the current density J are pre-defined in the 0 directory. The
field initial values which need to apply in the inner region of geometries, such as
electrical conductivity in the cases which use the A-V solver and current density J
in the cases which use the A-J solver cases, are defined in file system/setFieldDict.
Table 3.3 shows the boundary conditions and initial values of the single bar 3D
problem when using the A-V formulation. The naming of boundaries are shown in
Fig. 3.1.
Table 3.3 Boundary conditions and initial values of the Single Bar 3D problem when using the A-V
formulation
Boundary names
(see Fig. 3.1)
Air left and right
and back and
front
Air up and down
air metal
Fixed value 0
1.26e-6
Zero gradient
Zero gradient
1.26e-6
Zero gradient
1.26e-6
670 V
1.26e-6
0V
Conductor up
Zero gradient
Conductor down
Zero gradient
Regarding figure 3.1, the left, right, back and front patches of the air box should
be infinitely far from the cable where the magnetic vector potential A is set to 0. The
boundary condition zero gradient represents that the normal gradient of the field is
zero. Since the values of the electrical conductivity are almost zero in air but quite
large in metal, 1e-5 S/m is chosen to represent air and 2700 S/m for the metallic
material. Furthermore, since the magnetic permeability of air in some cases is
almost equal to the value of metal, such as copper, 1.26e-6 Wb / A m is chosen in
this simulation.
In order to have the same driving force with the A-J formulation as it is using the
A-V formulation (which applies 670V) a current density of 310500 A / m 2 is applied,
according to the calculation below:
U = I R = I
L
L
=J
S
i.e.
J=
U 670 2700
=
= 301500 A / m 2
L
6
where U is the voltage of the conductor, I is the current of the conductor, L is the
total length of the conductor, and S is the total area of the conductor.
Table 3.4 shows the boundary condition and initial values of the single bar 3D
simulation when the A-J formulation is applied.
Table 3.4 Boundary condition and initial values of the single bar 3D problem when
using the A-J formulation
A
Air left and right and
back and front
Air up and down
Conductor up
Fixed value 0
air copper
1.26e-6
Zero gradient
Zero gradient
1.26e-6
1.26e-6
Conductor down
Zero gradient
1.26e-6
J
Zero gradient
Zero gradient
Fixed value
(0,0,301500)
Fixed value
(0,0,301500)
A A + V = 0
(2.16)
The A-J formulation in the magneto static case is given as the following:
A A = J
Since
( F ) = ( F ) 2 F
(2.19)
(3.1)
the A-V formulation and the A-J formulation can be transformed into formula (3.2)
and (3.3):
2 A + V = 0
(3.2)
2 A = J
(3.3)
where muMag is the permeability (of air and the copper conductor).
Appendix 5 gives an example of the .C file of the A-J solver. This .C file
combined with the files in the Make directory (refer to figure 2.2) can be compiled
by the wmake command. After compilation, an executable file of the new solver is
generated. The executable file should be copied into the Case directory and it is
used for calculations with specific boundary conditions and initial values.
(5)Post-processing
Chapter 4.1 will show great similarity of calculation results from the OpenFOAM
A-V solver and A-J solvers applied to the same 3D geometry. Hence, if OpenFOAM
calculation results from the A-V solver show similarity to COMSOL calculation
results from this simulation test, the validity of the A-J solver in the 2D axissymmetric case should also be proved.
OpenFOAM needs 3D control volumes even for the 2D cases. However, the patch
type empty can be applied to generate 2D or even 1D geometries by specifying
special patches into empty in one dimension or in two dimensions. Meanwhile, the
patch type wedge can be used to specify a geometry as a wedge with a small angle
and being one cell thick along the plane of symmetry to obtain a 2D axi-symmetric
geometry, which is applied in the simulation below.
(1) Geometry properties for 2D axi-symmetric simulation applied in
OpenFOAM and COMSOL
R
R
Symmetry Axis
L
Air box
Cylindrical
conductor
y
x
r
r(inner)(m)
0.2
R(outer)(m)
2
L(m)
6
In OpenFOAM, the square meshes are still applied in the 2D axis symmetric
geometry; 10 cells are applied on the x direction and 30 cells are applied on the y
direction. Meanwhile, 1 cell unit is applied along the z direction to express the 2D
geometry. In COMSOL, the mesh generation is done by the default mesh type
Normal.
(3) Boundary Condition and Initial Field
As shown in table 3.6 and table 3.7, the current density is used as the driving force
in COMSOL and the equivalent voltage is applied for the A-V solver in
OpenFOAM. In COMSOL, the current density is set to 0.3 A / mm 2 along y direction.
Consequently, in OpenFOAM, a voltage value of 670V is applied to the cylindrical
conductor, in order to introduce a current density of 0.3 A / mm 2 .
Table 3.6 Initial values for single bar simulation with 2D axis symmetric geometry in COMSOL
copper
air
(H /m)
1.26e-6
1.26e-6
J ( A / mm 2 )
0.3
0
(S/m)
2700
1e-5
Table 3.7 Initial values for single bar simulation with 2D axis symmetric geometry in OpenFOAM
( H / m)
copper
Air
1.26e-6
1.26e-6
U(V)
High: 670, Low: 0
0
(S/m)
2700
1e-5
There is one difference between OpenFOAM and COMSOL when they are used
for 2D axis symmetric simulation. In COMSOL, when 2D simulation is taken, 2D
geometry is drawn at the first place. In OpenFOAM, though 2D simulation is needed
to be done, 3D geometry is always needed to be drawn first. Then 2D axis
symmetric simulation is needed, the boundary conditions of two surfaces which are
parallel to the rotating axis are treated as type wedge. Table 3.8 shows the
boundary conditions applied in this single bar 2D axis symmetric simulation.
Table 3.8 Boundary conditions for single bar simulation with 2D axis symmetric geometry in
OpenFOAM
A( V s / m )
Fixed value 0
U(V)
Zero gradient
Wedge
Zero gradient
Zero gradient
Wedge
Fixed value 670
Fixed value 0
The simulation results from both soft-wares are discussed in chapter 4.1.
left
steel bar
back
down
up
front
copper conductor
x
right
Figure 3.4 Geometry of a steel cuboid bar with a ring-shape copper conductor
The geometry data of a steel cuboid bar with a ring-shape copper conductor is
shown as the table 3.9 below.
Table 3.9 Geometry data table of single bar with copper ring case
x(m)
y(m)
z(m)
Air Box
Steel Bar
0.14
0.14
Inner R
Outer R
0.2
0.21
Copper Ring
The mesh on x-y plate, for the simulation with one copper conductor and one
cuboid steel bar, is generated as it is shown in figure 3.5.
22
21
16
15
14
17
20
23
12
18
19
25 26
11
33 32 13
34
0 1
38 40 39
27 35 36 37 5
10
4
24
29
2
3
28
9
y
x
30
6
31
As it is shown in figure 3.5 of the up surface in the x-y plane, 41 points are
defined. For the whole 3D geometry, 82 points and 35 blocks are defined to
illustrate the steel bar rounded by the copper tube which is surrounded by air. The
mesh definition approach is similar to the approach applied in the single bar case.
3.3.3 Boundary Conditions and Initial Values
The current density, J, in the ring shaped coil (with a magnitude of unity) can
be expressed as below:
J =
y
x2 + y2
x
x2 + y 2
,0
In order to set the coil different permeability values, the funkySetField utility is
used. Compared to the utilization of the utility setFieldDict used in previous cases,
this utility can set non-uniform field values and non-uniform field application regions
by C++ commands whereas the utility setFieldDict can only set uniform fields. The
OpenFOAM C++ commands used instead of the mathematic expressions to set the
fields are given in the text below. The commands are put in the file
system/funkySetFieldsDict.
currentdensity2
{
field J;
//field initialise
[0 -2 0 0 0 1 0];
In the code above, field is used to declare one field, expression is used to write
field values, condition will select a subset of the cells, and dimensions describes
the dimension of the field in SI units.
In this simulation, the relative permeability is assumed to be constant (thus
saturation effects are ignored), and relative reluctivity can be expressed as in
formula (3.4),
0 r =
1
1
=
0 r
(3.4)
//field to initialise
expression 1;
[0 0 0 0 0 0 0];
//field initialise
[0 0 0 0 0 0 0];
In order to show the permeability difference between the coil and the steel, a high
permeability value (2000) is chosen for the steel bar. Table 3.10.1 and table 3.10.2
show the initial values and boundary conditions after executing funkySetField utility.
Table 3.10.1 Initial values for a copper coil with steel bar geometry:
J( A / m 2 )(absolute value)
r ( r )
Air
1 (1)
Coil
301500
1 (1)
Steel bar
2000 (0.0005 )
Table 3.10.2 Boundary conditions for a copper coil with steel bar geometry:
A( V s / m )
Fixed value 0
J( A / m 2 )
Zero gradient
Zero gradient
Zero gradient
Zero gradient
Zero gradient
Zero gradient
Zero gradient
Since the copper conductor only have inner boundaries in this simulation, there
is no outer boundary definition for the conductor.
A A H C = J
(2.21)
2 A = J H C
(3.5)
( 0 r A) ( 0 r A) H C = J
(3.6)
The geometry data is given in table 3.11 and seen in figure 3.7.
Table 3.11 Geometry data of a single bar and copper ring and a permanent magnet
x(m)
y(m)
z(m)
Air box
Steel bar
0.14
0.14
0.14
0.14
0.2
Inner R
Outer R
0.2
0.21
Permenant
magnet
Copper
conductor
Steel bar
Air box
z
Figure 3.7 Geometry of a single bar with a copper ring and a permanent magnet
The same mesh generation method as in the previous case (which includes a
cuboid steel bar and a copper conductor) is applied in this simulation. But the mesh
along the z dimension is increased from 10 to 60 in order to obtain accurate enough
calculation results for the permanent magnet of which the thickness is 0.2m.
3.4.2 Boundary Condition and Initial Values
Compared to the the previous case (the cuboid steel bar and a copper conductor),
the coercive magnetic field value of a permanent magnet is added by the following
command.
magnetization4
{
field M;
//field initialise
expression vector(724000,0,0);
//coercive magnetic field in permanent
magnets
condition pos().x>=-0.14 && pos().x<=0.14 && pos().y>=-0.1 && pos().y<=0.1 &&
pos().z>=1.3 && pos().z<=1.5;
keepPatches 1;
dimensions
[0 -1 0 0 0 1 0];
}
Since the copper conductor and magnet block only have inner boundaries in this
simulation, there is no outer boundary definition for the conductor and magnet block.
3.4.3 Solver Compilation
Air box
Bar1
Bar2
X(m)
2
0.4
0,4
Y(m)
2
0,4
0,4
Z(m)
10
10
10
In this simulation, 1A and -1A current with opposite directions are applied into the
two parallel bars separately. The current densities are calculated using the formula
below, where J is the current density, I is the current flowing through one conductor
and S is the total area of one conductor.
J=
1A
I
=
= 6.25 A / m 2
2
S 0.4 0.4m
Table 3.14.1 and table 3.14.2 show the initial values and boundary conditions of
the force calculation case.
Table 3.14.1 Initial value for force calculation case
J( A / m 2 )
H c ( A / m)
0 ( N / A2 )
Air box
(0, 0, 0)
1.26e-6
Bar 1
(0, 0, 6.25)
0.999 994
1.26e-6
Bar 2
(0, 0, -6.25)
0.999 994
1.26e-6
A( V s / m )
Fixed value 0
J( A / m 2 )
Zero gradient
H C ( A/m)
Zero gradient
Zero gradient
Zero gradient
Zero gradient
Fixed value
(0,0,6.25)
Fixed value
(0,0,-6.25)
Zero gradient
Zero gradient
Zero gradient
Zero gradient
One header file Force.H is applied to calculate Lorenz Force by executing volume
integration of Lorenz Force density. This header file is called by the major solver
file (.C file). The formulation of Lorenz Force calculation in the header file is given
as below:
F= Jtot^B
A loop of the volume integration for Lorenz Force density on one of the bars in
file Force.H is calculated using the method below, where dz1 is the length of one
cell in the z direction.
After 3D calculation by the A-V solver, the maximum values of current density
Je, vector potential A and magnetic flux density B are given in the table 4.1. Figure
4.1 gives the surface plots of the three parameters for the single bar geometry.
Table 4.1 Maximum values from 3D calculation for single bar geometry by A-V solver
A( V s / m )
0,0269804
Je( A / m 2 )
301500
Surface plot of current density Je
vector potential A
B(T)
0,0411579
magnetic flux density B
Figure 4.1 Surface plots from 3D calculation for single bar geometry by A-V solver
For this simulation, the current density J is one of the applied initial values,
and the calculated fields are vector potential A and magnetic flux density B. The
maximum values of A and B are given in table 4.2 and the surface plots of current
density, vector potential and magnetic flux density are shown in figure 4.2.
Table 4.2 Maximum values from 3D calculation for single bar geometry by A-J solver
A (V s / m )
0,0269804
B (T)
0,0411579
vector potential A
Figure 4.2 Surface plots from 3D calculation for single bar geometry by A-J solver
As seen in the OpenFOAM calculation results above, the maximum values and
field distribution of magnetic flux density are the same no matter if voltage or
equivalent current density is applied as driving source. In both of the cases, the
maximum values of magnetic flux density B are obtained in close proximity to the
edge of the conductor.
4.1.3 OpenFOAM and COMSOL Simulation Comparisons on a 2D Axis
Symmetric Geometry
Next section shows the simulation results from OpenFOAM and COMSOL about
the single bar geometry by the 2D axis symmetric method.
As discussed in chapter 3, the same values of magnetic permeability and
electrical conductivity are applied both in OpenFOAM and COMSOL. Current
density is applied in COMSOL and the equivalent voltage as driving source is
applied in OpenFOAM.
(a) The maximum field values calculated by both softwares are given in table 4.3 and
table 4.4 separately. Figure 4.3 uses bar chart to illustrate the difference of
maximum values between the OpenFOAM calculation and the COMSOL
calculation.
Table 4.3 Maximum field values from OpenFOAM
J( A / m 2 )
3.015e5
B (T)
0.0370033011
A (V s / m )
0.0224271994
J( A / m 2 )
3e5
B (T)
0.0377
A (V s / m )
0.0212
As shown in table 4.3 and table 4.4, the results from OpenFOAM and COMSOL
show a 5.7% difference in the magnetic vector potential A and 1.8% for the magnetic
flux density B. This comparison shows great similarity between the calculations from
OpenFOAM and COMSOL.
(b) Figure 4.4 and figure 4.5 show the field distributions of the magnetic vector
potential A and magnetic field density B:
Figure 4.4 Surface plots of magnetic vector potential and magnetic field density in OpenFOAM
Figure 4.5 Surface plots of magnetic vector potential and magnetic field density in COMSOL
As one can see here, the surface plots on magnetic vector potential A and magnetic
field density B from OpenFOAM and COMSOL show great similarity, which can
prove the accuracy of OpenFOAM calculation.
(c) Magnetic field distributions obtained from OpenFOAM and COMSOL
calculations are further compared by using a line plot approach along a line which
is parallel to the x axis and passes through point (0, 3). Figure 4.5 and figure 4.6
show the line plots obtained by the two softwares:
Figure 4.6 Line plot of magnetic flux density B along the x coordinate when y=3 in OpenFOAM (length
VS magnetic flux density B)
Figure 4.7 Line plot of magnetic flux density B along the x coordinate when y=3 in COMSOL (length VS
magnetic flux density B)
According to the results above, when current density values are around 3e5 A / m 2
the maximum value of magnetic flux density B is about 0.37T from both
OpenFOAM and COMSOL calculations.
From the distribution plot of flux density B from COMSOL calculation in figure
4.7, the magnetic flux density is increased from zero obtained at x=0m to the
maximum value at the rim of the conductor when x=0.2m. Later on, the flux density
B is decreased again from the rim of the conductor (when x=0.2m) to the rim of the
air box (when x=2m). This distribution of flux density B also confirms the
distribution of B obtained from 3D simulations above.
Comparing the field distribution plot of flux density B from the OpenFOAM
calculation in figure 4.6 to figure 4.7, the main distribution trends are the same, but
one calculation difference can be observed. In figure 4.6, the magnetic flux density
B is started from a small non-zero value at x=0 and is ended at x=2.25m. The reason
to this phenomenon in OpenFOAM can be the unit square mesh application. Thus,
the calculation procedure cannot distinguish the starting calculation point within one
cell distance. This problem should be solved by applying a finer mesh in the future
work.
One problem should be pointed out that, in 2D axis-symmetric geometry, the A-V
solver can solve the electromagnetic problems with low electrical conductivity
values but not the cases with high electrical conductivity. The reason of this problem
has not been defined. But since the A-J solver will not involve any electrical
conductivity , it doesnt have this issue at all.
Line Plot of
Magnetic Flux Density B
when Relative
Reluctivity=5e-4 in steel bar
Figure 4.8 Line plots comparison between cases with and without high permeability steel bar
z
y
x
Figure 4.9 Surface plots from COMSOL
In both surface plots from COMSOL in figure 4.9 and OpenFOAM in figure 4.10,
the comparisons between the figures at the left side and the figures at the right side
show that, the magnetic flux density B is greatly increased after adding permanent
magnetic material.
Furthermore, the comparison between the two plots on the right side in figure 4.9
and figure 4.10 shows better computing ability in OpenFOAM than COMSOL for
this 3D simulation case including a permanent magnet. In OpenFOAM, a suitable
mesh approach can be decided according to the complexity of one specific
geometry. For this simulation done by the author, COMSOL can only solve this
problem with coarser mesh because of the limitation of computer memory, which
gives the worse field distribution stream line as plot of 4.9 (to the right).
In figure 4.11, the surface plots of force density distribution calculated from
OpenFOAM and COMSOL are compared. For the convenience of the following
discussion, the terms inner edge and outer edge are defined as shown in figure
4.11. Both plots illustrate that the forces acting on the inner edges of two parallel
cables are larger than the forces acting on the outer edges. Further the shorter
distance between the two bars, the larger the magnetic force is.
surface plot of Lorenz Force density
in OpenFOAM
inner
edge
outer
edge
x dimension
left
right
Figure 4.11 Surface plots of force density distribution from OpenFOAM and COMSOL calculations
4.4.2 Numerical Results of Lorenz Force Method and Maxwell Stress Tensor
Method
Both the Lorenz force method and the Maxwell stress tensor method are applied
in OpenFOAM. The numerical values of the magnetic force from the OpenFOAM
calculation are given in table 4.5.
Table 4.5 Numerical values comparison of magnetic forces calculated by two different methods in
OpenFOAM
(-1.57206e-07 -7.64703e-23 0)
(1.57206e-07 -7.41412e-23 0)
(1.20408e-08 -1.2851e-23 0)
Force on surface 2
(-1.61322e-07 5.35406e-23 0)
Force on surface 3
(-3.98597e-09 -3.3557e-09 0)
Force on surface 4
(-3.98597e-09 3.3557e-09 0)
(-1.57253e-07 6.4851e-23 0)
(1.61322e-07 -3.90662e-23 0)
Force on surface 2
(-1.20408e-08 -1.44123e-23 0)
Force on surface 3
(3.98597e-09 -3.3557e-09 0)
Force on surface 4
(3.98597e-09 3.3557e-09 0)
(1.57253e-07 -7.94093e-23 0)
-1.55704e-07
1.557316e-07
Force(N/m)
As the results shows in figure 4.12, with opposite current direction flowing in the two
cables, the left side cable undertakes a force pointing to the left direction and vice versa.
In other words, a repulsive force is applied to both of the cables. And when 1A of current
is flowing though two parallel cuboid cables, the force between them is around 1.57e7N/m according to the OpenFOAM calculation. Compared to the COMSOL case with the
same geometry and boundary condition and initial values, the repulsive force is around
1.55e-7N/m. Figure 4.12 shows that the force calculation difference which acts on the
right or left cable is less than 0.1%, which shows great calculation accuracy in
OpenFOAM.
Since the force calculation is based on the first calculation results of electromagnetic
fields, new errors will be introduced which can be quite large. For the Maxwell Stress
Tensor method, the solution also depends on the choice of the integration path which will
bring calculation results difference for different paths [10].
No.
solver
Case
function
A-V solver
SingleBarCase3D1
A-J solver
SingleBarCase3D2
A-V solver
SingleBarCase2D1
A-J solver
SingleBarCase2D2
tbmCase12
2barsCase17
In general, a more complex mesh takes more simulation time. From the perspective of
simulation speed, OpenFOAM has its great advantage to handle parallel computing,
which can divide one task into several processors in order to decrease the total simulation
time. Along with more and more complex mesh, parallel computing utility should be
another research division of OpenFOAM.
Appendix 1
Mesh generation in file constant/polyMesh/blockMeshDict by blockMesh utility of
SingleBarCase3D1
/*---------------------------------------------------------------------------*\
| =========
| \\
/ F ield
| \\
/ O peration
| \\ /
A nd
M anipulation |
\\/
| Web:
http://www.openfoam.org
|
\*---------------------------------------------------------------------------*/
FoamFile
{
version
2.0;
format
ascii;
class
dictionary;
object
blockMeshDict;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(-0.2 0.2 0) //0
(0.2 0.2 0)
//1
//5
(2 -2 0) //6
(-2 -2 0) //7
(-2 -0.2 0) //8
(-2 0.2 0) //9
(2 0.2 0)
//10
(2 -0.2 0) //11
(-0.2 2 0) //12
(0.2 2 0) //13
(0.2 -2 0) //14
(-0.2 -2 0) //15
//17
//21
(2 -2 6) //22
(-2 -2 6) //23
(-2 -0.2 6) //24
(-2 0.2 6) //25
(2 0.2 6)
//26
(2 -0.2 6) //27
(-0.2 2 6) //28
(0.2 2 6) //29
(0.2 -2 6) //30
(-0.2 -2 6) //31
);
blocks
(
hex (0 3 2 1 16 19 18 17) (20 20 10) simpleGrading (1 1 1) //0
hex (1 2 11 10 17 18 27 26) (20 10 10) simpleGrading (1 1 1) //1
hex (2 14 6 11 18 30 22 27) (10 10 10) simpleGrading (1 1 1) //2
hex (3 15 14 2 19 31 30 18) (10 20 10) simpleGrading (1 1 1) //3
hex (8 7 15 3 24 23 31 19) (10 10 10) simpleGrading (1 1 1) //4
hex (9 8 3 0 25 24 19 16) (20 10 10) simpleGrading (1 1 1) //5
hex (4 9 0 12 20 25 16 28) (10 10 10) simpleGrading (1 1 1) //6
hex (12 0 1 13 28 16 17 29) (10 20 10) simpleGrading (1 1 1) //7
hex (13 1 10 5 29 17 26 21) (10 10 10) simpleGrading (1 1 1) //8
);
edges
(
);
patches
(
patch airleft
//patch 0
(
(8 7 23 24) //face 0
(9 8 24 25) //face 1
(4 9 25 20) //2
)
patch airdown
//patch 2
(
(4 20 28 12) //6
(12 28 29 13) //7
(13 29 21 5) //8
)
patch airup
//patch 3
(
(14 6 22 30)
//9
(15 14 30 31)
//10
(7 15 31 23)
//11
)
patch rod
//patch 4
(
(0 1 2 3)
//12
(16 19 18 17)
//13
mergePatchPairs
(
);
// ************************************************************** //
Appendix 2
Boundary type definition
SingleBarCase3D2
in
constant
/polyMesh
/boundary
| \\
/ F ield
| \\
/ O peration
| \\ /
A nd
M anipulation |
\\/
| Web:
www.OpenFOAM.org
|
|
\*---------------------------------------------------------------------------*/
FoamFile
{
version
2.0;
format
ascii;
class
polyBoundaryMesh;
location
"constant/polyMesh";
object
boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
6
(
airleft
{
type
patch;
nFaces
400;
startFace
45600;
}
airright
{
type
patch;
nFaces
400;
startFace
46000;
}
airdown
{
type
patch;
nFaces
400;
startFace
46400;
file
of
}
airup
{
type
patch;
nFaces
400;
startFace
46800;
}
rod
{
type
patch;
nFaces
800;
startFace
47200;
}
backandfront
{
type
patch;
nFaces
2400;
startFace
48000;
}
)
// ********************************************************************* //
Appendix 3
Initial field setting file of 0 /A in SingleBarCase3D1
/*---------------------------------------------------------------------------*\
| =========
| \\
/ F ield
| \\
/ O peration
| \\ /
A nd
M anipulation |
\\/
| Web:
http://www.openfoam.org
|
\*---------------------------------------------------------------------------*/
// Field Dictionary
FoamFile
{
version
2.0;
format
ascii;
class
volVectorField;
object
A;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions
[1 1 -2 0 0 -1 0];
boundaryField
{
airleft
{
type
fixedValue;
value
uniform (0 0 0);
airright
{
type
value
fixedValue;
uniform (0 0 0);
airdown
{
type
fixedValue;
value
uniform (0 0 0);
airup
{
type
value
fixedValue;
uniform (0 0 0);
rod
{
type
zeroGradient;
backandfront
{
type
zeroGradient;
}
// ********************************************************************* //
Appendix 4
Inner field setting file of system /setFieldsDict in SingleBarCase3D2
/*---------------------------------------------------------------------------*\
| =========
| \\
/ F ield
| \\
/ O peration
| \\ /
A nd
M anipulation |
\\/
| Web:
http://www.openfoam.org
|
\*---------------------------------------------------------------------------*/
FoamFile
{
version
2.0;
format
ascii;
class
dictionary;
object
setFieldsDict;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
defaultFieldValues
(
volVectorFieldValue J (0 0 0) // air box surrounding
);
regions
(
boxToCell
{
box (-0.21 -0.21 -0.01) (0.21 0.21 6.1);
fieldValues
(
volVectorFieldValue J (0 0 301500) // copper conductor
);
}
);
// ***************************************************************** //
Appendix 5
.C file of A-J solver
#include "fvCFD.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solve ( fvm::laplacian(A)==-muMag*J);
B= fvc::curl(A);
# include "IeEqn.H"
runTime++;
J.write();
A.write();
B.write();
return(0);
}
// ********************************************************************* //
References
1
Group CEDRAT, Flux 10 user guide. K101-10-EN-07/09, 2009, volume 3, p. 26104 and volume 4, p. 105.
Sonja Lundmark, Review of force and torque calculations from numerical field
solutions.
http://www.microwaves101.com/encyclopedia/highpermeability.cfm,High
Permeability Materials.
10