Theoretical Formulation of a Time-Domain Finite Element Method for Nonlinear Magnetic Problems in Three Dimensions.

1. INTRODUCTION

Ferromagnetic compounds are widely used in data storage and processing, electric power generation, telecommunications, and transducers in smart systems. Magnetic hysteresis and constitutive nonlinearities are fundamental properties of all ferromagnetic materials, which exhibit nonlinear behavior where the permeability/susceptibility is a function of the magnetic field intensity. To achieve optimal material or device performance, the nonlinear hysteresis property of ferromagnetic materials must be modeled properly in the design process.

To describe the hysteresis property of ferromagnetic materials, many models have been proposed by physicists and material scientists in the past [1-3], among which the Preisach model [4] and the Jiles-Atherton (J-A) model [5,6] are very well known and commonly used. Different variations of these models have also been developed for different applications [7-10], including the generalizations from the static to the dynamic model, from the scalar to the vector model, and from the forward to the inverse model.

Even though the constitutive relations of ferromagnetic compounds have been well established, numerical modeling to predict the performance of a device made of such materials is still highly challenging because the overall performance of such a device is determined not only by the material it is made of, but also by many other factors such as the structure of the device, the external excitation, the interactions between different parts of the device and between the device and the environment it is installed in. To correctly model the device and generate reasonable prediction results, numerical methods have to combine the solution of Maxwell's equations with advanced hysteresis models. The capability of numerically simulating complicated structures made of nonlinear materials would provide an important analysis tool for a deeper insight into the physical and engineering behaviors.

During the past decades, extensive research was conducted and published on the modeling of the magnetic nonlinearities. For example, the nonlinear magnetic problems are analyzed using the finite element method in one dimension [11-13] and two dimensions [14,15], with the use of scalar hysteresis models. The vector hysteresis model is combined with a two-dimensional finite element analysis in [16] for the analysis of a T-joint region problem with a rotating flux excitation. In [17], a 2D vector hysteresis model is presented and included in a 2D eddy current finite element analysis. However, these research efforts have only dealt with either scalar hysteresis models or low-dimensional problems (1D or 2D) in magneto-static or magneto-quasi-static cases.

In order to model a real device, a vector hysteresis modeling combined with a three-dimensional dynamic finite element analysis is needed. To this end, a three-dimensional vector finite element method [18] in the frequency domain is developed to deal with nonlinear materials through the Newton-Raphson method with relaxation factors [19,20]. The method is further extended by using a domain decomposition method, the dual-primal finite element tearing and interconnecting (FETI-DP) method [21-23], to enhance its computational capability [24]. In the meantime, a preliminary study of the three-dimensional vector finite element method in the time domain is also conducted to deal with nonlinear materials [25].

The time-domain finite element analysis has two obvious advantages over the frequency-domain method. The first advantage comes from the convergence issue of the nonlinear solver. In a real application, there is no straightforward way to choose a good initial guess in the frequency domain that guarantees the convergence of the nonlinear solver. As a result, the nonlinear method could be unstable or divergent in solving complicated problems in the frequency domain. In a time- domain nonlinear solver, the numerical solution at a previous time step can be used as a natural choice for the initial guess at the successive time step. The convergence of the nonlinear method can always be achieved if the time-step size is small enough such that the solutions at two successive time steps are sufficiently close. The second advantage of the time-domain solver is the ability to capture all the physical phenomena in a real process, such as higher-order harmonic responses, magnetic remanence, and hysteresis losses. The understanding of these physical processes is beneficial, sometimes critical, to engineering design. However, there are also several major challenges in the time-domain modeling of nonlinear magnetic problems. One challenge comes from the fact that nonlinear magnetic devices often operate at very low frequencies. In such a frequency regime, the well-known low frequency breakdown problem becomes a big issue to the finite element analysis, which limits its solution accuracy and reduces its computational efficiency. Another challenge comes from the application of the widely used nonlinear method, the Newton-Raphson method, where the Jacobian matrix needs to be updated and solved at every Newton (nonlinear) iteration, which makes the solution very time consuming for large problems.

In this paper, the three-dimensional time-domain finite element method (TDFEM) [18] is employed and combined with the Newton-Raphson and Newton-like methods to solve magnetic- dynamic problems with nonlinear ferromagnetic materials modeled by the inverse vector J-A model [9,10] at relatively low frequencies. The second-order partial differential equation (PDE) in the time domain is first formulated, and the nonlinear hysteresis model is incorporated. The resulting nonlinear PDE is then discretized by employing the Newmark-[beta] scheme, and solved by applying the Newton-Raphson method. Different approaches of constructing the Newton-Raphson schemes are investigated and compared. The performances of the proposed methods are investigated and demonstrated through a number of numerical examples such as the simulation of a physical demagnetization process, the prediction of magnetic remanence in a ferromagnetic material, and the generation of higher- order harmonics.

2. FERROMAGNETICS AND HYSTERESIS MODELS

As one of the most widely used hysteresis models, the Jiles-Atherton (J-A) model [5,6] characterizes magnetic hysteresis through the reversible and irreversible dipole switching mechanisms and energy losses (domain wall losses) relative to the equilibrium anhysteretic magnetization. The construction of the J-A model thus includes the characterization of the anhysteretic magnetization [M.sub.an], the quantification of the irreversible magnetization [M.sub.irr], and the characterization of the reversible magnetization [M.sub.rev]. The total magnetization is then defined as M = [M.sub.irr] + [M.sub.rev].

2.1. Scalar Jiles-Atherton Model

In the J-A model, the magnetic hysteresis is physically attributed to the defects of material which cause a resistive force acting opposed to the domain wall motion [5]. In the ideal case, where the material is lossless, the magnetization M is a single-valued function of the magnetic field H. Such an magnetization is called anhysteretic magnetization [M.sub.an], and can be described by either the Ising relation [26]

[M.sub.an](H) = [M.sub.s] tanh ([H.sub.e]/a), (1)

or the Langevin relation [27,28]

[M.sub.an](H) = [M.sub.s][coth ([H.sub.e]/a) - a/[H.sub.e]]. (2)

In the above expressions, [M.sub.s] is the saturation magnetization of the material, a(T) = [H.sub.h]T /[[mu].sub.0][T.sub.c], where T is the temperature in Kelvin, [T.sub.c] is the Curie point of the ferromagnetic compound, [H.sub.h] is the bias magnetic field, and [H.sub.e] is the effective field experienced by each individual magnetic dipole moment, which is defined as

[H.sub.e] = H + [alpha]M, (3)

where [alpha] = [H.sub.h]/[[mu].sub.0][M.sub.s] is a material dependent constant.

In a real ferromagnetic compound, where loss is present, a change in [M.sub.irr] will introduce energy loss proportional to the magnitude of the change

d[epsilon] = [k.sub.p] [absolute value of d[M.sub.irr]], (4)

where [k.sub.p] is the irreversible loss constant in the domain wall models, which is related to the density of the pinning defects in the material. As a result, the magnetization energy increase comes from the total energy increase in the lossless case, minus the energy loss

[mathematical expression not reproducible], (5)

where

[mathematical expression not reproducible].

The irreversible magnetization can therefore be formulated as

d[M.sub.irr]/d[H.sub.e] = [M.sub.an] - [M.sub.irr]/[k.sub.p][delta]. (6)

However, in the case when [M.sub.an] > [M.sub.irr] and dH/dt < 0, Equation (6) leads to nonphysical negative slope d[M.sub.irr]/d[H.sub.e] < 0, which contradicts experimental observations. To eliminate such a discrepancy, it is suggested [29] to assume that when ([M.sub.an] - [M.sub.irr]) d[H.sub.e] < 0, there is no domain wall displacement and d[M.sub.irr] = 0. Defining

[mathematical expression not reproducible] (7)

we have

[d.sub.Mirr] = 1/[k.sub.p][delta] [[([M.sub.an] - [M.sub.irr]) d[H.sub.e]].sup.+]. (8)

The reversible magnetization [M.sub.rev] is proportional to the difference between [M.sub.an] and [M.sub.irr], which leads to the differential relation that resembles Hooke's law

d[M.sub.rev] = c(d[M.sub.an] - d[M.sub.irr]), (9)

where c is the Kelvin-Voigt damping coefficient, which quantifies the degree of reversibility. The total magnetization can then be determined as

[mathematical expression not reproducible]. (10)

Noting that

(1 - c)([M.sub.an] - [M.sub.irr]) = [M.sub.an] - [[M.sub.irr] + c([M.sub.an] - [M.sub.irr])] = [M.sub.an] - [[M.sub.irr] + [M.sub.rev]] = [M.sub.an] - M (11)

and defining [zeta] = [partial derivative][M.sub.an]/[partial derivative][H.sub.e], we have

dM = 1/[k.sub.p][delta] [[([M.sub.an] - M) d[H.sub.e]].sup.+] + c [zeta] d[H.sub.e]. (12)

Since [M.sub.an] is a function of [H.sub.e], thus a function of H, this yields

dM/dH = [[chi].sup.d](M,H, dH/[absolute value of dH]), (13)

where [[chi].sup.d] is the differential susceptibility. It can be found, after some mathematical manipulations, that

[[chi].sup.d] = [([M.sub.an] - M) /[k.sub.p][delta] + c [zeta]]/[1 - [alpha] [([M.sub.an] - M)/[k.sub.p][delta] + c [zeta]]] (14)

when ([M.sub.an] - M) d[H.sub.e] > 0, and

[[chi].sup.d] = c [zeta]/[1 - [alpha] c[zeta]] (15)

when ([M.sub.an] - M) d[H.sub.e] [less than or equal to] 0. From these differential relations, the magnetization M can be determined by the value and variation of the magnetic field H.

2.2. Vector Jiles-Atherton Model

The scalar J-A model is able to describe isotropic materials, and can be used in analyzing scalar problems. In order to characterize anisotropic materials and to be applied in the analysis of three-dimensional vector problems, a vector model is needed. In this work, the vector generalization of the J-A model developed in [9] is adopted.

By following a similar physical reasoning, the scalar variables used in the preceding section can be elevated to vectors and tensors. Specifically, the anhysteretic magnetization can be rewritten in vector form as

[M.sub.an](H) = [M.sub.an]([parallel][H.sub.e][parallel]) [H.sub.e]/[parallel][H.sub.e][parallel], (16)

in which [M.sub.an] can be described by different relations as in (1) and (2), and [H.sub.e] = H + [bar.[alpha]] x M, where [bar.[alpha]] is a tensor elevation of [alpha].

To obtain a model for the vector irreversible magnetization, an auxiliary vector can be defined to resemble the right-hand side of Eq. (6) as

[[chi]'.sub.f] = [[bar.k].sup.-1.sub.p] x ([M.sub.an] - [M.sub.irr]), (17)

where [[bar.k].sub.p] is the tensor correspondence of [k.sub.p]. In the anisotropic case, [[bar.k].sub.p] is a 3-by-3 matrix representing different loss constants in different directions, while in the isotropic case, it reduces to the scalar [k.sub.p]. If the direction of d[M.sub.irr] is assumed to be parallel to that of [[chi]'.sub.f], and its magnitude can be similarly determined by [[chi]'.sub.f] x d[H.sub.e], the change of the irreversible magnetization can be expressed as

d[M.sub.irr] = [[chi]'.sub.f] [[parallel][[chi]'.sub.f][parallel].sup.-1] [([[chi]'.sub.f] x d[H.sub.e]).sup.+] (18)

after taking into account the assumption that d[M.sub.irr] = 0 when [[chi]'.sub.f] x d[H.sub.e] < 0.

To model the vector reversible magnetization, a direct elevation of Eq. (9) can be used

d[M.sub.rev] = [bar.c] x (d[M.sub.an] - d[M.sub.irr]), (19)

where [bar.c] is the tensor correspondence of c. Finally, the total vector magnetization can be expressed as

[mathematical expression not reproducible], (20)

where [[chi].sub.f] = [[bar.k].sup.-1.sub.p] x ([M.sub.an] - M), and [bar.[zeta]] = [partial derivative][M.sub.an]/[partial derivative][H.sub.e].

Since [M.sub.an] is a function of [H.sub.e], thus a function of H, this yields

dM = [[bar.[chi]].sup.d](M, H, dH/[parallel]dH[parallel]) x dH, (21)

where [[bar.[chi]].sup.d] is the differential susceptibility tensor. It can be found, after some mathematical manipulations, that

[[bar.[chi]].sup.d] = [[[bar.I] - [[chi].sub.f][[parallel][[chi].sub.f][parallel].sup.-1][[chi].sub.f] x [alpha] - [bar.c] x [bar.[zeta]] x [bar.[alpha]].sup.-1] x [[[chi].sub.f][[parallel][[chi].sub.f][parallel].sup.-1] [[chi].sub.f] + [bar.c] x [bar.[zeta]]] (22)

when [[chi].sub.f] x d[H.sub.e] > 0, and

[[bar.[chi]].sup.d] = [[[bar.I] - [bar.c] x [bar.[zeta]] x [bar.[alpha]].sup.-1] x [[bar.c] x [bar.[zeta]]] (23)

when [[chi].sub.f] x d[H.sub.e] [less than or equal to] 0. From these differential relations, the magnetization M can be determined by the value and variation of the magnetic field H.

2.3. Inverse Vector Jiles-Atherton Model

In the finite element modeling of magnetic problems, the magnetic flux density B is usually employed as the unknown quantity. Therefore, it is more convenient to model the magnetization M and the magnetic field H as functions of the magnetic flux B. Since they are related by

B = [[mu].sub.0] (H + M), (24)

the hysteresis model (21) can be modified to relate the differential magnetization dM and the differential magnetic field dH to the differential magnetic flux dB by the inverse J-A model [10] as

dM = [[bar.[xi]].sup.d.sub.r](B, H) x dB/[[mu].sub.0] (25)

dH = [[bar.[upsilon]].sup.d.sub.r](B, H) x dB/[[mu].sub.0], (26)

where [[bar.[xi]].sup.d.sub.r] is the relative differential magnetizability tensor and [[bar.[upsilon]].sup.d.sub.r] is the relative differential reluctivity tensor. It can be found, after some mathematical manipulations, that

[mathematical expression not reproducible] (27)

when [[chi].sub.f] x d[H.sub.e] > 0, and

[[bar.[xi]].sup.d.sub.r] = [[bar.I] + [bar.c] x [bar.[zeta]] x [([bar.I] - [bar.[alpha]])].sup.-1] x [[bar.c] x [bar.[zeta]]] (28)

when [[chi].sub.f] x d[H.sub.e] [less than or equal to] 0. The differential reluctivity tensor can be obtained by

[[bar.[upsilon]].sup.d.sub.r] = [bar.I] - [[bar.[xi]].sup.d.sub.r]. (29)

From these differential relations, the magnetization M and the magnetic field H can be determined by the value and variation of the magnetic flux B.

3. NONLINEAR FINITE ELEMENT FORMULATION IN THE TIME DOMAIN

In this section, the TDFEM formulation for linear magnetic problems will first be derived, based on which the TDFEM formulations for nonlinear magnetic problems will then be constructed and discussed. The construction of the nonlinear TDFEM formulation will mainly be based on the most widely used numerical method for solving nonlinear equations, the Newton-Raphson method [30], due to its quadratic local convergence. Different interpretations to the Newton- Raphson method will lead to different formulation ideas. The resulting formulations will be compared and discussed. Moreover, the difference between the Newton-Raphson based nonlinear TDFEM formulations using the scalar and the vector hysteresis models will be investigated. At the end of this section, a simple fixed point method for solving the nonlinear TDFEM equation will be presented.

3.1. Time-Domain Finite Element Formulation for Linear Problems

To construct a TDFEM formulation for linear magnetic problems, we start from Maxwell's equations in the time domain

[nabla] x E = - [partial derivative]/[partial derivative]t B (30)

[nabla] x H = [partial derivative]/[partial derivative]t D + [sigma]E + [J.sub.imp]. (31)

By using the magnetic vector potential

B = [nabla] x A (32)

E = - [partial derivative]/[partial derivative]t A, (33)

the second-order curl-curl equation for the magnetic vector potential

[mathematical expression not reproducible] (34)

can be obtained, where [[epsilon].sub.r] and [[mu].sub.r] are the relative permittivity and permeability of the material, and [[eta].sub.0] and [c.sub.0] are the intrinsic impedance and the speed of light in free-space, respectively. Also, [J.sub.imp] is the impressed current source, which is usually solenoidal because it is generated by coils in magnetic problems. As a result, [J.sub.imp] can be expressed in terms of the impressed current vector potential [T.sub.imp] [31]

[J.sub.imp] = [nabla] x [T.sub.imp], (35)

which is equivalent to enforcing a gauge condition on A.

To discretize the curl-curl vector Equation (34), we first expand the magnetic vector potential in terms of basis functions [N.sub.j], which can be either low-order edge elements [32,33] or higher-order hierarchical basis functions [34]

A = [N.summation over (j=1)][a.sub.j](t) [N.sub.j](r). (36)

The semi-discrete system can be obtained as

[mathematical expression not reproducible], (37)

where

[mathematical expression not reproducible] (38)

[mathematical expression not reproducible] (39)

[mathematical expression not reproducible] (40)

[mathematical expression not reproducible] (41)

[mathematical expression not reproducible] (42)

In the above expressions, V is the solution domain, S = [partial derivative]V is the boundary of V, and [[eta].sub.r] is the relative impedance on the truncation boundary [35,36].

The semi-discrete Equation (37) is an ordinary differential equation (ODE), which can be solved using one of time integration methods. In this work, the well-known Newmark- [beta] method [37-39] is adopted because of its unconditional stability and second-order accuracy, when [beta] [greater than or equal to] 1/4. Specifically, in the Newmark-[beta] method, a time dependent variable is approximated by the weighted sum of its values at three consecutive time steps, and its time derivatives are approximated by central difference formulae

{u(t)} = [[beta]u.sup.n+1] + (1 - 2[beta]) [u.sup.n] + [[beta]u.sup.n-1] (43)

d/dt{u(t)} = 1/2[DELTA]t ([u.sup.n+1] - [u.sup.n-1]) (44)

[d.sup.2]/d[t.sup.2]{u (t)} = 1/[DELTA][t.sup.2] ([u.sup.n+1] - 2[u.sup.n] + [u.sup.n-1]). (45)

Application of the Newmark-[beta] method to the ODE in Eq. (37) yields the following fully discrete system (with [beta] = 1/4)

[mathematical expression not reproducible], (46)

which can be solved step by step from the known solutions at the (n - 1)th and the nth time steps to the (n + 1)th time step. Once the magnetic vector potential A is obtained at some time step [t.sub.n+1], the magnetic flux B and electric field E at the same time step can be calculated using Eqs. (32) and (33).

3.2. Newton-Raphson Method for Nonlinear Equations

For an N dimensional nonlinear problem

f(x) = 0, (47)

where f = [[[f.sub.1], [f.sub.2],...,[f.sub.N]].sup.T] contains N nonlinear functions, and x = [[[x.sub.1], [x.sub.2],...,[x.sub.N]].sup.T] is the N dimensional unknown vector, the Newton-Raphson method takes the following form

[x.sub.k+1] = [x.sub.k] - [[J ([x.sub.k])].sup.-1] f([x.sub.k]), (48)

where J is the Jacobian matrix of the nonlinear function f

[{J (x)}.sub.ij] = [partial derivative][f.sub.i](x)/[partial derivative][x.sub.j]. (49)

The nonlinear function f is also known as the residual since it vanishes at the true solution.

In order to solve the nonlinear Equation (47), an initial guess of solution [x.sub.0] is first given. At the kth nonlinear (Newton) step (k = 0, 1, 2,...), the Jacobian matrix equation, which is a set of linear equations, is solved for the incremental vector [s.sub.k]

J ([x.sub.k]) [s.sub.k] = -f ([x.sub.k]). (50)

The nonlinear solution at the (k + 1)th Newton step is updated by

[x.sub.k+1] = [x.sub.k] + [s.sub.k]. (51)

The Jacobian matrix equation is then solved again at the (k + 1)th step to find next update to the nonlinear solution until convergence is achieved. This procedure is known as the Newton iteration. Apparently, the Newton iteration Eq. (48) has a clear graphical meaning. Take a one-dimensional problem as an example, the meaning of the Jacobian matrix equation (50) is to draw a tangent line from the trial solution [x.sub.k], and find the next estimated solution [x.sub.k+1] from the intersection with the x axis. For a higher-dimensional problem, the Jacobian matrix simply stands for the higher-dimensional tangent surface.

Algebraically, the Newton-Raphson method can also be interpreted from the Taylor series expansion point of view. Expanding the nonlinear function f in the neighborhood of [x.sub.k] in terms of the Taylor series and keeping only the zeroth and first order terms, we have

f([x.sub.k+1]) = f([x.sub.k] + [s.sub.k]) = f([x.sub.k]) + [partial derivative]f/[partial derivative][x.sub.k] [s.sub.k] + O([s.sup.2.sub.k]). (52)

Since the objective is to make the nonlinear function f vanish, we have

f([x.sub.k]) + [partial derivative]f/[partial derivative][x.sub.k] ([x.sub.k+1] - [x.sub.k]) = 0, (53)

[x.sub.k+1] = [x.sub.k] - [([partial derivative]f/[partial derivative][x.sub.k]).sup.-1] f([x.sub.k]). (54)

Noting the definition of the Jacobian matrix equation (49), the above expression is the same as the Newton-Raphson method described in Eq. (48).

3.3. Direct Derivative Based Newton-Raphson Scheme for TDFEM

To construct a TDFEM formulation for nonlinear problems, we follow the Newton iteration defined in Eqs. (48)-(51), which requires the direct derivative of the nonlinear function with respect to the unknown vector. In this case, The nonlinear function, also known as the residual {r}, is expressed as

{r} = [[K].sup.n+1] [{a}.sup.n+1] + [[K].sup.n] [{a}.sup.n] + [[K].sup.n-1] [{a}.sup.n-1] - {[??]}, (55)

where

[mathematical expression not reproducible] (56)

[[K].sup.n] = 1/2 [c.sup.2.sub.0][DELTA][t.sup.2] [S] - 2[M] (57)

{[??]} = [[eta].sub.0][c.sub.0][DELTA][t.sup.2](1/4 [{b}.sup.n+1] + 1/2 [{b}.sup.n] + 1/4 [{b}.sup.n-1]). (58)

Since in the definition of the stiffness matrix [S]

[mathematical expression not reproducible], (59)

the relative reluctivity of the material [[upsilon].sub.r] = 1/[[mu].sub.r] is a function of the magnetic flux density, which is a function of the unknown vector [{a}.sup.n+1], the residual {r} in Eq. (55) is nonlinear.

3.3.1. Nonlinear TDFEM Formulation with a Scalar Constitutive Model

If the constitutive relation of the ferromagnetic material is described by a scalar relation, such as the one expressed in Eq. (13), the relative reluctivity is a function of B = [parallel]B[parallel]

[[upsilon].sub.r] = [[upsilon].sub.r](B). (60)

Since in the TDFEM formulation, the magnetic flux B is expressed by the magnetic vector potential A, which is expanded in terms of the basis functions in space and the Newmark- [beta] method in time

[mathematical expression not reproducible], (61)

from Eq. (55) it is clear that the partial derivative of [r.sub.i] with respect to [a.sup.n+1.sub.j] results in not only the finite element entry [K.sup.n+1.sub.ij], but also several terms involving the partial derivative of a stiffness matrix element [S.sub.ik] with respect to [a.sup.n+1.sub.j] as well. After some mathematical manipulations, the Jacobian matrix at (n + 1)th time step is obtained as

[mathematical expression not reproducible]. (62)

The partial derivative of the stiffness matrix can be further expressed as

[mathematical expression not reproducible], (63)

where

[partial derivative][[upsilon].sub.r](B)/[partial derivative][a.sup.n+1.sub.j] = [partial derivative][[upsilon].sub.r](B)/[partial derivative]B [partial derivative]B/[partial derivative][a.sup.n+1.sub.j] (64)

by using the chain rule of derivatives. In the above expression

[partial derivative]B/[partial derivative][a.sup.n+1.sub.j] = 1/4 ([nabla] x [N.sub.j]) x B/B (65)

[partial derivative][[upsilon].sub.r]/[partial derivative]B = [partial derivative]/[partial derivative]B ([[mu].sub.0]H/B) = 1/B ([[upsilon].sup.d.sub.r](B) - [[upsilon].sub.r](B)], (66)

where [[upsilon].sup.d.sub.r] is the relative differential reluctivity, and [[upsilon].sub.r] is the relative reluctivity. They can be obtained from the nonlinear constitutive relation of the material described in Section 2.1.

Substituting Eqs. (63)-(66) back into Eq. (62), we obtain

[mathematical expression not reproducible]. (67)

Interchanging the summation and the integration and noting

[mathematical expression not reproducible], (68)

Equation (67) is reduced to

[mathematical expression not reproducible]. (69)

Since the second term in Eq. (69) is related to the derivative of the stiffness matrix, we define the incremental stiffness matrix [[DELTA]S] as

[mathematical expression not reproducible] (70)

and the summation of [S] and [[DELTA]S] as

[[S.sup.J]] = [S] + [[DELTA]S]. (71)

The Jacobian matrix can be written as

[mathematical expression not reproducible], (72)

which is the final expression for the Jacobian matrix in the nonlinear TDFEM formulation if the scalar constitutive relation is used.

3.3.2. Nonlinear TDFEM Formulation with a Vector Constitutive Model

If the constitutive relation of the ferromagnetic material is described by a vector relation, such as those expressed in Eqs. (21), (25), and (26), the relative reluctivity in the stiffness matrix [S] becomes a tensor; hence,

[mathematical expression not reproducible]. (73)

The Jacobian matrix at the (n + 1)th time step is obtained similarly as in Eq. (62). The difference in the formulations between using the scalar and the vector constitutive relations occurs in the partial derivative of the stiffness matrix

[mathematical expression not reproducible]. (74)

Application of the chain rule of derivatives yields

[mathematical expression not reproducible]. (75)

In the above expression

[partial derivative]B/[partial derivative][a.sup.n+1.sub.j] = 1/4 [nabla] x [N.sup.j]. (76)

The Jacobian matrix therefore becomes

[mathematical expression not reproducible]. (77)

Unfortunately, [partial derivative][[bar.[upsilon]].sub.r]/[partial derivative]B cannot be obtained in a straightforward manner. To achieve its explicit expression, on one hand we expand H in terms of Taylor series (ignore the higher-order terms)

H([B.sub.2]) = H([B.sub.1] + [DELTA]B) = H([B.sub.1]) + [partial derivative]H/[partial derivative][B.sub.1] x [DELTA]B = H([B.sub.1]) + [[bar.[upsilon]].sup.d]([B.sub.1]) x [DELTA]B, (78)

where [[bar.[upsilon]].sup.d] = [[bar.[upsilon]].sup.d.sub.r]/[[mu].sub.0] = [[bar.[upsilon]].sup.d.sub.r][[upsilon].sub.0] is the differential reluctivity tensor. On the other hand, we can also express H as

[mathematical expression not reproducible], (79)

where [bar.[upsilon]] = [[bar.[upsilon]].sub.r]/[[mu].sub.0] = [[bar.[upsilon]].sub.r][[upsilon].sub.0] is the reluctivity tensor. Comparing Eqs. (78) and (79), it is clear that

[[bar.[upsilon]].sup.d.sub.r]([B.sub.1]) x [DELTA]B = [[bar.[upsilon]].sub.r]([B.sub.1]) x [DELTA]B + [partial derivative][[bar.[upsilon]].sub.r]/[partial derivative][B.sub.1] x [DELTA]B x [B.sub.1]. (80)

Noting that [DELTA]B is the incremental vector [s.sub.k] in the Jacobian matrix equation (50), which is expanded as

[DELTA]B = [summation over (j)] [nabla] x [N.sub.j] [DELTA][a.sub.j], (81)

we have

[partial derivative][[bar.[upsilon]].sub.r]/[partial derivative]B x ([nabla] x [N.sub.j]) x B = ([[bar.[upsilon]].sup.d.sub.r] - [[bar.[upsilon]].sub.r]) x ([nabla] x [N.sub.j]). (82)

Substituting Eq. (82) back into Eq. (77), the incremental stiffness matrix [[DELTA]S] can be written as

[mathematical expression not reproducible] (83)

and

[mathematical expression not reproducible]. (84)

The Jacobian matrix has the same form as that in Eq. (72). The relative differential reluctivity [[bar.[upsilon]].sup.d.sub.r] can be obtained from the nonlinear constitutive relation of the material described in Section 2.2.

3.3.3. Remarks

After the Newton-Raphson scheme of the TDFEM is constructed, the nonlinear equation can be solved iteratively. Specifically, there are several levels of iteration the solution process goes through in order to complete the simulation. The top level of iteration is the marching-on-in- time process. At each time step, the excitation from an external source, in this case [J.sub.imp] or [T.sub.imp], is specified, which results in a nonlinear equation given by Eq. (46). To solve such a nonlinear equation, the Newton-Raphson method is used, which leads to the second level of iteration, the Newton iteration. At each Newton step, an initial guess of solution is given, usually chosen the same as the solution from the preceding time step. From the initial guess, the nonlinear residual function can be calculated using Eq. (55), and the Jacobian matrix can be assembled using either Eq. (72) or Eq. (84). The Jacobian matrix equation is then solved for the incremental vector using either a direct or an iterative linear solver (the third level of iteration), and the estimated solution is updated using Eq. (51). The Newton iteration terminates upon the convergence of the estimated solution, and the solution process proceeds to the next time step.

It can be seen from the above description that the total computational cost comes from the multiplication of three parts, which are the total number of time steps, the average number of Newton iterations at each time step, and the computational cost for solving the linear problem (the Jacobian matrix equation). The most time-consuming part of the method is the construction and the solution of the Jacobian matrix equation. To construct the Jacobian matrix equation, the Jacobian matrix has to be assembled and solved at each Newton step, which involves one matrix assembly and one matrix solution. Besides the Jacobian matrix, the nonlinear residual function has to be updated as well. In the direct derivative based Newton-Raphson scheme, in order to update the nonlinear residual function, the stiffness matrix Eq. (59) has to be assembled repeatedly because the value of [[upsilon].sub.r] is constantly changing. The solution of the Jacobian matrix equation is even more computationally expensive. Since the Jacobian matrix is being updated at every Newton iteration, to solve it with a direct solver such as the LU decomposition [30], the matrix needs to be factorized repeatedly. Similarly, to solve the Jacobian matrix equation with an iterative solver such as GMRes [40] or BiCGstab [41], a good preconditioner such as the incomplete LU (ILU) [42] or the sparse approximate inverse (SAI) [43] preconditioner needs to be constructed repeatedly based on the latest Jacobian matrix.

Another comment needs to be made is that if the Jacobian matrix is constructed accurately, and the initial guess is chosen sufficiently close to the true solution, the convergence of the Newton-Raphson method is quadratic. However, it should also be noted that the objective here is to make the nonlinear residual function (55) vanish. Therefore, even if the Jacobian matrix is not constructed precisely, as long as it contains the information of the tangent plane, the nonlinear problem can still converge, but likely with a slower rate. This argument paves the way for other nonlinear solvers that use the approximation of the Jacobian matrix, such as Broyden's method, which will be discussed in a later paper.

3.4. Polarization Technique Based Newton-Raphson Scheme for TDFEM

Another approach of constructing the Newton-Raphson scheme for the TDFEM is to use the idea of Taylor series expansion Eqs. (52)-(54). To this end, the polarization method [44] is first adopted to separate the magnetic field into the linear and the nonlinear parts

[mathematical expression not reproducible], (85)

where [[upsilon].sup.opt.sub.r] is a constant, which can be chosen as ([[upsilon].sup.min.sub.r] + [[upsilon].sup.max.sub.r])/2 for instance. In this case, the curl-curl equation becomes

[mathematical expression not reproducible]. (86)

By treating R as an unknown quantity and going through the same discretization process as the one described in Section 3.1, the following nonlinear matrix equation can be obtained

[[[K.sub.p]].sup.n+1] [{a}.sup.n+1] + [[[K.sub.p]].sup.n] [{a}.sup.n] + [[[K.sub.p]].sup.n-1] [{a}.sup.n-1] + {R} = {[??]}, (87)

where [[[K.sub.p]].sup.n+1], [[[K.sub.p]].sup.n], and [[[K.sub.p]].sup.n-1] have the same form as those in Eqs. (56)- (57), with the difference being the definition of the stiffness matrix [[S.sub.p]]

[mathematical expression not reproducible]. (88)

Note that the [[K.sub.p]] matrices are now no longer functions of the unknown vector [{a}.sup.n+1]. The nonlinearity occurs in

[mathematical expression not reproducible], (89)

where R is a nonlinear function of the vector potential A. The temporal and spatial discretization yields

[mathematical expression not reproducible], (90)

and the nonlinear residual function can be expressed as

{r} = [[[K.sub.p]].sup.n+1] [{a}.sup.n+1] + [[[K.sub.p]].sup.n] [{a}.sup.n] + [[[K.sub.p]].sup.n-1] [{a}.sup.n-1] + {R}-{[??]}. (91)

3.4.1. Nonlinear TDFEM Formulation with a Scalar Constitutive Model

If the constitutive relation of the ferromagnetic material is described by a scalar relation [upsilon]r = [upsilon]r (B), the nonlinear part of the magnetic field can be written as

R = [[[upsilon].sub.r](B) - [[upsilon].sup.opt.sub.r][[upsilon].sub.0]B. (92)

The Taylor series expansion of Eq. (91) requires the Taylor series expansion of the nonlinear quantity R

[mathematical expression not reproducible], (93)

where

[mathematical expression not reproducible]. (94)

From Eqs. (65), (66), and (76), it can be obtained that

[mathematical expression not reproducible]. (95)

By substituting Eq. (95) back into Eq. (93), and then into Eq. (89), the Taylor series expansion of the nonlinear residual function {r} can be obtained as

[mathematical expression not reproducible]. (96)

By enforcing the left-hand side of the above expression to be zero, the Newton- Raphson scheme can be derived as

[J] {[delta]a} = -{r}, (97)

where

[mathematical expression not reproducible] (98)

[J.sub.ij] = [K.sup.n+1.sub.p,ij] + [partial derivative][R.sub.i]/[partial derivative][a.sup.n+1.sub.j] (99)

and

[mathematical expression not reproducible]. (100)

Using the definition of the stiffness matrix [S] in Eq. (59), the definition of the stiffness matrix [[S.sub.p]] in Eq. (88), and the definition of the incremental stiffness matrix [[DELTA]S] in Eq. (70)

[partial derivative][R.sub.i]/[partial derivative][a.sup.n+1.sub.j] = 1/4[c.sup.2.sub.0][DELTA][t.sup.2]([S.sub.ij] - [S.sub.p,ij] + [DELTA][S.sub.ij]), (101)

we have

[mathematical expression not reproducible], (102)

which is the same as Eq. (72). However, it should be noted that although the Jacobian matrix is the same as that in Section 3.3.1, the residual vectors {r} in these two methods are different, which will result in a different computational efficiency. This will be elaborated in Section 3.4.3.

3.4.2. Nonlinear TDFEM Formulation with a Vector Constitutive Model

If the constitutive relation of the ferromagnetic material is described by a vector relation, the nonlinear part of the magnetic field has to be written as

R = [[[bar.[upsilon]].sub.r](B) - [[upsilon].sup.opt.sub.r][bar.I]] x [[upsilon].sub.0]B, (103)

where [bar.I] is the unit tensor. In the Taylor series expansion of Eq. (93),

[mathematical expression not reproducible]. (104)

The application of Eq. (82) yields

[mathematical expression not reproducible]. (105)

Following the same process as described in the preceding section, we have

[mathematical expression not reproducible] (106)

and hence,

[mathematical expression not reproducible], (107)

which is the same as that in Section 3.3.2. However, as mentioned in the previous section, the residual vectors {r} in these two methods are expressed differently, which will result in a different computational efficiency. This will be elaborated in Section 3.4.3.

Another approach to deriving the Jacobian matrix is to expand the nonlinear term R directly into Taylor series

[mathematical expression not reproducible]. (108)

Since B = [nabla] x A, we have [DELTA]B = [nabla] x [DELTA]A; therefore,

R(A + [DELTA]A) [approximately equal to] R(A) + ([[bar.[upsilon]].sup.d.sub.r] - [[upsilon].sup.opt.sub.r][bar.I]) x [[upsilon].sub.0] [nabla] x [DELTA]A. (109)

With the same small incremental field [DELTA]A, the magnetic field can be expressed as

[mathematical expression not reproducible]. (110)

As a result,

[nabla] x H(A + [DELTA]A) = [[upsilon].sub.0][nabla] x ([[bar.[upsilon]].sup.d.sub.r] [nabla] x [DELTA]A) + [[upsilon].sub.0][nabla] x ([[upsilon].sup.opt.sub.r] [nabla] x A) + [nabla] x R, (111)

which can be substituted back into the curl-curl Equation (34) to obtain the Newton-Raphson scheme.

3.4.3. Remarks

The solution process of the nonlinear TDFEM is the same as that described in Section 3.3.3. The difference here, as mentioned before, is the way by which the residual vector {r} is evaluated. In the direct derivative based Newton-Raphson scheme, the residual is evaluated as

{r} = [[K].sup.n+1] [{a}.sup.n+1] + [K[].sup.n] [{a}.sup.n] + [[K].sup.n-1] [{a}.sup.n-1] - {[??]}, (112)

whereas in the polarization technique based Newton-Raphson scheme, it is evaluated as

{r} = [[[K.sub.p]].sup.n+1] [{a}.sup.n+1] + [[[K.sub.p]].sup.n] [{a}.sup.n] + [[[K.sub.p]].sup.n-1] [{a}.sup.n-1] + {R} - {[??]}. (113)

Apparently, in order to evaluate {r} using Eq. (112), the [K] matrices need to be updated, which involve the evaluation of the stiffness matrix Eq. (59) or Eq. (73). The matrix-vector products [[K].sup.n+1] [{a}.sup.n+1], [[K].sup.n] [{a}.sup.n], and [[K].sup.n-1] [{a}.sup.n-1] are then performed and {r} is finally calculated. At every Newton step in the Newton iteration process, this evaluation usually has to be done more than once, if the relaxation method or the safe-guarding method is adopted to accelerate the convergence of the Newton-Raphson method. In contrast, to evaluate {r} using Eq. (113), only the vector {R} needs to be updated and one matrix-vector product [[[K.sub.p]].sup.n+1] [{a}.sup.n+1] needs to be performed because the stiffness matrix [[S.sub.p]] Eq. (88) involved in the expression is not a function of the unknown vector [{a}.sup.n+1]. As a result, the polarization technique based method is expected to have a better computational efficiency than the direct derivative based method. If the residuals expressed in Eqs. (112) and (113) are examined more carefully, it can be seen that since

[mathematical expression not reproducible], (114)

the stiffness matrix related terms in Eq. (112) can be rewritten as

[mathematical expression not reproducible], (115)

which are those in Eq. (113). From the above derivation, it is obvious that although the residual vectors {r} are evaluated in different ways in both formulations, they are actually the same. Since the Jacobian matrices in both formulations are also the same, the Newton-Raphson schemes for the TDFEM formulations introduced in Sections 3.3 and 3.4 are identical. Furthermore, it is also clear that the value of [[upsilon].sup.opt.sub.r] does not affect the final formulation, and therefore, it can be safely chosen to be 0 in the Newton-Raphson formulation.

3.5. Polarization Technique Based Fixed-Point Scheme for TDFEM

A simple fixed-point method can be constructed based on the polarization technique. Applying Eq. (85), Maxwell's Equation (31) becomes

[mathematical expression not reproducible], (116)

where A can be solved with a fixed-point iteration

[mathematical expression not reproducible]. (117)

The above formulation is almost the same as that in Eq. (86), and the only difference is that the nonlinear term -[[mu].sub.0][nabla] x R is now moved to the right-hand side. Such a difference will result in a different solution strategy. Specifically, to solve Eq. (117) with a fixed-point iteration, the following steps are performed.

(i) At time step n, set up the initial guess [A.sup.n.sub.0] = [A.sup.n-1];

(ii) At the kth fixed-point iteration (k = 0, 1, 2,...), calculate [B.sup.n.sub.k] = [nabla] x [A.sup.n.sub.k], and evaluate [H.sup.n.sub.k] using the constitutive model described in Section 2;

(iii) Calculate the nonlinear term with [R.sup.n.sub.k] = [H.sup.n.sub.k] - [[upsilon].sup.opt.sub.r] [[upsilon].sub.0] [B.sup.n.sub.k], and update right-hand side of the FEM system;

(iv) Solve the linear FEM system for [A.sup.n.sub.k+1];

(v) If converged, n := n + 1, go to step 1; otherwise, k := k + 1, go to step 2.

Similar to the Newton-Raphson scheme, the total computational cost of the fixed- point scheme also comes from the multiplication of three parts, which are the total number of time steps, the average number of fixed-point iterations at each time step, and the computational cost for solving the linear fixed-point problem. Since the system matrix resulting from Eq. (117) is not a function of the unknown vector, it only needs to be factorized once during the entire solution process. However, the convergence of the fixed-point iteration is approximately linear, which would result in a higher computational cost if strong nonlinearity is involved.

4. NUMERICAL EXAMPLES

In this section, several examples are given to demonstrate the capability of the proposed methods. The applications in the simulation of magnetic hysteresis and the prediction of higher-order harmonic generation are presented. The accurate and efficient numerical solution of the nonlinear formulation will be discussed separately in another paper. The models used in this section are adopted from several benchmark problems from the "testing electromagnetic analysis methods" (TEAM) workshop [45].

4.1. TEAM Problem 10

The TEAM workshop problem 10 is first considered. Shown in Fig. 1 is the problem setting, where a racetrack coil is placed between two steel channels and a flat steel plate is inserted in between of these two channels. The structure of the coil is shown in Fig. 1(a) with its geometrical sizes labeled. The geometry and dimension of the steel plate and channels are shown in Fig. 1(b). It should be noted that there is a very small air gap between the center plate and the channels, which is only 0.5 mm. Such an air gap would allow some leakage of magnetic flux, and therefore, in the discretization of the geometry, the gap should be discretized with sufficient spatial resolution in order for the leakage to be captured. The steel channels and the center plate are made of ferromagnetic material, which have a relative permittivity of [[epsilon].sub.r] = 1.0, a conductivity of [sigma] = 7.505 x [10.sup.6] S/m, and their permeability is given by the nonlinear B-H relation through measurement data as shown in Fig. 2(a), where the dots are the discrete data from measurement, and the curve is a continuous model obtained from curve fitting.

The racetrack coil has 162 turns in total. When it is turned on, the current flowing in the coil increases gradually from 0 A to a maximum value of [I.sub.m], which can be expressed mathematically as

[mathematical expression not reproducible] (118)

where [tau] is a parameter that controls how fast the current increases. The impressed current is shown in Fig. 2(b) as a function of time with [I.sub.m] = 5.64 A, and [tau] = 0.05 s.

The magnetic flux density B is recorded along line [L.sub.1] indicated in Fig. 3(a), as a function of position. The proposed nonlinear solver is used to advance time from 0 to that when a steady state is achieved. The problem is solved by setting [I.sub.m] = 5.64 A and [I.sub.m] = 18.52 A, respectively, and the results are presented in Fig. 3(b) and compared with those obtained by a static solver [24]. Excellent agreement is achieved. Next, the magnitude of the averaged magnetic flux in areas [S.sub.1], [S.sub.2], and [S.sub.3] are shown as a function of time in Figs. 3(c) and 3(d). The numerical results are obtained by using the Newton-Raphson method as well as the fixed-point method, as discussed in Section 3. Compared with the measurement results, good agreement is achieved.

4.2. Hysteresis Modeling

To investigate the modeling of the magnetic hysteresis using the inverse Jiles- Atherton vector model introduced in Section 2.3, we consider the magnetic problem involved in the demagnetization process and solve it with the Newton-Raphson scheme proposed in Section 3.4.

4.2.1. Demagnetization Process and Magnetic Remanence

Sometimes it is necessary to demagnetize or re-magnetize a magnetic device if its magnetic property has been distorted by heat, large magnetic fields, or stress. Failure in demagnetization would cause the remanent magnetic field, which may affect the normal operation of the device, or interfere with other electronic equipments. One way to remove the magnetization from a device is to apply an oscillating field with a reducing amplitude to it, causing the magnetic dipoles in the magnet to flip back and forth and eventually get neutralized.

To demonstrate the capability of the proposed method in solving magnetic hysteresis problems, the nonlinear magnetization and demagnetization process is simulated. The model is adopted from TEAM workshop problem 10 [46], with an excitation current designed as

[mathematical expression not reproducible], (119)

where [I.sub.m] = 2.5 A, [f.sub.0] = 10 Hz, [T.sub.m] = 0.4 s, and the coil is assumed to have 100 turns. The coil current is shown in Fig. 4(a), from which it can be seen that the current is turned on from 0 A, increasing to its maximum where the structure is magnetized. It is then oscillating at a frequency of 10 Hz, while the amplitude of its envelope is decreasing linearly. After the current fails back to 0 A at 0.4 s, it stays at 0 A for another 0.1 s for the magnetic response of the structure to be stabilized.

To model the ferromagnetic material, the inverse J-A model with its anhysteretic property described by the Langevin relation is used. The model parameters are chosen as [M.sub.s] = 1.3 x [10.sup.6] A/m, a = 25.3 A/m, [k.sub.p] = 66.6 A/m, [alpha] = 0, and c = 0.2. In Fig. 4(b), a typical hysteresis loop of a demagnetization process using this set of parameters is shown.

The magnetization M and magnetic field H at three different observation points ([P.sub.1], [P.sub.2], and [P.sub.3]) in Fig. 5(a) are recorded, and the demagnetization process and magnetic remanence are shown in Figs. 5(b)-5(d). From the figures, it can be seen that among these three sampling locations, only the steel at point [P.sub.1] located at the center of the center plate has been fully saturated by the externally applied current. The steel at the other two locations are not fully saturated. Clearly, the hysteresis histories and the magnetic remanence at different points are different due to the structure of the channel and the way it is excited. This example demonstrates successfully the capability and importance of the proposed method in the simulation and design of a ferromagnetic device. From this simulation and the results shown in Figs. 5(b)-5(d) and their inserts, it can be suggested that in order to obtain a better demagnetization result, it is better for the entire structure to be fully saturated first, then decrease the applied field with a slower rate such that more cycles are applied to the structure and smaller minor hysteresis loops can be formed, which eventually results in a smaller remanence.

4.2.2. Generation of Higher-order Harmonics

A well-known consequence of nonlinearity is the generation of higher-order harmonics in a dynamic problem. Here, we investigate such a phenomenon using the model adopted from TEAM workshop problem 32 [47], where a three-limbed ferromagnetic core is considered. As shown in Fig. 6, the core is made of five layers of 0.48-mm-thick Fe-Si 3.2%wt laminations, having a conductivity of [sigma] = 1.78 x [10.sup.6] S/m. Two 90-turn windings are placed on the external limbs to supply excitations. In this problem, the excitation current flowing in each winding is a sinusoidal function with a magnitude of 1.15 A and a frequency of 10 Hz. The Fe-Si material is modeled by the inverse J-A model with the Langevin anhysteretic relation. The model parameters are chosen as [M.sub.s] = 1.168 x [10.sup.6] A/m, a = 60 A/m, [k.sub.p] = 130 A/m, [alpha] = [10.sup.-4], and c = 0.2.

With the 10-Hz sinusoidal current excitation, the magnetic flux density is solved and shown in Fig. 7. From this figure, it can be seen clearly that the magnetic flux flowing in the core follows the right-hand rule. The field singularities can also be observed at the corners of the structure. The magnetization M and magnetic field H at two observation points are recorded. One is at the center of the center limb and the other one is at the center of the right limb. The hysteresis loops at these two locations are shown in Fig. 8, from which different hysteresis behavior can be observed. The magnetic flux B at these two observation locations are also recorded and shown in Fig. 9 as a function of time.

To demonstrate the generation of higher-order harmonics, the magnetic flux results in the time domain are converted into the frequency domain using the Fourier transform. The power spectrum of the excitation current and the corresponding magnetic flux are shown in Fig. 10. From these figures, it is very clear that the input current has a power concentrating at 10 Hz, as expected. The magnetic flux, as a nonlinear response, contains a power spectrum distributing not only at the fundamental frequency of 10 Hz, but also at higher-order harmonics of 30 Hz, 50 Hz, 70 Hz, etc., with a decreasing power magnitude, which indicates that the higher-order harmonics are captured accurately. This example demonstrates successfully the capability of the proposed nonlinear time domain solver in predicting the nonlinear behavior of a ferromagnetic device.

Note that for time-harmonic nonlinear problems, the harmonic balance method [48] can be employed to solve for higher-order harmonics with a frequency-domain solver. However, when the problem is not time harmonic, such a method cannot be applied because the response will have a continuous spectrum distribution. One example would be the demagnetization process discussed in the preceding section. As a result, the time-domain nonlinear solver is much more versatile.

5. CONCLUSION

In this work, a numerical solution to nonlinear magnetic problems was formulated with a full-wave time-domain finite element method in three dimensions. Nonlinear magnetic hysteresis models were first introduced. The vector curl-curl equation was then formulated in the time domain with the magnetic vector potential employed as the unknown function. Both Newton-Raphson and fixed-point methods were explored and the corresponding formulations were constructed. The capability of these formulations were demonstrated through several well-acknowledged numerical examples. The applications in the simulation of magnetic hysteresis and demagnetization process were shown, and the prediction of the harmonic generation was also demonstrated.

ACKNOWLEDGMENT

This work was supported in part by a grant from Sandia National Laboratories and in part by the National Natural Science Foundation of China (NSFC) under Contract No. 61201012 and the Fundamental Research Funds for the Central Universities under Contract No. ZYGX2012J016.

REFERENCES

[1.] Mayergoyz, I. D., Mathematical Models of Hysteresis, Springer-Verlag, New York, 1991.

[2.] Smith, R. C., Smart Material Systems: Model Development, SIAM, Philadelphia, 2005.

[3.] Dupre, L. and J. Melkebeek, "Electromagnetic hysteresis modelling: From material science to finite element analysis of devices," International Compumag Society Newsletter, Vol. 10, No. 3, 4-15, 2003.

[4.] Preisach, F., "Uber die magnetische nachwirkung," Zeitschrift fur Physik, Vol. 94, 277-302, 1935.

[5.] Jiles, D. C. and D. L. Atherton, "Theory of the magnetisation process in ferromagnetics and its application to the magnetomechanical effect," J. Phys. D: Appl. Phys., Vol. 17, No. 6, 1265-1281, Jun. 1984.

[6.] Jiles, D. C. and D. L. Atherton, "Theory of the magnetisation process in ferromagnetics and its application to the magnetomechanical effect," J. Phys. D: Appl. Phys., Vol. 17, No. 6, 1265-1281, Jun. 1984.

[7.] Mayergoyz, I. D., "Dynamic preisach models of hysteresis," IEEE Trans. Magn., Vol. 24, No. 6, 2925-2927, Nov. 1988.

[8.] Bertotti, G., "Dynamic generalization of the scalar Preisach model of hysteresis," IEEE Trans. Magn., Vol. 28, No. 5, 2599-2601, Sep. 1992.

[9.] Bergqvist, A. J., "A simple vector generalization of the Jiles-Atherton model of hysteresis," IEEE Trans. Magn., Vol. 32, No. 5, 4213-4215, Sep. 1996.

[10.] Leite, J. V., N. Sadowski, P. Kuo-Peng, N. J. Batistela, J. P. A. Bastos, and A. A. de Espindola, "Inverse Jiles-Atherton vector hysteresis model," IEEE Trans. Magn., Vol. 40, No. 4, 1769-1775, Jul. 2004.

[11.] Vecchio, R. D., "An efficient procedure for modelling complex hysteresis processes in ferromagnetic materials," IEEE Trans. Magn., Vol. 16, 809-811, 1980.

[12.] Miano, G., C. Serpico, L. Verolino, and C. Visone, "Comparison of different hysteresis models in FE analysis of magnetic field diffusion," IEEE Trans. Magn., Vol. 31, 1789-1792, 1995.

[13.] Dupre, L., O. Bottauscio, M. Chiampi, M. Repetto, and J. Melkebeek, "Modelling of electromagnetic phenomena in soft magnetic materials under unidirectional time periodic flux excitations," IEEE Trans. Magn., Vol. 35, 4147-4184, 1999.

[14.] Takahashi, N., S. Miyabara, and K. Fujiwara, "Problems in practical finite element analysis using Preisach model," IEEE Trans. Magn., Vol. 35, 1243-1246, 1999.

[15.] Park, G., S. Hahn, S. Lee, and H. Jung, "Implementation of hysteresis characteristics using the Preisach model with MB-variables," IEEE Trans. Magn., Vol. 29, 1542-1545, 1993.

[16.] Dupre, L., R. V. Keer, and J. Melkebeek, "Complementary 2D finite element procedures for the magnetic field analysis using a vector hysteresis model," Intern. J. for Num. Meth. in Eng., Vol. 42, 1005-1023, 1998.

[17.] Matsuo, T., Y. Osaka, and M. Shimasaki, "Eddy current analysis using vector hysteresis models with play and stop hysterons," IEEE Trans. Magn., Vol. 36, 1172-1177, 2000.

[18.] Jin, J.-M., The Finite Element Method in Electromagnetics, 3rd edition, Wiley, Hoboken, NJ, 2014.

[19.] Kuczmann, M., "Using the Newton-Raphson method in the polarization technique to solve nonlinear static magnetic field problems," IEEE Trans. Magn., Vol. 46, No. 3, 875-879, 2010.

[20.] Fujiwara, K., T. Nakata, N. Takahashi, and K. Muramatsu, "Method for determining relaxation factor for modified Newton-Raphson method," IEEE Trans. Magn., Vol. 29, 1962- 1965, Mar. 1993.

[21.] Li, Y. and J.-M. Jin, "A vector dual-primal finite element tearing and interconnecting method for solving 3-D large-scale electromagnetic problems," IEEE Trans. Antennas Propag., Vol. 54, No. 10, 3000-3009, Oct. 2006.

[22.] Li, Y.-J. and J.-M. Jin, "Parallel implementation of the FETI-DPEM algorithm for general 3D EM simulations," J. Comput. Phys., Vol. 228, No. 9, 3255-3267, 2009.

[23.] Li, Y.-J. and J.-M. Jin, "A new dual-primal domain decomposition approach for finite element simulation of 3-D large-scale electromagnetic problems," IEEE Trans. Antennas Propag., Vol. 55, No. 10, 2803- 2810, Oct. 2007.

[24.] Yao, W., J.-M. Jin, and P. T. Krein, "An efficient domain decomposition method for 3-D finite element analysis of nonlinear electric machine problems," 2013 IEEE International Electric Machines & Drives Conference (IEMDC), 709-715, May 2013.

[25.] Yan, S. and J.-M. Jin, "Analysis of nonlinear electromagnetic problems using time-domain finite element method," Proc. IEEE Antennas Propag. Symp., Orlando, FL, Jul. 2013.

[26.] Ising, E., "Beitrag zur theorie des ferromagnetismus," Zeitschrift fur Physik, Vol. 31, No. 1, 253-258, 1925.

[27.] Jiles, D. C., Introduction to Magnetism and Magnetic Materials, Chapman and Hall, New York, 1991.

[28.] Chikazumi, S., Physics of Ferromagnetism, 2nd edition, Clarendon Press, Oxford, 1997, English edition prepared with the assistance of C. D. Graham, Jr.

[29.] Jiles, D. C., J. B. Thoelke, and M. K. Devine, "Numerical determination of hysteresis parameters for the modeling of magnetic properties using the theory of ferromagnetic hysteresis," IEEE Trans. Magn., Vol. 28, No. 1, 27-35, Jan. 1992.

[30.] Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd edition, Cambridge University Press, New York, 2007.

[31.] Ren, Z., "Influence of the R.H.S. on the convergence behaviour of the curl-curl equation," IEEE Trans. Magn., Vol. 32, No. 3, 655-658, May 1996.

[32.] Whitney, H., Geometric Integration Theory, Princeton University Press, Princeton, NJ, 1957.

[33.] Nedelec, J. C., "Mixed finite elements in R3," Numer. Meth., Vol. 35, 315- 341, 1980.

[34.] Webb, J. P., "Hierarchal vector basis functions of arbitrary order for triangular and tetrahedral finite elements," IEEE Trans. Antennas Propag., Vol. 47, No. 8, 1244-1253, Aug. 1999.

[35.] Peterson, A. F., "Absorbing boundary conditions for the vector wave equation," Microw. Opt. Tech. Lett., Vol. 1, No. 2, 62-64, 1988.

[36.] Webb, J. P. and V. N. Kanellopoulos, "Absorbing boundary conditions for the finite element solution of the vector wave equation," Microw. Opt. Tech. Lett., Vol. 2, No. 10, 370-372, 1989.

[37.] Newmark, N. M., "A method of computation for structural dynamics," J. Engineering Mechanics Division. ASCE, Vol. 85, 67-94, Jul. 1959.

[38.] Zienkiewicz, O. C., "A new look at the Newmark, Houboult and other time stepping formulas: A weighted residual approach," Earthquake Engineering and Structural Dynamics, Vol. 5, 413-418, 1977.

[39.] Gedney, S. D. and U. Navsariwala, "An unconditionally stable finite element time-domain solution of the vector wave equation," IEEE Microw. Guided Wave Lett., Vol. 5, No. 10, 332-334, Oct. 1995.

[40.] Saad, Y. and M. H. Schultz, "GMRes: A generalized minimal residual algorithm for solving nonsymmetric linear systems," SIAM J. Sci. Stat. Comput., Vol. 7, No. 3, 856- 869, Jul. 1986.

[41.] G. L. G. Sleijpen and D. R. Fokkema, "BiCGstab(l) for linear equations involving unsymmetric matrices with complex spectrum," Electronic Trans. Numer. Anal., Vol. 1, 11-32, Sep. 1993.

[42.] Carpentieri, B., I. S. Duff, and L. Giraud, "Experiments with sparse preconditioning of dense problems from electromagnetic applications," CERFACS, Tech. Rep. TR/PA/00/04, Toulouse, France, 2000.

[43.] Alleon, G. M. Benzi, and L. Giraud, "Sparse approximate inverse preconditioning for dense linear systems arising in computational electromagnetics," Numer. Algorithms, Vol. 16, 1-15, 1997.

[44.] Hantila, F. I., G. Preda, and M. Vasiliu, "Polarization method for static fields," IEEE Trans. Magn., Vol. 36, No. 4, 672-675, Jul. 2000.

[45.] International Compumag Society, "Testing electromagnetic analysis methods (T.E.A.M.)," http://www.compumag.org/jsite/team.

[46.] Nakata, T., N. Takahashi, and K. Fujiwara, "Summary of results for benchmark problem 10 (steel plates around a coil)," Compel, Vol. 14, No. 2/3, 103-112, Sep. 1995.

[47.] Bottauscio, O., M. Chiampi, C. Ragusa, L. Rege, and M. Repetto, "A test- case for validation of magnetic field analysis with vector hysteresis," IEEE Trans. Magn., Vol. 38, No. 2, 893-896, Mar. 2002.

[48.] Yamada, S., K. Bessho, and J. Lu, "Harmonic balance finite element method applied to nonlinear AC magnetic analysis," IEEE Trans. Magn., Vol. 24, No. 4, 2971-2973, Jul. 1989.

Su Yan (1, 2), * and Jian-Ming Jin (1) (Invited Paper)

Received 10 September 2015, Accepted 12 October 2015, Scheduled 14 October 2015

Invited paper for the Commemorative Collection on the 150-Year Anniversary of Maxwell's Equations.

* Corresponding author: Su Yan (suyan@illinois.edu).

(1) Center for Computational Electromagnetics, Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801-2991, USA. (2) Department of Microwave Engineering, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, China.

Caption: Figure 1. Problem setting and dimensions of TEAM problem 10. (a) Racetrack coil. (b) Steel channels and center plate.

Caption: Figure 2. Problem setup for TEAM problem 10. (a) Constitutive relation curve. (b) Current excitation.

Caption: Figure 3. Numerical results of TEAM problem 10 with a static current excitation. (a) Observation areas on the steel channels. (b) Magnetic flux density along the observation line [L.sub.1] at steady state. (c) Average magnetic flux density versus time in the observation areas [S.sub.1] and [S.sub.3]. (d) Average magnetic flux density versus time in the observation area [S.sub.2].

Caption: Figure 4. Problem setting of the simulation of a demagnetization process. (a) The time profile of the coil current excitation. (b) A typical hysteresis loop of a demagnetization process using Langevin relation.

Caption: Figure 5. Numerical results of the nonlinear demagnetization problem solved by the nonlinear TDFEM with a Newton-Raphson scheme. (a) Observation points on the steel channels. (b)- (d) Demagnetization process and magnetic remanence (small inserts) at the three observation points [P.sub.1]-[P.sub.3].

Caption: Figure 6. Problem setting and dimensions of the ferromagnetic core in TEAM problem 32.

Caption: Figure 7. Magnetic flux density distribution on the ferromagnetic core (only one fourth is shown due to symmetry). The magnitude and three components of the B vector are shown.

Caption: Figure 8. Numerical results of TEAM problem 32 solved by the nonlinear TDFEM with a Newton-Raphson scheme and inverse J-A hysteresis model. (a) Hysteresis loop observed at the center of the center limb. (b) Hysteresis loop observed at the center of the right limb.

Caption: Figure 9. Numerical results of TEAM problem 32 solved by the nonlinear TDFEM with a Newton-Raphson scheme and inverse J-A hysteresis model. (a) Magnetic flux density versus time at the center of the center limb. (b) Magnetic flux density versus time at the center of the right limb.

Caption: Figure 10. Power spectrum distribution of the excitation current and the responding magnetic flux at (a) the center limb and (b) the right limb.