Analysis of Electric Field
Distribution Induced by
50Hz Magnetic Fields
Utilizing FastMultipole
SurfaceCharge
SimulationMethod for
VoxelModels
Shoji Hamada,
Tetsuo Kobayashi
(Kyoto Univ., Japan)
Main features
Objective
To develop a high speed, high capacity and high stability
Calculation method for induced ELF field in tissues.
(1) It requires O(D2)memory capacity and operation cost.
(2) It guarantees the solution globally satisfies Gauss’law.
Whole
region
D =
The number
of division
of a side
Superior to
O(D3) :FDM, FEM, FIT
O(D4 6
) : normal BEM,
SCM, MoM
Demonstration of the performance
Field analysis in a high resolution human model (8M voxels)
100 minutes on a Windows PC (2.66GHz, Core2duo)
We applied the diagonal form FMM to Surface
Chargesimulation Method (SCM) forvoxelmodels
The solution is globally stabilized.
Basic equations ( the effect of permittivity is neglected )
E: induced electric field (vector)
j: the imaginary unit
:angular frequency
A:vector potential equivalent to an applied magnetic field
:unknown scalar potential
=(1 j ) electric scalar potential
satisfiesLaplace equation
J: current density (vector)
: electrical conductivity (scalar)
Magnetically induced, low frequency, faint current in tissues
Boundary equations
The continuity of normal component of current density,
and the continuity of scalar potential
, ,
,
SurfaceChargesimulationMethod (SCM)
SCM = indirectBEM = a kind of Method of Moment
A numerical solver for Laplacian field
dS:small surface area
腅
d =
d( )=dE=
r
4 r
xdS
4 r2
xdS
r
r
By solving NBEQsfor Nunknown charge densities (x1
xN).
dS, x
x :apparent charge density on dS
,
=臧s d
= 臧s dE
SCM by definition satisfies . 膩don’t care
we can calculate and at arbitrary position.
SCM represents and as effects of surface charges
SCM only requires as BEQ.
xdS: a point charge
Inside
conductivity
Outside
conductivity
Voxelmodel Surface model
Surface
element
of SCM
Concept of SCM for voxelmodels
On each
element
Continuity of flux
through each element
=
Governing
linear system
Plane square shape
Constant surface charge density
Boundary matching Method of moment
(or a weighted residual method)
Element
Application of FMM to SCM for voxelmodels
m 芘m 芘m
Voxels
=
1 leaf cell
Leaf
cell
Sourcecharge element j( )
to flux through targetelement i( )
Octal tree
structure
for
the FMM
臞菓
i
= 臧Sj
Fluxi
=臧Si腩臧
Sj
腅nidSi
4 r2
xj
dSj
r
r
4 r2
xj
dSj
r
r
Leaf
Near region
= xj
= xj
We can
construct
whole tree
structure
Near part direct calculation procedure
(Analytical formula of square
charge field Is also available.)
dSi
dS芊
Index that indicates
relative position
of two elements
j
rel
(i,j)
Elem.
leafMn
m
leafLn
m
Flux
leafMn
m
leafLn
m
Sourcecharge element j( )
to Multipolecoefficients leafM
n
m
Local coefficients leafLn
m
to flux through targetelement i( )
FMM
leafMn
m
=臧Sj
TranslationQ M
xjdSj
= xj n
m
( j ) = xj
(
n
m)j
rel
Fluxi
=臧Si
腩TranslationL 腅nidSi
= ({Re/Im}leafLn
m)k k(i)
Far part indirect calculation procedure
= ({Re/Im}leafLn
m)k
( k)i
rel
Index that indicates relative
position in leaf cell
Index that indicates relative
position in leaf cell
We canpreliminarily calculate the proportional coefficients
j
rel, ( n
m)j
rel, and ( k)i
rel
by numerical integration.
It may require several minutes, but we have to do it once.
Elem. Flux
Elem. leafM
n
m
leafLn
m
Flux
negligibly small memory requirement
The required memory capacity for storing the coefficients
in arrays is only several Megabyte.
Using these coefficients skips numerical integration
during the iteration of the iterative solver.
by using a linear system solver, BiCGSTAB2, with FMM.
Now, we can solve the BEQs for charge densities x1
xN
fast speed performance
The interface with FMM is summarized as follows.
Representative field at the center of eachvoxel
we calculate six fluxes inside the voxelwall
and define Ecenter
as follows.
fx

fx
+
fz

fz
+
f
y

f
y
+
k
2
j
2
i
2
E
center
S
ff
S
ff
S
ff
zzyyxx
−
+
−
+
−
+
+
+
+
+
+
≡
Caused by
staircase shape
approximation
of voxelwall
Voxel
(Microscopically
natural behavior)
(Macroscopically
natural distribution)
Smooth
the field
near
boundaries
It represents induced field in good accuracy
except nearboundary region.
In nearboundary region
Sharp field fluctuation
One practical example of smoothing procedure
Sum of w6
for the voxelssatisfying both (a) and (b)
(a) treatedvoxelitself or
surrounding 3x3x31voxels
(b) the conductivity is the same
as that of treatedvoxel
(1) Weight value w6
( ) for each voxel
The number of ownvoxelsurfaces
not treated as boundary elements
large w6
= far from boundary = reliable
(2) Weight value w162
( ) for each voxel
(3) When w162 is neither 162 nor 0, we replace the Ecenter
by the following weighted average field.
: according to the rule (2)(a) and (2)(b)
w6=5
w6=4
w6=6
This is effective only
near the boundaries.
–0.100.1
–10
0
10
E
y
菊Vm
–1
xm
x
zuniform B
z
–0.100.1
0
1
2
3
E
y
菊Vm
–1
xm
x
z
magnetic
dipole
source M
z
Examples (after the near boundary smoothing)
Red lines: true values, black dots: calculated values.
Radius of
outer sphere
=0.1m
Voxelsize
=0.001m
Radius of
outer sphere
=0.1m
Voxelsize
=0.001m
‘Taro’model
Licensed by National Institute of
Information and communications
Technology (NICT)
Japanese adult male model
7,977,906 voxelsexcept air region
2mm 2mm 2mm cubic voxel
50 kinds of tissues
3,921,953 surface elements
Applied
50Hz
Homogeneous
1 T
B
x , B y
, B z
magnetic fields
Table 1. Number of tissue, conductivity (Sm
1), and voxels
Specifications of applied FMM
Maximum degree of the M, L expansions = 10
Maximum order of exponential expansion = 18
Leaf cell size = 5x5x5 voxels
Convergence criterion of the relative residual norm = 10 7
Regularization of the linear system: total charge = 0
0204060
10−6
10−4
10−2
100
Iterations
Relative residual norms
10
−1
10
0
10
1
Time/sec
Eq. (4)
Eq. (5)
Eq. (6)
M to M
M to L
L to L
Convergence characteristics
of the BiCGSTAB2
Time per FMM calculation
Elem.膨Flux
Elem.膨M
n
m
L
n
m
膨Flux
One iterationstep requires
two FMM calculations.
Bz
application
73
21.2s
14.4s
0.45s
0.69s
0.72s
0.45s
(37.9s)
16.7s
Average number of sources
Average M to L tarns. / cell =125
/ target for Eq. (4) =1650
B x (lefttoright)B
y (backtofront)B
z (foottohead)
100min. 15sec
74 iterations
Emax=147.0 Vm1
Eave= 11.1 Vm1
115min. 27sec
86 iterations
Emax=218.3 Vm1
Eave= 14.6 Vm1
99min. 07sec
73 iterations
Emax=171.1 Vm1
Eave= 9.79 Vm1
Emax
(zslice), Eave
(zslice), E
max
(tissue) and Eave
(tissue)
Smoothing significantly reduces maximum field
doesn’t change average field significantly
B x (lefttoright)B
y (backtofront)B
z (foottohead)
0
100
200
300
01000
0
10
20
30
0
100
200
300
0204060
0
10
20
30
zmm
Number of tissues
E
max
(z)菊Vm
−1
E
ave
(z)菊Vm
−1
E
max
(tissue) E
ave
(tissue)
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
0
100
200
300
01000
0
10
20
30
0
100
200
300
0204060
0
10
20
30
zmm
Number of tissues
E
max
(z)菊Vm
−1
E
ave
(z)菊Vm
−1
E
max
(tissue) E
ave
(tissue)
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After
smoothing
: Before
smoothing
: After
smoothing
0
100
200
300
01000
0
10
20
30
0
100
200
300
0204060
0
10
20
30
zmm
Number of tissues
E
max
(z)菊Vm
−1
E
ave
(z)菊Vm
−1
E
max
(tissue) E
ave
(tissue)
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
Conclusion
(4)The combination of
FMM, SCM, and voxelmodel
brings great advantages.
(3)Field smoothing procedure
near boundaries eliminates
sharp field fluctuations.
(2) Magnetically induced ELF
electric fields are successfully
analyzed in a human model
composed of 8 million voxels,
4 million surface elements.
(1)We applied the FMM to
an SCM for voxelmodels
to develop a high speed,
high capacity and high stability
field calculation method.
Givenvoxel shape
Q. Can we guess the original shape?
But sometimes we tend to prefer this one, which is
a kind of low pass filtered, that is smoothed, shape.
We also tend to prefer smoothed E & J near boundaries.
We should also be aware there is no physical reason.
腅腅腅
E臠
E臠0
Ans. We can’t.
f2
f1
f1
f1
f1
f0
f0: definite value
f1: values difficult to
integrate numerically
f2: easily integrated value
f1=0.25(f0f2)
f2
f1
f1
f1
f1
f0
f4
f3
f3
f3
f3
f2=4f3 +f
4
, f1=0.25(f04f
3f
4
)
better accuracy
Accuracy refinement technique
Voxel identification
X
Y
Z
ix
iy
iz
(ix,i
y,i
z)
Surface identification
X
Y
Z
ix
iy
iz
(ix,i
y, iz,1)
X
Y
Z
ix
iy
iz
(ix,i
y, iz,2)
X
Y
Z
ix
iy
iz
(ix,i
y, iz,3)
Identification of relative position of two surfaces
X
Z
Y
(ix
bix
a, iy
biy
a, iz
biz
a, is
b,i
s
a)
Xrel
Zrel
Yrel
(ix
a, iy
a
, iz
a, i
s
a)
(ix
b, iy
b, iz
b, is
b)
Voxel identification
X
Y
Z
ix
iy
iz
(ix,i
y,i
z)
Surface identification (alternative definition)
X
Y
Z
ix
iy
iz
(ix,i
y, iz,1)
X
Y
Z
ix
iy
iz
(ix,i
y, iz,2)
X
Y
Z
ix
iy
iz
(ix,i
y, iz,3)
Identification of relative position of two surfaces (alternativedefinition)
X
Z
(ix
a, iy
a
, iz
a, i
s
a)
Y
(ix
b, iy
b, iz
b, is
b)
(ix
bix
a, iy
biy
a, iz
biz
a, is
b,i
s
a)
Xrel
Zrel
Yrel
Examples ( : sharp field fluctuations near boundaries)
Red lines: true values, black dots: calculated values.
Radius of
outer sphere
=0.1m
Voxelsize
=0.001m
Radius of
outer sphere
=0.1m
Voxelsize
=0.001m
–0.100.1
–10
0
10
E
y
菊Vm
–1
xm
x
zuniform B
z
–0.100.1
0
1
2
3
E
y
菊Vm
–1
xm
x
z
magnetic
dipole
source M
z
The worst case of the field fluctuation near boundaries
(a)Five surfaces are used as boundary elements.
(b)The outside conductivity is zero (air).
0
00
0
0
0
Five fluxes are by definition zero.
The Gauss’lawset the sixth
flux to be zero.
The Ecenter
is set to zero.
It appears 100 % error.
A pit
This kind of problem is shared by numerical methods
based on voxelmodels & Gauss’law.
To understand the characteristics of each numerical method
inform us how to interpret and utilize the solution
in a sense of engineering.
B x
(left
to
right)
0
100
200
300
01000
0
10
20
30
0
100
200
300
0204060
0
10
20
30
zmm
Number of tissues
E
max
(z)菊Vm
−1
E
ave
(z)菊Vm
−1
E
max
(tissue) E
ave
(tissue)
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
0
100
200
300
01000
0
10
20
30
0
100
200
300
0204060
0
10
20
30
zmm
Number of tissues
E
max
(z)菊Vm
−1
E
ave
(z)菊Vm
−1
E
max
(tissue) E
ave
(tissue)
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After
smoothing
: Before
smoothing
: After
smoothing
B y
(back
to
front)
0
100
200
300
01000
0
10
20
30
0
100
200
300
0204060
0
10
20
30
zmm
Number of tissues
E
max
(z)菊Vm
−1
E
ave
(z)菊Vm
−1
E
max
(tissue) E
ave
(tissue)
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
: Before smoothing
: After smoothing
B z
(foot
to
head)
Comments 0
Log in to post a comment