Efficient frequency-transient co-simulation of coupled heat-electromagnetic problems

• A submitted manuscript is the author's version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
With the recent advent of inductive power charging systems and wireless power transmission in consumer and mobile phone technology, [], all major automotive manufacturers develop concepts to wirelessly charge electric vehicles, both plug-in and pure electric vehicles (EV). For example the prototype from the Leopold Kostal GmbH & Co. KG of such an inductive charging station is shown in Figure . The necessity to charge EVs with today's battery technology after every prolonged use -at least every night -is seen as one of the major drawbacks in the usability of EVs. A system to automate the charging process would reduce the burden on the driver; it could increase the acceptance of EVs, and, in the case of plug-in hybrid EVs, it could help to further reduce the CO  footprint since the battery of the plug-in hybrid could always be considered fully charged. This is important for the calculation of the fleet CO  emission.
A future inductive charging system does not necessarily exhibit a lower efficiency than a comparable conductive charging system, since there are only a few additional compo-http://www.mathematicsinindustry.com/content/4/1/1 (a) Charge station.
(b) Charging of a electric vehicle.

Figure 1 Prototype of an inductive charging station that charges the vehicle through its number plate (Images from Leopold Kostal GmbH & Co. KG).
nents; in a simplified view, the inductive charging system could be considered as a conductive charging system that has been cut in half in the middle of the transformer. There are, however, certain aspects that require attention and detailed design and optimization. These include positioning tolerances of the stationary ('primary') and car-mounted ('secondary') coils, magnetic stray fields, and thermal aspects. The efficiency of both conductive and inductive charging systems is aimed well above %, measured from the primary AC connection to the drive train battery. But even with this high efficiency, at . kW of first generation systems charging power there is a nonnegligible amount of heat to be dissipated. Later generations with even higher power will further increase the heat load on the components. This heat load is the result of several different processes, namely resistance losses due to DC resistance and proximity effects, ferrite core losses and switching losses in the active semiconductor switching components. The main effects appear at the same frequency range as the magnetic field, which is of the order of - kHz. The resulting temperature however changes on much slower timescales, in the order of minutes, determined by the heat capacity and the (relatively large) mass of the involved components. This electromagnetic-thermal problem is fully coupled, as many of the material parameters show a significant thermal dependence. Typical ferrite core losses, for instance, are minimal at temperatures around  • C and increase below and above this temperature. This drives the equilibrium temperature of the ferrite material always close to this temperature, if the dissipated power is small enough, or makes the system thermally unstable, if the heat power is too high.
Engineering samples of such systems are expensive, heavy, possibly dangerous to operate, and take a lot of time to build and optimize. Virtual prototyping using efficient simulation methods accelerates this process. There are different methodologies and models available, [].
The paper is structured as follows: In Section , we propose a particular frequencytransient model for electromagnetic-thermal problems. The electromagnetic (EM) field is considered at high frequencies, where the time scale of the heat conduction is significantly lower than the time scale of the EM field. This can be exploited in the modeling. We propose a co-simulation scheme, similar to dynamic iteration [], and analyze its con-http://www.mathematicsinindustry.com/content/4/1/1 vergence properties in Section . The analysis exhibits interesting results, especially for high frequencies. Possibilities to generalize this model are also discussed. Section  validates the results by a numerical simulation of a simplified model of the inductive charging system shown in Figure .
In contrast to [], where the different ways of co-simulation are discussed, we focus here on modeling and analysis of the frequency-transient model.

Modeling
In this section, we derive a model, which describes the electromagnetic field coupled to the temperature in the materials. For that, in Section . Maxwell's equations are introduced with temperature dependent material parameters. The conduction of the heat is described in Section . by the heat equation together with an electromagnetic power term as source to describe the Joule losses of the EM field. Finally, in Section . assumptions and approximations lead to the frequency-transient model [].

Maxwell's equations
Maxwell's equations read where E and H are the electric and magnetic field strength, D and B are the electric and magnetic flux densities, ρ and J are the electric charge and current densities. These laws are supported by the constitutive relations where the permittivity ε and permeability μ parameters depend only on space while the conductivity σ may also depend on the temperature T. In this way the EM field solution is a function of the temperature distribution (parameter coupling). However, the source current density J src describes an external excitation [, ]; it is assumed to be independent of the temperature. Now to reduce the unknowns in Maxwell's equations, we introduce the magnetic vector potential A and the electric scalar potential ϕ as By using these potentials, () and () are fulfilled automatically. From () we get With Buchholz gauge transformation (∇ϕ = ) this reduces to (  ) http://www.mathematicsinindustry.com/content/4/1/1 The so-called curl-curl equation will be treated in the following on a finite domain with adequate boundary conditions (BC) and initial values (IV). For the low frequency range, where inductive effects dominate, usually the displacement current ∂D/∂t= ω  εÂ can be disregarded. This is called magnetoquasistatic formulation. For details on these formulations we refer to []. Here we are interested in the high frequency range. Therefore our model is based on the full Maxwell formulation ()-().

Heat equation
The classical heat equation describes conduction of heat in materials: with heat conductivity k, density ρ and specific heat capacity c, all constant in time and again on a finite domain with BC and IV. The term Q represents the Joule losses of the EM field. Hence we use a source coupling to connect the heat equation with the EM curlcurl equation (). For simplicity we consider only eddy-current and Joule losses and thus obtain The power term is further simplified in the next section.

Frequency-transient model
Now we aim at a model which allows an efficient simulation. The model consisting of () and () with () is defined in the time domain. A multirate co-simulation scheme could simulate both equations with different time steps. However, for a fast varying source current density the main part of the computational costs is caused by the simulation of (). A discussion of single-rate and multirate schemes can be found in []. We will reduce these costs further by refining the model. Since the temperature is only slowly varying in comparison to the EM field in () the temperature T, where σ is evaluated at, can be averaged bỹ and then () can be approximated by However, this is still in time domain. To allow a solution in frequency domain, we assume a time harmonic source current density whereĴ src is a complex phasor. It follows for μ and ε independent of A that where the complex phasorÂ c is the solution forĴ src . This means the amplitude is constant Let us look at () again, but now with the approximation σ (T(t)) ≈ σ (T i ). The dot product is the usual real inner product.
whereÂ c =Â c (T i ). Now we are interested in the mean power loss in the time interval [τ i , τ i+ ]: and thus all parts with e ±jωt vanish (for interval length being a multiple of the half period length, i.e. τ i+τ i = c π ω , where c ∈ N). Then we are left with The curl-curl equation is formulated with constant material parameters in frequency domain. Thus, only a linear, complex system has to be solved once for each time window, instead for each time step of the curl-curl equation in time domain.

Algorithm
We will now discuss the algorithm to simulate heat-EM problems with the frequencytransient model. After discretization, the model is solved in a Gauss-Seidel scheme, which can be interpreted as co-simulation. It is comparable with a dynamic iteration for time integration. Section . will briefly discuss the co-simulation scheme. In Section . the convergence analysis for the iteration is proved.

Method
The co-simulation scheme uses () and () in a discretized form. For simplicity we assume Finite Integration Technique (FIT) [-] for spatial discretization and then for time discretization the implicit Euler method. As an alternative to FIT discretization, one could use the (lowest order) Finite Element or Boundary Element Methods [, ] with the http://www.mathematicsinindustry.com/content/4/1/1 drawback of non-diagonal material matrices. This would complicate the following derivations. However, in FIT-notation the curl-curl equation () becomes with (diagonal) positive semi-definite matrix for electric conductivity M σ , (diagonal) positive definite matrices for permittivity and reluctivity M , M ν , discrete curl operators C and C . Note, that M l σ := M σ (T l ), where T l denotes the discretized temperature at the lth iteration step. Let n denote the number of nodes in the computational domain, then the matrices are from R n×n ; the discretized (facet-integrated) source current density j s ∈ R n and the (edge-integrated) magnetic vector potential a ∈ R n are ordered in xyzdirection, e.g., a := [ a x , a y , a z ], []. For the heat equation, () gives From Algorithm  it can be seen that time steps are chosen according to the time constant of the heat equation and only one (complex) linear system is solved per iteration for the curl-curl equation in frequency domain. Thus the time step size is independent of the excitation frequency ω. In contrast, a time domain solution would require many time steps per period π/ω and consequently the solution of a large number of real-valued linear systems. This is the reason, why the model is very efficient for high frequencies since there is multirate behavior.

Convergence analysis
We first verify that the Maxwell operator is bounded: Lemma Let the EM problem () be given with adequate BC, IV for a frequency ω > ω  .

Then the inverse of the discrete Maxwell operator exists and is bounded
for some frequency dependent upper bound C  (ω) that is independent of the temperature.
Proof Beforehand we introduce some abbreviations: For ω > ω  for some ω  the matrix X = X(ω) is real, symmetric positive definite and M σ = M σ (T) is real, diagonal, positive semi-definite. Then where X / := U / U - for an eigendecomposition X = U U - . Now let A := X -/ M σ X -/ , which is real and symmetric positive semi-definite. It follows that the eigenvalue decomposition is A = Q - Q, with Q = . Then because A is semidefinite.
The Lemma guarantees solvability of Maxwell's equations independent of the temperature. However, generalizations are possible but we focus here on the high frequency case because it exhibits a distinct multirate potential. In practice the time-harmonic approach can be used over a wide range of possible excitation frequencies. In particular for lowfrequencies where the displacement currents are often disregarded [, ]. Furthermore for other choices of gauging similar results are found. Also in pure transient simulation [], a wide range of excitation frequencies is covered in practice.
Theorem Let the coupled problem ()-() be given with adequate BC, IV for a frequency ω > ω  . We assume for nonlinear materials, e.g., metals, positivity and differentiability for σ (T) w.r.t. temperature T and σ < . Let the exact (monolithic) solution be denoted by a ∞ and T ∞ , then the iteration is convergent for h i small enough with Proof Consider the inner loop of Algorithm  (steps  and ). It consists of () and (). We prove convergence for this inner loop and use the same abbreviations as introduced in () with the shortcut Z l := Z(T l ). Then we deduce from () Z l a l+ = j s ⇔ a l+ = Z l - j s () due to the Lemma above. We define N :=SM kS , use () with T ∞ and a ∞ subtract it from () Next, by adding and subtracting h i (   ) http://www.mathematicsinindustry.com/content/4/1/1

Figure 3 Mimetic discretizations of Maxwell's equations (e.g. by FIT) use a primal and dual grid.
The discrete temperatures T k are located at primary nodes, the line-integrated electric field and the magnetic vector potentials e k = a k are defined at primary edges and finally the facet-integrated source current density j k at dual surfaces. The material laws relate the quantities on both grids, e.g., the conductivity connects by Ohm's Law the currents and voltages, i.e. j k = σ k (T k ) e k . The figure shows the case for a current in x-direction.

Now, let us consider the term
The conductivity matrix M σ is a diagonal matrix for any iteration step. We assume, that the kth component σ k ( ≤ k ≤ n) depends only on the neighboring temperature T k := (P T) k , see Figure . Thus the kth component of the term can be written as Now the mean value theorem can be applied component-wise and yields for the kth component where ζ k between T l+ k and T ∞ k , and σ k is by assumption non-positive. Hence, in matrix vector notation with (ζ ) k = ζ k where A l := diag( a l ) and the diagonal matrix M σ (ζ ) has only non-positive elements, i.e., it is negative semi-definite. It follows Now, we need estimates for (a) the linear term a l+ -a ∞ and (b) the quadratic term a l+ • a l+ -a ∞ • a ∞ . http://www.mathematicsinindustry.com/content/4/1/1 (a) We start with the linear term: Consider () for the difference of exact solution and (l + )th iterate: where we have exploited the Lemma to obtain c(ω) := C   (ω)ω|σ max | j s . (b) We now consider the quadratic term in () Because a • b ≤ a b , this can be written as where the Lemma gives us Using the estimates () and () in () gives Now consider the asymptotic behavior of this for large ω and small h i . It holds Then, for fixed h i and ω large enough, it follows that with some C(ω) ∼  ω  . http://www.mathematicsinindustry.com/content/4/1/1

Generalization
The frequency-transient model can be generalized in different ways. To enhance versatility, one can use a multi-frequency excitation of the EM subsystem. An application would be steel hardening of gears [, ], where two frequencies are necessary to obtain a homogeneous heating of the surface. Also to approximate periodic signals other than sinusoidal, multi-frequency excitation can be used. This also allows for a Harmonic Balance approach, which enables usage of a nonlinear permeability μ.

.. Time dependent phasor
The model can be improved by weakening the assumptions made. In [] it is suggested to consider the complex phasor not as constant within one time window. This means () becomes (   ) http://www.mathematicsinindustry.com/content/4/1/1 (a) d view and d cut view on geometry. From left to right: ferrite (gray), primary coil (blue), air, secondary coil (blue), ferrite (gray), air, steel slice (red). The left coil represents the charging station and the right coil the coil behind the number plate in the car.
(b) Temperature distribution after a simulation time of  min. The maximum temperature is . K. This is the cross section as shown in (a).

Frequency-transient co-simulation example
In this section the frequency-transient model is applied to an inductive charging system for electric vehicles. The charging system is modeled in a d-axisymmetric way by a primary coil for the station and a secondary coil as receiver in the car. Both coils include ferrite. Additionally there is a steel slice to model the part of the car body behind the secondary coil. The electric conductivity of this steel slice depends on the temperature, i.e., the system is mutually coupled. The geometry and simulation results are shown in Figure . For more details of the set up we refer to []. A constant conductivity would have lead to a single way coupled system and thus the magnetic vector potential would by a constant vector phasor. Models with a constant conductivity (according to the initial temperature) will systematically underestimate the heating. In this case the obtained maximum temperature would be  K below the corresponding solution of the mutually coupled problem, [].
For mutually coupled problems the frequency-transient model has proved to be very efficient. This is expected, since the main part of the computational effort -the time integration of the EM subsystem -is avoided. In this numerical example ω was chosen to be  · π ·  kHz. In fact the proposed co-simulation algorithm reached the end time of t end =  min by using only  time steps with n ≤  iterations (on averagen ≈ . iterations). This coarse time grid sufficiently resolves the dynamics of the heat equation. For comparison, in a naive monolithic simulation with  time steps per period  million time steps would be necessary to resolve the dynamics of the EM subsystems. This underlines the computational gain.

Conclusions
A frequency-transient model tailored for heat-electromagnetic problems was derived. The time step size of the coupled system is determined by the heat subsystem only. The con-http://www.mathematicsinindustry.com/content/4/1/1 vergence analysis is presented in detail. Convergence for high frequencies is guaranteed. The error bound for the iteration decreases quadratically with higher frequencies. This result also applies to approaches by Driesen and Hameyer [] and similar implementations in Comsol []. Thus the approach fits perfectly for applications where inductive heating either appears as losses or is intended.