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
|J||=||free charge current density vector (amp/ m2)|
|E||=||electric field intensity (V/m)|
|H||=||magnetic field intensity (amp/m)|
|B||=||Magnetic flux density (webers/ m2)|
|ρ||=||free charge density (coulombs/ m3)|
|D||=||electric flux density vector (coulombs/ m2)|
J is in defined in conductors to be
In a conductor, we can substitute this relationship into Euquation 2 to give us
Then taking the divergence of this relationship.
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
In the absence of free charges, it is appropriate to set V = 0, which allows us to neglect the last equation, ∇2V = 0. This then gives us
Knowing that ∇2A = –μJ and then solving this equatio for A would alallow the determination of B from
E could also be determined from
At present, we do not really need B because we are seeking the power dissipated in the crucible, which can be determined from
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)
In the coil, we assume that J is known (prescribed and time varying)
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
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 Ar = 0 and Az = 0. This reduces the previous set of equations to
Because we know that Ar and Az are both 0 we can conclude that
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 Br ande Bz that were just
Recall that we were left with
but, this formulation is not ameanable to Galerkin weak formulation, so we let
This remarkably looks like ∇⋅ (K∇ψ), but do not forget that we are in cylindrical coordinates. In cylindrical coordinates, ∇⋅ (K∇ψ), would look like
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 time-dependce of the problem by seeking the quasi-steady (periodic or time-harmonic) solution associated with the periodic driving force – the alternating current in the coil.
In the conductor we have J = σE and recall that
This leads to the following equation in the coil
In the conductors we have
And in the dielectrics
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
In addition, recall that
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.
Which reduces to
The penetration depth, or “skin” depth of the heating is defined as
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.
For a steady problem that does not consider velocity for the source term Rn we have
Recall the governing equations for the coil, conductors, and all other locations
This results in a situation where the source term qcn within FIDAP has to be set to equal the right-hand-side of the governing equations. The diffusivity αn is set equal to K which is defined as 1/r – a variable.
 Gresho, Philip M., and Derby, Jeffrey, J., 1987, “Solution of the Time-Harmonic Semi-Maxwell Equations for Induction Heating Using FIDAP,” First FIDAP User’s Group Meeting, Summer, Evanston, Illinois.
 Atherton, L.J., and Martin, R.W., 1988 “Modeling Induction Heating and 3-D Heat Transfer for Growth of Rectangular Crystals Using FIDAP,” Second FIDAP User’s Group Meeting, Summer, Evanston, Illinois.