Induction Heating¶
Problem and setup description¶
In this tutorial, a 3D model representing a single-turn coil wrapped around a cylindrical iron core is simulated as a quarter model, leveraging symmetry around the rotational axis. Positioning the iron core at the center of the coil ensures that all parts of the core are equally exposed to the magnetic field. The surrounding air domain is chosen to be sufficiently large to prevent boundary conditions at its edges from influencing the simulation results. The dimensions of the simulated model are illustrated in the following sketch.
The coil is excited by an alternating current, generating a time-varying magnetic field around it. This magnetic field penetrates the conductor inside the coil, in this case, the iron core, inducing eddy currents within the material. These eddy currents circulate through the core, producing heat due to Joule heating. In this tutorial, we focus on analyzing heat conduction within the core region and how to define the coupled magneto-thermal problem.
Defining the problem in openCFS¶
There are multiple ways to define the problem in openCFS, each offering different levels of coupling between the physical phenomena involved. One straightforward approach is forward coupling via sequence steps, which requires two sequence steps. In this case, the magnetic simulation is performed in sequence step 1, where we solve the magnetic field problem and compute the resulting Joule losses. Then, in sequence step 2, the heat conduction problem is solved using the previously calculated Joule losses as a source term on the right-hand side.
Another possibility is to use just one sequence step and couple via couplingList. This has the advantage that you do not have to think about when (at which timesteps) you evaluation the Joule losses, both PDE's use the same timestep in this approach. This approach also works in the fully coupled case, discussed next.
For a fully coupled magneto-thermal simulation, both physics need to be handled within a single sequence step. This is achieved by defining and managing the coupling through a couplingList, which follows the pdeList. We will explore this setup later in the tutorial.
With this approach, we can specify which physical quantities should be updated. The coupling can either involve direct interactions between the magnetics and thermal problem, such as adapting permeability based on temperature, or it can be limited to internal dependencies within a single PDE, such as temperature-dependent thermal conductivity.
Additionally, the simulation can account for convection, allowing for a more realistic thermal analysis.
In this tutorial, we will cover several different cases, including: - Linear transient forward coupling with/without initial temperature - Linear transient forward coupling incorporating convection - Linear harmonic forward coupling - Nonlinear coupling with temperature-dependent B-H curves, where permeability and other parameters can vary with temperature and/or magnetic field.
Hint
For an introduction into the theory behind the magnetic-thermal coupling, visit our explanations for the MagHeatPDE
Linear transient pure forward coupling¶
The first example in this tutorial is the linear transient simulation with a pure forward coupling (first magnetic then heat) using linear_transient_no_init_temp.xml (without an initial temperature, so assuming to start at 0K).
If you are not familiar with magnetic simulations in openCFS yet, visitMagneticPDE.
In our example we have only one coil and we prescribe a current source density in order to excite the system (have a look at our Coils explanation):
<coilList>
<!-- we only have one coil -->
<coil id="coil1">
<!-- prescribe a current source -->
<source type="current" value="i*N*sin(2*pi*f_coil*t)"/>
<part id="1">
<regionList>
<region name="V_coil"/>
</regionList>
<direction>
<!-- use the cylindrical coordinate system that we defined in the domain-tag -->
<analytic coordSysId="cyl1">
<comp dof="phi" value="1.0"/>
</analytic>
</direction>
<wireCrossSection area="A"/>
</part>
</coil>
</coilList>
The structure of the xml file linear_transient_no_init_temp.xml simply consists of two sequence steps
<sequenceStep index="1">
<analysis>
<transient>
<numSteps>100</numSteps>
<deltaT>1e-5</deltaT>
<allowPostProc>yes</allowPostProc>
</transient>
</analysis>
<pdeList>
<magneticEdge>
...
<elemResult type="magJouleLossPowerDensity">
<allRegions outputIds="hdf5"/>
</elemResult>
</storeResults>
</magneticEdge>
</pdeList>
</sequenceStep>
<sequenceStep index="2">
<analysis>
<transient>
<numSteps>100</numSteps>
<deltaT>1e-5</deltaT>
</transient>
</analysis>
<pdeList>
<heatConduction>
...
<bcsAndLoads>
<!-- this is where we get the Joule loss density from the magnetic PDE -->
<heatSourceDensity volumeRegion="V_core" name="V_core">
<sequenceStep index="1">
<quantity name="magJouleLossPowerDensity" pdeName="magneticEdge" />
<timeFreqMapping>
<continuous interpolation="linear"/>
</timeFreqMapping>
</sequenceStep>
</heatSourceDensity>
</bcsAndLoads>
...
</heatConduction>
</pdeList>
</sequenceStep>
Hint
You could also use different timesteps for the second sequence step but keep in mind that the Joule losses are then not consistent anymore!
Hint
Instead of the two-sequence-step method, you can also use the couplingList method from the examples down below - but you will not observe any iterations because the subproblems as well as the global problem is linear.
We observe heat exchange with the environment due to the heat transport boundary condition (for a more detailed explanation on boundary conditions in heat conduction PDE, visit the explanations here for the HeatPDE)
<heatTransport name="S_core_circum" volumeRegion="V_core" bulkTemperature="10" heatTransferCoefficient="15"/>.
Visualization in ParaView¶
Observing the plots down below, we can first of all check the correct specification of boundary conditios of the magnetic problem. We can observe a concentration of the magnetic flux density in the vicinity of the coil - as we expected.
Include an initial temperature¶
The inclusion of an initial temperature is straight forward. We just need to add another sequence step to the xml file (see linear_transient_with_init_temp.xml) including another heat conduction PDE and define the initialValues tag.
...
<sequenceStep index="1">
<analysis>
<static/>
</analysis>
<pdeList>
<heatConduction>
...
<bcsAndLoads>
<temperature name="S_core_circum" value="25"/>
</bcsAndLoads>
...
</heatConduction>
</pdeList>
</sequenceStep>
<sequenceStep index="2">
<analysis>
<transient>
<numSteps>100</numSteps>
<deltaT>1e-5</deltaT>
<allowPostProc>yes</allowPostProc>
</transient>
</analysis>
<pdeList>
<magneticEdge>
...
<storeResults>
<elemResult type="magJouleLossPowerDensity">
<allRegions outputIds="hdf5"/>
</elemResult>
</storeResults>
</magneticEdge>
</pdeList>
</sequenceStep>
<sequenceStep index="3">
<analysis>
<transient>
<numSteps>100</numSteps>
<deltaT>1e-5</deltaT>
</transient>
</analysis>
<pdeList>
<heatConduction>
...
<initialValues>
<initialState>
<sequenceStep index="1"/>
</initialState>
</initialValues>
<bcsAndLoads>
<!-- this is where we get the Joule loss density from the magnetic PDE -->
<heatSourceDensity volumeRegion="V_core" name="V_core">
<sequenceStep index="2">
<quantity name="magJouleLossPowerDensity" pdeName="magneticEdge" />
<timeFreqMapping>
<continuous interpolation="linear"/>
</timeFreqMapping>
</sequenceStep>
</heatSourceDensity>
...
</bcsAndLoads>
...
</heatConduction>
</pdeList>
</sequenceStep>
It can be seen that we now start heating at a temperature of 25°C, which is now set in the simulation file as the initial temperature.
Include convective heat transfer¶
Now we deal with the iron core moving in positive z direction .
In order to simulate the core moving in z direction, a velocity has been added additionally to the heat conduction PDE as follows:
<velocityList>
<velocity name="vel1">
<comp dof="z" value="3" /> <!-- in m/s -->
</velocity>
</velocityList>
After defining the pdeList, we need to set the storage scheme to non-symmetric due to the convective term in the heat conduction PDE. That is done by including:
<linearSystems>
<system>
<solutionStrategy>
<standard>
<matrix storage="sparseNonSym"/>
</standard>
</solutionStrategy>
<solverList>
<pardiso>
<IterRefineSteps>20</IterRefineSteps>
</pardiso>
</solverList>
</system>
</linearSystems>
From the visualisation of this model that was done in ParaView, it can be seen how the core moves along the z axis.
Linear harmonic pure forward coupling with initial temperature¶
If you are not interested in the transient behavior of the heating process, you can also switch to a harmonic magnetic simulation. The Joule losses are evaluated in a consistent way (averaged over one period - thinking in time domain) and represent the steady state of the magnetic system.
Hint
(Linear) harmonic simulations can only handle constant (no solution-dependent) material parameters. Then you would need a multiharmonic simulation, which is out of scope for this tutorial.
This approach has several benefits compared to the transient approach. You could, for example, couple the harmonic magnetic simulation with a steady state heat conduction simulation (represents the steady state), where you can also include convection. Or you can use a transient simulation and observe the transient heat behavior (of course then it is assumed that the magnetic time constants are much smaller compared to the thermal ones). If this sounds like gibberish, have a look at the explanations in MagHeatPDE.
In openCFS, harmonic simulations are typically used to study the steady-state behavior of electromagnetic devices operating under sinusoidal excitation. By assuming a sinusoidal excitation, the simulation can be simplified and the analysis can be focused on the behavior of the system at a specific frequency of interest.
When defining an analysis type in openCFS we are also defining the wanted frequency as follows:
Therefore, it is sufficient to specify only the amplitude of the current or voltage for the coil excitation, since the frequency and phase of the sinusoidal excitation are known.
<coilList>
<!-- we only have one coil -->
<coil id="coil1">
<!-- prescribe a current source -->
<!-- only specify the amplitude - it's a harmonic analysis -->
<source type="current" value="i*N"/>
<part id="1">
<regionList>
<region name="V_coil"/>
</regionList>
<direction>
<!-- use the cylindrical coordinate system that we defined in the domain-tag -->
<analytic coordSysId="cyl1">
<comp dof="phi" value="1.0"/>
</analytic>
</direction>
<wireCrossSection area="A"/>
</part>
</coil>
</coilList>
After defining necessary boundary conditions for magneticEdge PDE, we define heat conduction problem in a new sequence step as a static (or transient) type of analysis and read the initial state from sequence step 1. Setting the boundary conditions and storing results is analogous to the previous examples.
Nonlinear coupling with temperature-dependent parameters¶
Let us start with discussing the solution dependence of the material parameters involved and how they can be defined in the material xml file.
Involved parameters: - Magnetic parameters - magnetic permeability - electric conductivity - Heat parameters - heat conductivity - heat capacity - mass density
In general, each parameter in the material xml file can have a linear and a nonlinear definition. For example, let us have a look at the linear and nonlinear definition of the heat conductivity.
<heatConductivity>
<linear>
<isotropic>
<real>
17
</real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency>temperature</dependency>
<approxType>LinInterpolate</approxType>
<measAccuracy>0.05</measAccuracy>
<maxApproxVal>30</maxApproxVal>
<dataName>Cond_OTS_NC52.fnc</dataName>
</isotropic>
</nonlinear>
</heatConductivity>
Cond_OTS_NC52.fnc is just a plain text file with two columns, the first one lists the temperature and the second one the heat conductivity:
This explanation can be used as a template for electric conductivity, heat conductivity, heat capacity and mass density BUT not for the magnetic permeability. The definition of the magnetic permeability is special because it can depend on the magnetic flux density (classical BH curve) but also on the temperature \mu=\mu(\mathbf{B}, T). If there is no temperature dependence of the BH curve, we can use one of the following two defintions:
Describing an analytic reluctivity as a function of B:
<permeability>
<linear>
<isotropic>
<real> 0.055386 </real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency></dependency>
<approxType>analytic</approxType>
<measAccuracy>0</measAccuracy>
<maxApproxVal>0</maxApproxVal>
<dataName></dataName>
<nuExpr>39.2004*exp(0.00094279*B_R^(11.5695))-21.1455</nuExpr>
<nuDerivExpr>39.2004*0.00094279*11.5695*B_R^(10.5695)*exp(0.00094279*B_R^(11.5695))</nuDerivExpr>
<angle>0</angle>
<zScaling>0.0</zScaling>
</isotropic>
</nonlinear>
</permeability>
Using B-H datapoints from a textfile:
<permeability>
<linear>
<isotropic>
<real>20</real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency>magFluxDensity</dependency>
<approxType>smoothSplines</approxType>
<measAccuracy>1e-2</measAccuracy>
<maxApproxVal>2.0</maxApproxVal>
<dataName>BH100.fnc</dataName>
</isotropic>
</nonlinear>
</permeability>
Now, for the temperature dependence of the permeability, you can extend the previous definitions by cascading multiple permeability definitions with different temperature tags
Describing multiple analytic reluctivities as a function of B and temperature:
<permeability temperature="100">
<linear>
<isotropic>
<real> 0.055386 </real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency></dependency>
<approxType>analytic</approxType>
<measAccuracy>0</measAccuracy>
<maxApproxVal>0</maxApproxVal>
<dataName></dataName>
<nuExpr>39.2004*exp(0.00094279*B_R^(11.5695))-21.1455</nuExpr>
<nuDerivExpr>39.2004*0.00094279*11.5695*B_R^(10.5695)*exp(0.00094279*B_R^(11.5695))</nuDerivExpr>
<angle>0</angle>
<zScaling>0.0</zScaling>
</isotropic>
</nonlinear>
</permeability>
<permeability temperature="200">
<linear>
<isotropic>
<real> 0.055386 </real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency></dependency>
<approxType>analytic</approxType>
<measAccuracy>0</measAccuracy>
<maxApproxVal>0</maxApproxVal>
<dataName></dataName>
<nuExpr>39.2004*exp(0.00094279*B_R^(11.5695))-21.1455</nuExpr>
<nuDerivExpr>39.2004*0.00094279*11.5695*B_R^(10.5695)*exp(0.00094279*B_R^(11.5695))</nuDerivExpr>
<angle>0</angle>
<zScaling>0.0</zScaling>
</isotropic>
</nonlinear>
</permeability>
Describing multiple B-H functions for different temperature:
<permeability temperature="100">
<linear>
<isotropic>
<real>20</real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency>magFluxDensity</dependency>
<approxType>smoothSplines</approxType>
<measAccuracy>1e-2</measAccuracy>
<maxApproxVal>2.0</maxApproxVal>
<dataName>BH100.fnc</dataName>
</isotropic>
</nonlinear>
</permeability>
<permeability temperature="200">
<linear>
<isotropic>
<real>20</real>
</isotropic>
</linear>
<nonlinear>
<isotropic>
<dependency>magFluxDensity</dependency>
<approxType>smoothSplines</approxType>
<measAccuracy>1e-2</measAccuracy>
<maxApproxVal>2.0</maxApproxVal>
<dataName>BH200.fnc</dataName>
</isotropic>
</nonlinear>
</permeability>
Necessary adaptions in the simulation xml file¶
Since se are now dealing with potentially muliple nonlinearities and solution dependencies, we need to take special care about the definitions in the xml file.
Let us first of all carry out a simulation with temperature dependent B-H curves linear_transient_with_initial_tempdepend.xml.
Having a look at the xml file, we first of all, observe that we are now using the couplingList approach, so both magnetic and thermal pde inside of one sequence step.
The most relevant part is to tell openCFS about the considered nonlinearities by defining them in the magnetic PDE and assigning them to the regions where we want to use the nonlinear parameters
...
<regionList>
<region name="V_core" matDependIds="tempperm" nonLinIds="nlperm" />
<region name="V_coil"/>
<region name="V_air"/>
</regionList>
<nonLinList>
<permeability id="nlperm" />
</nonLinList>
<matDependencyList>
<permeability id="tempperm">
<coupling pdeName="heatConduction">
<quantity name="heatTemperature" />
</coupling>
</permeability>
</matDependencyList>
...
nonLinList which handles the BH curve interpolated for a specific temperature (the saturation characteristic of the material). Then we can tell openCFS which quantity shall be used to make the BH curves dependent on. In our case it's the temperature from the heatConduction PDE.
The next step is to handle the numerics that shall be used via couplingList and linearSystems
<!-- define the coupling and the stopping criterion, in this case: iterate until the
flux density of two consecutive iterations changes less than 1e-6 -->
<couplingList>
<iterative>
<!-- logging: write the iteration information to the info.xml file
maxNumIters: maximum number of iterations per timestep. if the stopping criterion was not met, stopOnDivergence defines if
we continue with the next timestep or if we stop the simulation -->
<convergence logging="yes" maxNumIters="20" stopOnDivergence="yes">
<quantity name="magFluxDensity" value="1e-6" normType="rel"/>
</convergence>
</iterative>
</couplingList>
<linearSystems>
<system >
<solutionStrategy>
<standard>
<nonLinear logging="yes" method="fixPoint" >
<lineSearch type="Armijo"/>
<incStopCrit> 1e-1</incStopCrit>
<resStopCrit> 1e-1</resStopCrit>
<maxNumIters> 20 </maxNumIters>
</nonLinear>
</standard>
</solutionStrategy>
<solverList>
<directLU></directLU>
<!-- <pardiso/> -->
</solverList>
</system>
</linearSystems>
The couplingList defines the stopping rule for the iteration between the magnetic and heat PDE, whereas the linearSystems handles the nonlinear magnetic problem itself by using a fixed point rule and Armijo's rule for finding a good linesearch parameter.
Running the simulation linear_transient_with_initial_tempdepend.xml, we can observe several things in the linear_transient_with_initial_tempdepend.info.xml file, which is written automatically for every openCFS simulation.
In the <PDE> section, we can see the number of iterations needed for the coupling of magnetics to heat:
<PDE>
<iterCoupledPDE>
<convergence maxNumIters="20" logging="yes" stopOnDivergence="yes" justNorm="no">
<step number="1" numIters="1" converged="yes">
<quantity name="magFluxDensity">
<iteration count="1" norm="0" converged="yes"/>
</quantity>
</step>
<step number="2" numIters="1" converged="yes">
<quantity name="magFluxDensity">
<iteration count="1" norm="0" converged="yes"/>
</quantity>
</step>
...
In the next step we can also observe the iterations and convergence of the particular magnetic problem in the <nonlinearConvergence> section
<nonlinearConvergence>
<timer wall-clock="0s" wall="0" cpu="0" calls="0" sub="true"/>
<setup method="fixPoint" lineSearch="Armijo" maxIter="20" incStopCrit="0.1" residualStopCrit="0.1" logging="yes"/>
<solStep value="1" analysis="step:1_time:1e-05_coupleIter:0_:magneticEdge">
<couplingStep value="0">
<iteration pdeName="magneticEdge" nr="1" residualErr="0.025245" incrementalErr="1" eta_linesearch="1"/>
<iteration pdeName="magneticEdge" nr="2" residualErr="0.00021273" incrementalErr="0.011307" eta_linesearch="1"/>
</couplingStep>
</solStep>
<solStep value="2" analysis="step:2_time:2e-05_coupleIter:0_:magneticEdge">
<couplingStep value="0">
<iteration pdeName="magneticEdge" nr="1" residualErr="0.50004" incrementalErr="1" eta_linesearch="1"/>
<iteration pdeName="magneticEdge" nr="2" residualErr="0.00047693" incrementalErr="0.39795" eta_linesearch="1"/>
<iteration pdeName="magneticEdge" nr="3" residualErr="5.1602e-07" incrementalErr="0.00014648" eta_linesearch="1"/>
</couplingStep>
</solStep>
...
Let us now have a look at the temperature in paraview and investigate if there is indeed a change in permeability depending on the temperature and the flux density.
It is hard to isolate the influence of the temperature and flux density on the displayed permeability values but it is obvious that the permeability is not constant anymore.
You can find all the used simulation files in the zip folder InductionHeating.zip