I was recently asked about developing models for induction heating. So, I decided that I would post the derivation for one method of approaching the heating process. As we know, induction heating is a surface phenomena and can be difficult to solve numerically. So, without further ado – here is that background information.
Development of a model for induction heating
1 Derivation
The development of this derivation was taken from References [1] to [6]. The development leads to a set of equations that can be used to solve for the amount of induction heating within a conductor.
The following equations are a version of Maxwell’s equations. Details of the intial assumptions to reach this state can be found in References [2] and [5].
Where
J  =  free charge current density vector (amp/ m^{2}) 
E  =  electric field intensity (V/m) 
H  =  magnetic field intensity (amp/m) 
B  =  Magnetic flux density (webers/ m^{2}) 
ρ  =  free charge density (coulombs/ m^{3}) 
D  =  electric flux density vector (coulombs/ m^{2}) 
μ  =  permeability (vs/am) 
σ  =  conductivity (mhos/m) 
J is in defined in conductors to be
(5) 
In a conductor, we can substitute this relationship into Euquation 2 to give us
(6) 
Then taking the divergence of this relationship.
which implies
(9) 
Going one step farther with the assumption of constant σ, we have
Next, we introduce a vector potential A and a scalar potential V as defined here.
Substituting in the new vector potential A B in Equation 2 and then by identity we have
Recalling that ∇⋅ E = 0 and substituting in the new scalar potential for E gives
Requiring the vector potential A to be conserved, ∇⋅ A = 0, gives us from Equation 15 ∇^{2}A = –μJ and in Equation 18 we have ∇^{2}V = 0.
In the absence of free charges, it is appropriate to set V = 0, which allows us to neglect the last equation, ∇^{2}V = 0. This then gives us
(19) 
Knowing that ∇^{2}A = –μJ and then solving this equatio for A would alallow the determination of B from
(20) 
E could also be determined from
(21) 
At present, we do not really need B because we are seeking the power dissipated in the crucible, which can be determined from
(22) 
Now, IF J were independet of time and recalling that
would imply that A≠F(t). In addition, this would imply that E=0, which can not be true because if this was true there would be no induction heating. But since we have alternating current power, E≠0.
Taking the equations above for a conductor and rearranging them gives
In an insulator (σ = 0)
(30) 
In the coil, we assume that J is known (prescribed and time varying)
(31) 
Up to this point, the solution has been discussed in general cartesian coordinates. If we were to assume axisymmetric behavior of the fields, the equations can be simplified even further. The coordinate system would depend on r, ϕ, and z. Assuming that the relevant quantities do not vary with ϕ allows us to define A as
(32) 
Now, going back to the coil equation (Equation 31) and recalling that A is a vector in cylindrical coordinates gives us
Eliminating all of the derivatives with respect to ϕ leaves
By neglecting the helicity (pitch) of the coils as well as 3D effects, along with appropriate boundary conditions gives A_{r} = 0 and A_{z} = 0. This reduces the previous set of equations to

(39) 
Because we know that A_{r} and A_{z} are both 0 we can conclude that
(40) 
Also
Once again eliminating all derivatives with respect to ϕ leaves us with
It can be shown that this implies that ∇⋅ B = 0 by examining the definition of this relationship in conjuction with the relationships for B_{r} ande B_{z} that were just
found.
Recall that we were left with
(52) 
but, this formulation is not ameanable to Galerkin weak formulation, so we let
(53) 
which gives
(54) 
This remarkably looks like ∇⋅ (K∇ψ), but do not forget that we are in cylindrical coordinates. In cylindrical coordinates, ∇⋅ (K∇ψ), would look like
(55) 
where, K = 1/r
We will return to the “form” of the equation later as it relates to being solved within a computational fluid dynamics package.
Now, assuming that a solution for ψ(r,z,t) could be obtained, we could post process the solution to obtain
where, E_{ϕ} is the only item needed to solve for the power density in a crucible.
For the type of engineering problem to be solved here, we want to remove the timedependce of the problem by seeking the quasisteady (periodic or timeharmonic) solution associated with the periodic driving force – the alternating current in the coil.
(61) 
In the conductor we have J = σE and recall that
and that
(64) 
This leads to the following equation in the coil
(65) 
In the conductors we have
And in the dielectrics
(68) 
But ψ is complex and we to revert to real arithmetic we let
This results in two real equations everywhere within the solution domain. For the coil we have
In the conductors this results in
Everywhere else we end up with
Rewriting these equations in tensor notation for cartesian coordinates and a variable K = 1/r gives us (coil, conductors, and elsewhere)
If we solve these equations with the appropriate boundary conditions, they will yield C(r,z) and S(r,z). Using S and C, we must now determine the power deposited in the conductive components. Recall that power is
(83) 
In addition, recall that
Making substitutions
Determining the time derivative of ψ gives
Substituting back into the power relationship.
To determine the power deposition rate, we need to integrate this over one period.
(97) 
Which reduces to
(98) 
The penetration depth, or “skin” depth of the heating is defined as
(99) 
2 FIDAP Implementation
As shown above, there are two coupled equations that must be solved in order to determine the power deposition. The numerical solution of these equations is needed in order to apply a source term within the energy equations.
These equations have previously solved in FIDAP. That implementation used modified versions of the momentum and energy equations to provide a mechanism for the solution of two coupled equations. Currently, we want to solve for the induction heating field in addition to the flow field and the energy equation. In order to do this, a mechanism has to be defined within FIDAP to solve these equations.
It is possible to solve for the transport of 15 different species within FIDAP. The goal is to be
able to use two of the species transport equations to determine C and S, which would allow the calculation of a source term within the energy equation.
The species equations within FIDAP, take the following form in vector notation.
(100) 
For a steady problem that does not consider velocity for the source term R_{n} we have
(101) 
Recall the governing equations for the coil, conductors, and all other locations
This results in a situation where the source term q_{cn} within FIDAP has to be set to equal the righthandside of the governing equations. The diffusivity α_{n} is set equal to K which is defined as 1/r – a variable.
References
[1] Gresho, Philip M., and Derby, Jeffrey, J., 1987, “Solution of the TimeHarmonic SemiMaxwell Equations for Induction Heating Using FIDAP,” First FIDAP User’s Group Meeting, Summer, Evanston, Illinois.
[2] Gresho, Philip M., 1987, “Model Development for Induction Heating,” Lawrence Livermore National Laboratory Internal Report UCID21154, August.
[3] Gresho, P. M., and Derby, J. J., “A finite element model for induction heating of a metal crucible,” Journal of Crystal Growth, 85:4048.
[4] Atherton, L.J., and Martin, R.W., 1988 “Modeling Induction Heating and 3D Heat Transfer for Growth of Rectangular Crystals Using FIDAP,” Second FIDAP User’s Group Meeting, Summer, Evanston, Illinois.
[5] Jackson, J.D., 1975, Classical Electrodynamics, John Wiley and Sons, New York, Second Edition.
[6] Derby, J. J., Atherton, L. J., and Gresho, P. M., 1989, “An integrated process model for the growth of oxide crystals by the Czochralski method,” Journal of Crystal Growth, 97:792826.