Simulation of Coupled Electromagnetic/Thermal Systems using CAE Software
Executive Summary
Electromagnetic induction heating occurs in electrically conducting materials when they are in the presence of time varying magnetic fields. This phenomenon is utilized in a wide range of industrial applications.
Simulating induction heating effects using Computer Aided Engineering (CAE) software requires the solution of both electromagnetic and thermal field distributions. Often a transient solution is required to determine the temperature distribution as a function of time, and this can be complicated by the characteristic wide disparity between the time constants of the electrical and thermal systems.
This paper will illustrate examples of coupled electromagnetic/thermal field simulations for both steady-state and transient conditions.
Overview
We will consider the simplest cases of coupled systems where the material properties will be assumed to be essentially constant over the range of operating temperatures. By using simple geometric models, we will be able to compare numerical simulation results with analytic solutions.
In order to develop an understanding of coupled electromagnetic/thermal simulations, it is easier to first consider the component systems separately. We will begin by studying the characteristics of the individual thermal and electromagnetic solutions before undertaking combined simulations.
Our thermal model will consist of a copper sphere with a diameter of 2 cm. The
picture shows a cutaway view with one quarter of the sphere removed. The exposed surfaces show contour plots of the temperature variation within the sphere for a uniformly distributed volume heat source.
We will assume initially that the sphere is uniformly heated and examine both steady-state and transient solutions for convection heat transfer from the sphere surface.
Our electromagnetic model will utilize the same sphere with the addition of a surrounding coil carrying 60 Hz current.
The picture shows a cutaway view of both the sphere and coil. The cones inside the coil indicate current flow produced by an external power source. The J Field contour plot on the exposed sphere surfaces shows the density of the resulting induced currents within the sphere.
We will use the electromagnetic solution to obtain the steady-state electrical power transferred to the sphere and also examine the transient start-up conditions that occur when the coil is first energized.
Finally, we’ll turn to coupled models and again perform both steady-state and transient simulations. In addition to transient results for convection heat transfer, results will be presented for radiation heat transfer and for systems where both modes are present.
Rotationally Symmetric Models
Since both the sphere and surrounding coil are solids of revolution about the same axis, a full 3D analysis is unnecessary because the models are Rotationally Symmetric (abbreviated as RS).
At right we show the equivalent RS model for the steady-state induced current solution. Note that the semicircle represents a full 3D sphere, and that the rectangle represents the coil. The J Field 
contour patterns for the RS model are the same as those shown previously for the 3D model.
Because of the Rotational Symmetry, the field solutions are equivalent to full 3D models. We will use this fact to our advantage and display simulation results for RS models throughout the rest of this paper.
Thermal Models
Model Parameters
The physical properties for copper will be assumed to have constant values of 401 W/m*K for the thermal conductivity, 384 J/kg*K for the specific heat, and 8960 kg/m^3 for the mass density. We will assume that the sphere is initially at a constant ambient temperature of 20 C, and the only heat transfer mode is surface convection with a local convection coefficient h = 10 W/m^2 *K.
We will model the temperature rise produced by a volume heat source of 2.5 W and examine both steady-state and transient solutions.
Static Thermal Simulations
The thermal simulations shown in this paper were obtained using the KELVIN (2D/RS) CAE software package developed by INTEGRATED Engineering Software.
If a uniform volume heat source of 2.5 W is applied, it can be analytically calculated that the surface temperature should rise to an equilibrium value of about 218.944 C. The temperature distribution inside the sphere can also be calculated analytically, and it will reach a maximum value of about 218.968 C at the center.
Shown at right are contour plots of the temperature distribution within the sphere obtained from a KELVIN static field solution. The results were obtained from Finite Element Method (FEM) analysis using Quadratic basis elements, and are in excellent agreement with the analytic predictions.
The difference between the center and surface temperatures is only 0.0248 C due to the small size of the sphere and the high thermal conductivity of the copper.
Transient Thermal Simulations
It was noted from the static solution shown in the previous section that the temperature within the sphere is extremely uniform. This fact will lead to a simplification of the analytic transient solution as we will explain in this section.
A completely general thermal transient would allow for temperature gradients within the sphere. An analytic solution is possible, but it is complex and can only be expressed as the sum of an infinite series.
Fortunately, when temperature gradients are small, the transients can be approximated using the lumped capacitance method [1] which reduces the analysis to essentially the same mathematical form as that of a common RC electrical circuit transient. For thermal systems, the corresponding resistance is equal to the product of the local convection coefficient and the boundary surface area. Similarly, the corresponding thermal capacitance is equal to the product of the mass and the specific heat.
A quantitative test used to assess the validity of lumped capacitance analysis for any particular case is obtained by calculating the Biot number (abbreviated as Bi) based on the convection, conduction and dimensional properties of the thermal system. For our sphere, the Biot number is calculated as:
Bi = (local convection coefficient*characteristic length)/thermal conductivity
= 2.494E-04
Here we have used the radius of the sphere for the characteristic length.
It is generally accepted that the lumped capacitance approximation will be sufficiently accurate for Biot numbers less than 0.1, so our system will more than satisfy this criteria.
Using the lumped capacitance approach, we can calculate the thermal time constant (τ) of our system using:
τ= (mass*specific heat)/(local convection coefficient*surface area)
= 1146.88 (s)
= 0.3186 (hrs.)
Note that we have shown the time constant in units of both seconds and hours. Because of the extremely long thermal time constant, we will find it more convenient to display results using hours for the time unit.
The graph below shows a comparison of the transient results obtained from KELVIN (using a time step of 0.012 hours) and analytic results calculated using the lumped capacitance method. For all practical purposes the results are identical.
Electromagnetic Models
Model Parameters
The electrical conductivity for copper will be assumed to be a constant value of 5.84E+07 S/m. The surrounding coil has an inner diameter of 5 cm, an outer diameter of 11 cm, and a height of 5 cm. The coil has 1000 turns with effective inductance of 50.90 mH and winding resistance of 5 ohms. As a result the L/R time constant for the coil is 10.18 ms.
Steady-State 60 Hz Simulations
The electromagnetic simulations shown in this paper were obtained using the OERSTED (2D/RS) CAE software package developed by INTEGRATED Engineering Software.
The 60 Hz steady-state solution was obtained by setting OERSTED to Harmonic/Single Frequency mode. The field solution for this model was obtained using Boundary Element Method (BEM) analysis.
To produce the same 2.5 W power in the sphere, it was found that a 60 Hz coil voltage of 
128.272V RMS was required (and this is commensurate with a coil current of 6.468 A RMS).
However, the power density is not at all uniform within the sphere. The contour plot at right shows that the power density increases with distance from the RS axis.
Transient Start-Up Simulations
Previous calculations showed the time constant of the thermal system was 1146.88 s which is more than five orders of magnitude greater than the electrical time constant of 10.18 ms. As such the time step needed to accurately calculate the electrical transient must be greatly reduced. The transient analysis results were obtained using the Finite Element Method (FEM) analysis with Quadratic basis elements.
The graph below shows the transient results obtained from OERSTED for a 60 Hz sinusoidal applied voltage with zero initial value using a time step of 0.2 ms. Both coil current and power transferred to the sphere reach their steady-state conditions in about 50 ms.
Note that the steady-state power loss is a time varying signal. In fact it is the sum of a constant average value plus a 120 Hz sinusoidal signal. It can be shown mathematically that the power signal for linear systems will always oscillate at twice the frequency of the current (or voltage) signals.
Coupled Electromagnetic/Thermal Models
Having gained an understanding of the different characteristics of the separate thermal and electromagnetic solutions, we are now ready to investigate the behavior of coupled systems. We will perform coupled simulations using the INDUCTO (2D/RS) CAE software package developed by INTEGRATED Engineering Software. INDUCTO contains both the OERSTED and KELVIN field solvers, and also provides coupling between the solutions electromagnetic and thermal field solutions.
Static Thermal/Time Harmonic Eddy-Current Simulations
For our first coupled simulation, we will use the induced eddy current power loss from a steady-state 60 Hz electromagnetic solution as a heat source for a static thermal analysis. INDUCTO automatically calculates the time average power density from the eddy-current solution and transfer this into a heat distribution for the thermal analysis.
Recall that in an earlier section we had seen that the eddy currents induced in the sphere produce a non-uniform power density distribution that increases with distance from the RS axis. We would expect that the resulting steady-state temperature distribution would be different from the earlier KELVIN model that was solved for a uniformly distributed volume heat source.
The first picture below left shows the conductor loss power density contours obtained from the eddy current solution that we had seen previously from our OERSTED model. The resulting steady-state temperature contours calculated from an INDUCTO simulation are shown in the center picture. For comparison, the last picture below right shows the contours produced by a uniformly distributed volume heat as we had seen in our KELVIN model.
The temperature contours produced by the uniform volume heat source are simple concentric spheres, but the non-uniform eddy current heating produces a more complex distribution.
Note though that because of the high thermal conductivity of copper the temperature throughout the sphere in both cases is essentially uniform.
Transient Thermal/Time Harmonic Eddy-Current Simulations
The previous section showed that the steady-state non-uniform eddy current heating results are indeed different from the uniform volume heat results, but for practical purposes both produce an extremely constant temperature distribution.
Not surprisingly, the transient results are also similar as shown in the plot below.
Here the transients are computed assuming only convection heat transfer.
Comparison of Transient Results for Different Heat Transfer Modes
For all the cases considered so far, the only mode by which the sphere could transfer heat to its surrounding has been convection. However at higher temperatures radiation can be the dominant mode, and in vacuum devices convection cannot occur, so radiation would be the only mode of heat transfer.
In order to investigate the effects of heat transfer modes on transient behavior, models were constructed that would have the same steady-state temperature for equivalent heat sources regardless of heat transfer mode. To do this, a radiation only model was constructed with a surface emissivity (ε) of 0.9. The amount of heat input required to raise the sphere temperature from 20 C to 220 C was calculated to be 3.3194 W. The same power was used as the input for a second model which had only convection (h=13.21) heat transfer, and a third which had equal convection and radiation (h= 6.604; ε=0.45).
The graph below compares INDUCTO and lumped capacitance results for all three cases. Here we use the abbreviation “LC” to indicate a lumped capacitance calculation.
It is interesting that heat transfer by radiation produces the fastest transient.
Summary
Coupled thermal/electromagnetic software provides significant benefits for designers of induction heating systems. Most significant among those benefits are (1) it circumvents the possibility of geometry errors in having to port geometry between programs, and (2) the time saved in having to work through a problem setup twice. The exercises in this white paper should give the reader a heightened confidence in the software's accuracy. These types of problems require very tight meshes especially when the skin depth is small. INTEGRATED prides itself in its auto-meshing capability in problems of this nature. In addition, the dual code formulation always gives the user the ability to check the fidelity of global quantities like thermal dissipation without software changes simply by changing from boundary element to finite element analysis for the eddy current solution. High frequency induction heaters represent a class of devices that require a great disparity in cell size to accurately predict the electromagnetic solution. The thermal dissipation is the all important quantity required by the thermal solver. The calculations here show the integrity of thermal solver's performance both in radiative and convective conditions. Finally, note that the analytical solution was available in this problem because no thermal gradients existed in the source. But the numerical solver is set up to readily solve for any source or media profile.
References
[1] Incropera, F. P., and D. P. DeWitt, Fundamentals of Heat and Mass Transfer, 3rd ed., John Wiley & Sons, 1990.
About INTEGRATED Engineering Software
Since its inception in 1984, INTEGRATED Engineering Software has created simulation tools that reflect the inspiration of our customers: thousands of engineers and scientists who, everyday, push the boundaries to envision what is possible. They take their ideas from a realm that is almost science fiction and bring them to reality.
As the name of our company suggests, all our programs are seamlessly integrated, starting from a concept, through entry of the geometry and physics of the problem, to the selection of type of solver and the problem's solution. Once the problem has been solved, a vast number of parameters can be calculated or the field quantities displayed.
INTEGRATED Engineering Software is a leading developer of hybrid simulation tools for electromagnetic and particle trajectory analysis. We provide a complete line of fully integrated 2 and 3 dimensional simulation software.
Since the creation of our company, our focus has always been here and our experience has grown hand-by-hand with a great recognition in our market.
INTEGRATED is staffed with leading R&D engineers in areas such as electrical engineering, magnetics, and high frequency applications. Our tools are used in a wide variety of industries, including manufacturing, automotive, medical, telecommunications, power, health care and aerospace markets, as well as universities and research laboratories.
INTEGRATED products allow engineers and scientists to reduce design cycles, save time and money and deliver more efficient products to the market faster than ever before.
INTEGRATED empowers engineers and scientists with many options to choose from:
The best solvers for each specific application: Boundary Element Method (BEM), Finite Element Method (FEM) or Finite Difference Time Domain (FDTD) solvers. The best optimization tool for each particular design: parametric analysis, scripting or application programming interface (API).
INTEGRATED’s commitment is to provide designers with the most sophisticated analysis tools to assist them in the creation of the future.
Contact us for an evaluation
Send us your model, whatever the level of complexity. We will show you how to get results from your exact design – no packaged demos.
Contact us for an evaluation and start improving productivity today. A live demo is also available.
Phone: +1.204.632.5636
Fax: +1.204.633.7780
Email: info@integratedsoft.com
Website: www.integratedsoft.com
