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

Randy Clarksean, Ph.D., P.E.

 

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].


∇ × E = - ∂B-- (1) ∂t ∇ × B = μJ (2) ∇ ⋅ B = 0 (3) ∇ ⋅ E = 0 (4)

Where

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)
μ = permeability (vs/a-m)
σ = conductivity (mhos/m)

J is in defined in conductors to be



J = σE
(5)

In a conductor, we can substitute this relationship into Euquation 2 to give us



∇ × B = μJ = μ (σE )
(6)

Then taking the divergence of this relationship.


∇ ⋅ (∇ × B) = ∇ ⋅ (μσE ) (7) 0 = μ σ(∇ ⋅ E ) (8)

which implies



∇ ⋅ E = 0
(9)

Going one step farther with the assumption of constant σ, we have


∇ ⋅ J = ∇ ⋅ (σE ) = σ(∇E ) (10) ∴ ∇ ⋅ J = 0 (11)

Next, we introduce a vector potential A and a scalar potential V as defined here.


 B ≡ ∇ × A (12) ∂A E ≡ - ∇V - --- (13) ∂t

Substituting in the new vector potential A B in Equation 2 and then by identity we have


 ∇ × (∇ × A ) = μJ (14) ∇ (∇ ⋅ A ) - ∇2A = μJ (15)

Recalling that ∇⋅ E = 0 and substituting in the new scalar potential for E gives


 ∇ ⋅ (- ∇V - ∂A-) = 0 (16) ∂t 2 ∂(∇ ⋅ A ) - ∇ V - --------) = 0 (17) ∂t ∇2V = - ∂(∇-⋅-A)) (18) ∂t

Requiring the vector potential A to be conserved, ∇⋅ A = 0, gives us from Equation 15 2A = μJ and in Equation 18 we have 2V = 0.

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



 ∂A- E ≡ - ∂t
(19)

Knowing that 2A = μJ and then solving this equatio for A would alallow the determination of B from



B ≡ ∇ × A
(20)

E could also be determined from



 ∂A E = - --- ∂t
(21)

At present, we do not really need B because we are seeking the power dissipated in the crucible, which can be determined from



P = J ⋅ E
(22)

Now, IF J were independet of time and recalling that


 J = σE (23) ∂A E = - --- (24) ∂t ∂A- J = - σ ∂t (25) ∇2A = - μJ (26)

would imply that AF(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, E0.

Taking the equations above for a conductor and rearranging them gives


 J = σE (27) 1 - ∂A - -∇2A = σ ----- (28) μ ∂t 2 -1-∂A- ∇ A = σ μ ∂t (29)

In an insulator (σ = 0)



∇2A = 0
(30)

In the coil, we assume that J is known (prescribed and time varying)



∇2A = - μJ
(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



A = erAr (r,z,t) + eϕA ϕ(r,z,t) + ezAz (r,z,t)
(32)

Now, going back to the coil equation (Equation 31) and recalling that A is a vector in cylindrical coordinates gives us


-∂- 1--∂-- -1 ∂2Ar- 2-∂A-ϕ- ∂2Ar- ∂r (∂r ∂r(rAr )) + r2 ∂ ϕ2 - r2 ∂ϕ + ∂z = - μJr (33) 2 2 ∂-( 1--∂-(rA )) + -1 ∂-A-θ+ 2-∂Ar- + ∂-A-ϕ = - μJ (34) ∂r ∂r∂r θ r2 ∂ ϕ2 r2 ∂ϕ ∂z2 ϕ ∂ ∂A 1 ∂2A ∂2A ---(r---z) + -2----z2 + ---2z = - μJ ϕ (35) ∂r ∂r r ∂ϕ ∂z

Eliminating all of the derivatives with respect to ϕ leaves


 ∂ 1 ∂ ∂2Ar --(------(rAr)) + ----- = - μJr (36) ∂r ∂r ∂r ∂2z -∂- 1--∂-- ∂-Aϕ- ∂r (∂r ∂r(rA ϕ)) + ∂z2 = - μJϕ (37) ∂ ∂Az 1 ∂2Az ∂2Az ---(r----) + -2 ---2-+ ---2-= - μJϕ (38) ∂r ∂r r ∂ ϕ ∂z

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

 


∂---1--∂- ∂2A-ϕ ∂r(∂r ∂r (rAϕ)) + ∂z2 = - μJ ϕ
(39)

Because we know that Ar and Az are both 0 we can conclude that



 ∂A ∂A ϕ E = - --- = - ----- ∂t ∂t
(40)

Also


 B ≡ ∇ × A (41) 1∂A ∂A Br = ----z - ---ϕ- (42) r ∂ϕ ∂z ∂Ar- ∂Az- Bϕ = ∂z - ∂r (43) 1 ∂ 1 ∂A Bz = -----(rA ϕ) - -----r (44) r ∂r r ∂ ϕ

Once again eliminating all derivatives with respect to ϕ leaves us with


 ∂A ϕ Br = - ----- (45) ∂z B = 1-∂-(rA ) (46) z r∂r ϕ

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
found.


 ∇ ⋅ B = 1--∂-(rB ) + 1-∂B-ϕ+ ∂Bz- (47) r ∂r r r ∂ ϕ ∂z 1 ∂ ∂A ∂ 1 ∂ = ----(- r---ϕ) + ---(------(rA ϕ)) (48) r ∂r ∂z ∂z ∂r ∂r = --1 ∂-(r∂A-ϕ-) + 1-∂-( ∂-(rA )) (49) r ∂r ∂z r∂z ∂r ϕ - 1 ∂ ∂A ϕ 1 ∂ ∂ = --- --(r-----) + ----(---(rA ϕ)) (50) r ∂r ∂z r∂r ∂z = 0 (51)

Recall that we were left with



∂ 1 ∂ ∂2A ϕ --(------(rAϕ)) + ---2- = - μJ ϕ ∂r ∂r ∂r ∂z
(52)

but, this formulation is not ameanable to Galerkin weak formulation, so we let



ψ(r,z,t) ≡ rA ϕ(r,z,t)
(53)

which gives



∂-(-1-(∂ψ-) + -∂-(1∂-ψ) = - μJ ∂r ∂r ∂r ∂z r ∂z ϕ
(54)

This remarkably looks like ∇⋅ (Kψ), but do not forget that we are in cylindrical coordinates. In cylindrical coordinates, ∇⋅ (Kψ), would look like



 1 ∂ ∂ψ ∂ ∂ψ ∇ ⋅ (K ∇ ψ) =------(rK ---) + ---(K ---) ∂r ∂r ∂r ∂z ∂z
(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


 ψ- A ϕ = r (56) 1 ∂ψ Eϕ = - ----- (57) r ∂t B = - 1-∂ψ-= - K ∂ψ- (58) r r ∂z ∂z 1∂-ψ ∂ψ- Bz = r ∂r = K ∂r (59) (60)

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.



J ϕ = Jo(r,z)expiωt
(61)

In the conductor we have J = σE and recall that


 ∂--1-∂ψ- ∂r(r ∂z ) = - μJϕ (62) ∇ ⋅ (K ∇ψ ) = - μJϕ (63)

and that



 1-∂ψ- E ϕ = - r ∂t
(64)

This leads to the following equation in the coil



 iωt ∇ ⋅ (K ∇ψ ) = - μJoexp
(65)

In the conductors we have


∇ ⋅ (K ∇ ψ) = - μ(σE ϕ) = μσ (- 1∂ψ-) (66) r ∂t μσ-∂ψ- ∇ ⋅ (K ∇ ψ ) = r ∂t (67)

And in the dielectrics



∇ ⋅ (K ∇ ψ) = 0
(68)

But ψ is complex and we to revert to real arithmetic we let


 Jϕ = Jocos ωt (69) ψ (r,z,t) = S(r,z)sinωt + C (r,z)cos ωt (70)

This results in two real equations everywhere within the solution domain. For the coil we have


∂ 1 ∂C ∂ 1 ∂C --(---(---) + ---(-----) = - μJo (71) ∂r ∂r ∂r ∂z r ∂z -∂- -1- ∂S- ∂--1-∂S- ∂r (∂r (∂r ) + ∂z(r ∂z ) = 0 (72)

In the conductors this results in


 ∂--( 1-(∂C-) + ∂-(1-∂C-) = μ-σω-S (73) ∂r ∂r ∂r ∂z r ∂z r -∂- -1- ∂S- ∂-- 1∂S- μσ-ω- ∂r (∂r (∂r ) + ∂z( r∂z ) = - r C (74)

Everywhere else we end up with


∂---1- ∂C-- -∂- 1-∂C-- ∂r(∂r (∂r ) + ∂z (r ∂z ) = 0 (75) ∂ 1 ∂S ∂ 1 ∂S ---(---(---) + --(-----) = 0 (76) ∂r ∂r ∂r ∂z r ∂z

Rewriting these equations in tensor notation for cartesian coordinates and a variable K = 1/r gives us (coil, conductors, and elsewhere)


∇ ⋅ (K ∇C ) = - μJo (77) ∇ ⋅ (K ∇S ) = 0 (78)


 μσω ∇ ⋅ (K ∇C ) = ----S (79) μrσω ∇ ⋅ (K ∇S ) = - ----C (80) r


∇ ⋅ (K ∇C ) = 0 (81) ∇ ⋅ (K ∇S ) = 0 (82)

If we solve these equations with the appropriate boundary conditions, they will yield C(r,zand S(r,z). Using S and C, we must now determine the power deposited in the conductive components. Recall that power is



P = J ⋅ E
(83)

In addition, recall that


 J = σE (84) 1∂ ψ E ϕ = - ---- (85) r ∂t ψ (r,z,t) = S (r,z)sin (ωt) + C (r,z)cos(ωt ) (86)

Making substitutions


 J = σE (87) P = J ⋅ E = J E = (σE )E = σE2 (88) ϕ ϕ ϕ ϕ ϕ = σ (- 1-∂ψ-)2 (89) r ∂t σ ∂ψ 2 = -2(---) (90) r ∂t

Determining the time derivative of ψ gives


∂ψ- -∂- ∂t = ∂t ( S (r,z)sin(ωt) + C(r,z) cos(ωt )) (91) = Sω cos(ωt) - C ω sin(ωt ) (92)

Substituting back into the power relationship.


P = σ-( S ω cos(ωt) - Cω sin(ωt))2 (93) r2 σω2 = --2-( S cos(ωt) - C sin(ωt))2 (94) r

To determine the power deposition rate, we need to integrate this over one period.


 ω ∫ 2ωπ Q (r,z) ≡ --- P (r,z,t)dt (95) 2π 0 ω σ ω2 ∫ 2πω 2 2 2 2 = ----2- ( S cos (ωt ) - SC cos(ωt)sin(ωt) + C sin (ωt ))dt (96) 2π r 0

 ω σω2 S2 t 2πω 1 2 2ωπ = 2π--r2-((2ω-sin(ωt)cos(ωt) + 2-)|0 - SC (2-ω sin (ωt))|0 1 t 2π- + C2(- ---cos(ωt )sin (ωt ) + -)|0ω ) 2ω 2 (97)

Which reduces to



 2 Q (r,z) = σω--(S2 + C2) 2r2
(98)

The penetration depth, or “skin” depth of the heating is defined as



 ∘ ----- -2--- δ = μσω
(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.



ρo(∂cn-+ ⃗u ⋅ ∇cn ) = ρo∇ ⋅ (αn∇cn ) + qcn + Rn ∂t
(100)

For a steady problem that does not consider velocity for the source term Rn we have



∇ ⋅ (αn ∇cn ) = - qcn ρo
(101)

Recall the governing equations for the coil, conductors, and all other locations


∇ ⋅ (K ∇C ) = - μJo (102) ∇ ⋅ (K ∇S ) = 0 (103)


 μσω ∇ ⋅ (K ∇C ) = ----S (104) μrσω ∇ ⋅ (K ∇S ) = - ----C (105) r


∇ ⋅ (K ∇C ) = 0 (106) ∇ ⋅ (K ∇S ) = 0 (107)

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.

References


[1]   
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.


[2]   
Gresho, Philip M., 1987, “Model Development for Induction Heating,” Lawrence Livermore National Laboratory Internal Report UCID-21154, 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:40-48.


[4]   
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.


[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:792-826.

Modeling Induction Heating

Randy Clarksean


Randy is a Ph.D., P.E., CFEI, CFII Mechanical Engineer with over 30 years of experience in failure analysis, fires, and forensic engineering. In addition he has expertise in areas of technical due diligence consulting, heat transfer, thermal systems, management, and general technical consulting services.


Post navigation