Electromagnetism

In CalculiX, certain types of electromagnetic calculations are possible. These include:

• electrostatic calculations: these can be performed as a special case of thermal calculations, cf. Section 6.9.13.
• magnetostatic calculations. Due to the absence of time derivatives the interaction between electric and magnetic fields drops out and the magnetic equations can be considered on their own.
• magnetic induction calculations. These are calculations without fixed electric charges and in which the displacement current can be neglected. The major industrial application for this type of calculation is inductive heating.

In this section only the last two applications are treated. The governing Maxwell equations run like (the displacement current term was dropped in Equation (735)):

$$\nabla \cdot \boldsymbol{D} = \rho$$ (732)

$$\nabla \times \boldsymbol{E} = - \frac {\partial \boldsymbol{B}}{\partial t}$$ (733)

$$\nabla \cdot \boldsymbol{B} = 0$$ (734)

$$\nabla \times \boldsymbol{H} = \boldsymbol{j}$$ (735)

where $\boldsymbol{E}$ is the electric field, $\boldsymbol{ D}$ is the electric displacement field, $\boldsymbol{B}$ is the magnetic field, $\boldsymbol{H}$ is the magnetic intensity, $\boldsymbol{ j}$ is the electric current density and $\rho$ is the electric charge density. These fields are connected by the following constitutive equations:

$$\boldsymbol{D} = \epsilon \boldsymbol{E}$$ (736)

$$\boldsymbol{B} = \mu \boldsymbol{H}$$ (737)

and

$$\boldsymbol{j} = \sigma \boldsymbol{E}.$$ (738)

Here, $\epsilon$ is the permittivity, $\mu$ is the magnetic permeability and $\sigma$ is the electrical conductivity. For the present applications $\epsilon$ and $\boldsymbol{ D}$ are not needed and Equation (732) can be discarded. It will be assumed that these relationships are linear and isotropic, the material parameters, however, can be temperature dependent. So no hysteresis is considered, which basically means that only paramagnetic and diamagnetic materials are considered. So far, no ferromagnetic materials are allowed.

Due to electromagnetism, an additional basic unit is needed, the Ampère (A). All other quantities can be written using the SI-units A, m, s, kg and K, however, frequently derived units are used. An overview of these units is given in Table 17 (V=Volt, C=Coulomb, T=Tesla, F=Farad, S=Siemens).

Table 17: Frequently used units electromagnetic applications.

symbol	meaning	unit
I	current	A
$\boldsymbol{E}$	electric field	$\frac{\text{V}}{\text{m}} = \frac{\text{kg}.\text{m}}{\text{A}.\text{s}^3}$
$\boldsymbol{ D}$	electric displacement field	$\frac{\text{C}}{\text{m}^2} = \frac{\text{A}.\text{s}}{\text{m}^2}$
$\boldsymbol{B}$	magnetic field	${\text{T}} = \frac{\text{kg}}{\text{A}.\text{s}^2}$
$\boldsymbol{H}$	magnetic intensity	$\frac{\text{A}}{\text{m}}$
$\boldsymbol{ j}$	current density	$\frac{\text{A}}{\text{m}^2}$
$\epsilon$	permittivity	$\frac{\text{F}}{\text{m}} = \frac{\text{A}^2\text{s}^4}{\text{kg}\text{m}^3}$
$\mu$	magnetic permeability	$\frac{\text{kg}.\text{m}}{\text{A}^2.\text{s}^2}$
$\sigma$	electrical conductivity	$\frac{\text{S}}{\text{m}}=\frac{\text{A}^2.\text{s}^3}{\text{kg}.\text{m}^3}$
P	magnetic scalar potential	A
V	electric scalar potential	V$= \frac{\text{kg}.\text{m}^2}{\text{A}.\text{s}^3}$
$\boldsymbol{A}$	magnetic vector potential	$\frac{\text{kg}.\text{m}}{\text{A}.\text{s}^2}$

In what follows the references [89] and [46] have been used. In inductive heating applications the domain of interest consists of the objects to be heated (= workpiece), the surrounding air and the coils providing the current leading to the induction. It will be assumed that the coils can be considered seperately as a driving force without feedback from the system. This requires the coils to be equiped with a regulating system counteracting any external influence trying to modify the current as intended by the user.

Figure 151: Electromagnetic setup for simply connected bodies

Let us first try to understand what happens physically. In the simplest case the volume to be analyzed consists of a simply connected body surrounded by air, Figure 151. A body is simply connected if any fictitious closed loop within the body can be reduced to a point without leaving the body. For instance, a sphere is simply connected, a ring is not. The coil providing the current is located within the air. Turning on the current leads to a magnetic intensity field through Equation (735) and a magnetic field through Equation (737) everywhere, in the air and in the body. If the current is not changing in time, this constitutes the solution to the problem.

If the current is changing in time, so is the magnetic field, and through Equation (733) one obtains an electric field everywhere. This electric field generates a current by Equation (738) (called Eddy current) in any part which is electrically conductive, i.e. generally in the body, but not in the air. This current generates a magnetic intensity field by virtue of Equation (735), in a direction which is opposite to the original magnetic intensity field. Thus, the Eddy currents oppose the generation of the magnetic field in the body. Practically, this means that the magnetic field in the body is not built up at once. Rather, it is built up gradually, in the same way in which the temperature in a body due to heat transfer can only change gradually. As a matter of fact, both phenomena are described by first order differential equations in time. The Ohm-losses of the Eddy currents are the source of the heat generation used in industrial heat induction applications.

From these considerations one realizes that in the body (domain 2, cf. Figure 151; notice that domain 1 and 2 are interchanged compared to [89]) both the electric and the magnetic field have to be calculated, while in the air it is sufficient to consider the magnetic field only (domain 1). Therefore, in the air it is sufficient to use a scalar magnetic potential P satisfying:

$$\boldsymbol{H}=\boldsymbol{T_0} - \nabla P.$$ (739)

Here, $\boldsymbol{T_0}$ is the magnetic intensity due to the coil current in infinite free space. $\boldsymbol{T_0}$ can be calculated using the Biot-Savart relationship [29]:

$$\boldsymbol{T_0}(\boldsymbol{X}) = \frac{1}{4 \pi } \int_{\Omega ... ...\xi } \Vert^3} d \Omega(\boldsymbol{\xi }), \;\;\; \boldsymbol{X} \in \Omega_1.$$ (740)

This integration is computationally very demanding, therefore parallellization is of utmost importance.

The body fields can be described using a vector magnetic potential $\boldsymbol{A}$ and a scalar electric potential V satisfying:

$$\boldsymbol{B}= \nabla \times \boldsymbol{A},$$ (741)

$$\boldsymbol{E}= -\frac{\partial \boldsymbol{A}}{\partial t} - \nabla V.$$ (742)

In practice, it is convenient to set $V=\frac{\partial v}{\partial t}$, leading to

$$\boldsymbol{E}= -\frac{\partial \boldsymbol{A}}{\partial t} -\frac{\partial \nabla v}{\partial t}.$$ (743)

This guarantees that the resulting matrices will be symmetric.

Figure 152: Electromagnetic setup for multiply connected bodies

If the body is multiply connected, the calculational domain consists of three domains. The body (or bodies) still consist of domain 2 governed by the unknowns $\boldsymbol{A}$ and $V$. The air, however, has to split into two parts: one part which is such that, if added to the bodies, makes them simply connected. This is domain 3 and it is described by the vector magnetic potential $\boldsymbol{A}$. It is assumed that there are no current conducting coils in domain 3. The remaining air is domain 1 described by the scalar magnetic potential P.

In the different domains, different equations have to be solved. In domain 1 the electric field is not important, since there is no conductance. Therefore, it is sufficient to calculate the magnetic field, and only Equations (734) and (735) have to be satisfied. Using the ansatz in Equation (739), Equation (735) is automatically satiesfied, since it is satisfied by $\boldsymbol{T_0}$ and the curl of the gradient vanishes. The only equation left is (734). One arrives at the equation

$$\nabla \cdot \mu (\boldsymbol{T_0} - \nabla P) = 0.$$ (744)

In domain 2, Equations (733), (734) and (735) have to be satisfied, using the approach of Equations (741) and (742). Taking the curl of Equation (742) yields Equation (733). Taking the divergence of Equation (741) yields Equation (734). Substituting Equations (741) and (742) into Equation (735) leads to:

$$\nabla \times \frac{1}{\mu } (\nabla \times \boldsymbol{A}) + \sigma \frac{\partial \boldsymbol{A} }{\partial t} + \sigma \nabla V = 0.$$ (745)

The magnetic vector potential $\boldsymbol{A}$ is not uniquely defined by Equation (741). The divergence of $\boldsymbol{A}$ can still be freely defined. Here, we take the Coulomb gauge, which amounts to setting

$$\nabla \cdot \boldsymbol{A}=0.$$ (746)

Notice that the fulfillment of Equation (735) automatically satisfies the conservation of charge, which runs in domain 2 as

$$\nabla \cdot \boldsymbol{j} =0,$$ (747)

since there is no concentrated charge. Thus, for a simply connected body we arrive at the Equations (744) (domain 1), (745) (domain 2) and (746) (domain 2). In practice, Equations (745) and (746) are frequently combined to yield

$$\nabla \times \frac{1}{\mu } (\nabla \times \boldsymbol{A}) - \na... ...{A} + \sigma \frac{\partial \boldsymbol{A} }{\partial t} + \sigma \nabla V = 0.$$ (748)

This, however, is not any more equivalent to the solution of Equation (735) and consequently the satisfaction of Equation (747) has now to be requested explicitly:

$$\nabla \cdot \sigma ( \frac{\partial \boldsymbol{A} }{\partial t} + \nabla V) = 0.$$ (749)

Consequently, the equations to be solved are now Equations (748) (domain 2), (749) (domain 2), and (744) (domain 1).

In domain 3, only Equations (734) and (735) with $\boldsymbol{j}=0$ have to be satisfied (the coils are supposed to be in domain 1). Using the ansatz from Equation (741), Equation (734) is automatically satisfied and Equation (735) now amounts to

$$\nabla \times \frac{1}{\mu } (\nabla \times \boldsymbol{A}) - \nabla \frac{1}{\mu } \nabla \cdot \boldsymbol{A} = 0.$$ (750)

The boundary conditions on the interface amount to:

• continuity of the normal component of $\boldsymbol{B}$
• continuity of the tangential component of $\boldsymbol{H}$, and
• no current flow orthogonal to the boundary, or ( $\boldsymbol{n_i}$ is the normal on domain i, pointing away from the domain):

$$\boldsymbol{B_1} \cdot \boldsymbol{n_1} + \boldsymbol{B_2} \cdot \boldsymbol{n_2} = 0,$$ (751)

$$\boldsymbol{H_1} \times \boldsymbol{n_1} + \boldsymbol{H_2} \times \boldsymbol{n_2} = 0$$ (752)

and

$$\boldsymbol{j_2} \cdot \boldsymbol{n_2}= 0$$ (753)

all of which have to be satisfied on $\Gamma _{12}$. In terms of the magnetic vector potential $\boldsymbol{A}$, electric scalar potential V and magnetic scalar potential P this amounts to:

$$\mu (\boldsymbol{T_0} - \nabla P) \cdot \boldsymbol{n_1} + (\nabla \times \boldsymbol{A}) \cdot \boldsymbol{n_2} =0,$$ (754)

$$(\boldsymbol{T_0} - \nabla P) \times \boldsymbol{n_1} + \frac{1}{\mu _1} ( \nabla \times \boldsymbol{A}) \times \boldsymbol{n_2} = 0$$ (755)

and

$$(\frac{\partial \boldsymbol{A} }{\partial t} + \nabla V)_2 \cdot \boldsymbol{n_2} = 0$$ (756)

on $\Gamma _{12}$. For uniqueness, the electric potential has to be fixed in one node and the normal component of $\boldsymbol{A}$ has to vanish along $\Gamma _{12}$ [89]:

$$\boldsymbol{A} \cdot \boldsymbol{n_2}=0.$$ (757)

To obtain the weak formulation of the above equations they are multiplied with trial functions and integrated. The trial functions will be denoted by $\delta \boldsymbol{A}, \delta V$ and $\delta P$. Starting with Equation(748) one obtains after multiplication with $\boldsymbol{\delta A}$ and taking the vector identies

$$\nabla \cdot (\boldsymbol{a} \times \boldsymbol{b})=(\nabla \time... ...{a}) \cdot \boldsymbol{b} - \boldsymbol{a} \cdot (\nabla \times \boldsymbol{b})$$ (758)

$$\nabla \cdot ( \alpha \boldsymbol{a} ) = \nabla \alpha \cdot \boldsymbol{a} + \alpha \nabla \cdot \boldsymbol{a}$$ (759)

into account (set $\boldsymbol{b}= (\nabla \times \boldsymbol{A})/\mu$ in the first vector identity):

$$\frac{1}{\mu } (\nabla \times \delta \boldsymbol{A}) \cdot (\nabl... ...rac{1}{\mu } (\nabla \cdot \boldsymbol{A}) \nabla \cdot (\delta \boldsymbol{A})$$

$$- \nabla \cdot [ \frac{1}{\mu } (\nabla \cdot \boldsymbol{A} ) \d... ...igma \left( \frac{\partial \boldsymbol{A} }{\partial t} + \nabla V \right) = 0.$$ (760)

Integrating one obtains, using Gauss' theorem (it is assumed that $\Omega_2$ has no free boundary, i.e. no boundary not connected to $\Omega_1$):

$$\int _{\Omega _2} \frac{1}{\mu } (\nabla \times \delta \boldsymbo... ...\delta \boldsymbol{A} \times \nabla \boldsymbol{A}) \cdot \boldsymbol{n_2} dS +$$

$$\int _{\Omega _2} \frac{1}{\mu } ( \nabla \cdot \boldsymbol{A}) (... ...(\nabla \cdot \boldsymbol{A}) \delta \boldsymbol{A} \cdot \boldsymbol{n_2} dS +$$

$$\int _{\Omega _2} \delta \boldsymbol{A} \cdot \sigma \left( \frac{\partial \boldsymbol{A} }{\partial t } + \nabla V \right) d \Omega =0$$ (761)

The trial functions also have to satisfy the kinematic constraints. Therefore, $\delta \boldsymbol{A} \cdot \boldsymbol{n_2}=0$ and the second surface integral is zero.

Applying the vector identity

$$(\boldsymbol{a} \times \boldsymbol{b}) \cdot \boldsymbol{n} = \boldsymbol{a} \cdot (\boldsymbol{b} \times \boldsymbol{n})$$ (762)

and the boundary condition from Equation (755), the integrand of the first surface integral can be written as:

$$\frac{1}{\mu } ( \delta \boldsymbol{A} \times \nabla \times \boldsymbol{A} ) \cdot \boldsymbol{n_2} =$$

$$\frac{1}{\mu } \left ( \delta \boldsymbol{A} \cdot [ ( \nabla \times \boldsymbol{A}) \times \boldsymbol{n_2}] \right) =$$

$$- \delta \boldsymbol{A} \cdot \left[ ( \boldsymbol{T_0} - \nabla \boldsymbol{P}) \times \boldsymbol{n_1} \right].$$ (763)

Consequently, the integral now amounts to:

$$\int _{\Gamma _{12}} \left[ - \delta \boldsymbol{A} \cdot ( \bold... ...) + \delta \boldsymbol{A} \cdot ( \nabla P \times \boldsymbol{n_2}) \right] dS.$$ (764)

Applying the same vector identity from above one further arrives at:

$$\int _{\Gamma _{12}} \left[ - \delta \boldsymbol{A} \cdot ( \bold... ...) + \boldsymbol{n_2} \cdot ( \delta \boldsymbol{A} \times \nabla P) \right] dS.$$ (765)

Finally, using the vector identity:

$$\boldsymbol{n_2} \cdot (\delta \boldsymbol{A} \times \nabla P) = ... ...ldsymbol{A})]-\boldsymbol{n_2} \cdot [ \nabla \times (P \delta \boldsymbol{A})]$$ (766)

one obtains

$$\int _{\Gamma _{12}} \left\{ - \delta \boldsymbol{A} \cdot ( \bol... ...\boldsymbol{n_2} \cdot [ \nabla \times (P \delta \boldsymbol{A}) ] \right\} dS.$$ (767)

The last integral vanishes if the surface is closed due to Stokes' Theorem.

Now the second equation, Equation (749), is being looked at. After multiplication with $\delta V$ it can be rewritten as:

$$\nabla \cdot \left[ \delta V \sigma \left( \frac{\partial \boldsy... ...igma \left( \frac{\partial \boldsymbol{A} }{\partial t} + \nabla V \right) = 0.$$ (768)

After integration and application of Gauss' theorem one ends up with the last term only, due to the boundary condition from Equation (756).

Analogously, the third equation, Equation (744) leads to:

$$\nabla \cdot ( \delta P \mu [ \boldsymbol{T_0} - \nabla P])-\nabla \delta P \cdot \mu (\boldsymbol{T_0} - \mu \nabla P) =0.$$ (769)

After integration this leads to (on external faces of $\Omega_1$, i.e. faces not connected to $\Omega_2$ or $\Omega_3$ the condition $\boldsymbol{B} \cdot \boldsymbol{n_1}=0$ is applied) :

$$\int _{\Omega _1} \nabla \delta P \cdot \mu ( \boldsymbol{T_0} - ... ...{12}} \delta P \mu ( \boldsymbol{T_0} - \nabla P] \cdot \boldsymbol{n_1} dS =0.$$ (770)

Applying the boundary condition from Equation (754) leads to:

$$\int _{\Omega _1} \nabla \delta P \cdot \mu ( \boldsymbol{T_0} - ... ...ma _{12}} \delta P (\nabla \times \boldsymbol{A}) \cdot \boldsymbol{n_2} dS =0.$$ (771)

So one finally obtains for the governing equations :

$$\int _ {\Omega_2} \frac{1}{\mu } (\nabla \times \delta \boldsymbo... ...} (\nabla \cdot \delta \boldsymbol{A} ) (\nabla \cdot \boldsymbol{A} d \Omega +$$

$$\int _{\Omega _2} (\delta \boldsymbol{A}) \cdot \sigma ( \frac{\p... ...amma _{12}} P (\nabla \times \delta \boldsymbol{A} ) \cdot \boldsymbol{n}_2 d S$$

$$= \int _{\Gamma _{12}} \delta \boldsymbol{A} \cdot (\boldsymbol{T_0 } \times \boldsymbol{n }_2) d S$$ (772)

$$\int _ {\Omega _2} \nabla \delta V \cdot \sigma \left ( \frac{\pa... ...mbol{A} }{\partial t} + \nabla \frac{\partial v}{\partial t}\right )d \Omega =0$$ (773)

$$- \int _{\Omega _1} \mu \nabla \delta P \cdot \nabla P d \Omega +... ...mma _{12}} (\delta P) (\nabla \times \boldsymbol{A}) \cdot \boldsymbol{n}_2 d S$$

$$= \int _{\Omega _1} \mu \nabla \delta P \cdot \boldsymbol{T_0} d \Omega.$$ (774)

Using the standard shape functions one arrives at (cf. Chapter 2 in [24]):

$$\sum _e \sum _{i=1}^N \sum _{j=1}^N \left\{ \left[ \int _{{V_{0e}}_2} \varphi _{i,L} \frac{1}{\mu } \varphi _{j,L} \delta _{KM} d V_e \right] \delta A_{iK} A_{jM} - \right.$$

$$\left[ \int _{{V_{0e}}_2} \varphi _{i,M} \frac{1}{\mu } \varphi _{j,K} d V_e \right] \delta A_{iK} A_{jM} +$$

$$\left[ \int _{{V_{0e}}_2} \varphi _{i,K} \frac{1}{\mu } \varphi _{j,M} d V_e \right] \delta A_{iK} A_{jM} +$$

$$\left[ \int _{{V_{0e}}_2} \varphi _{i} \sigma \varphi _j dV_e \delta _{KM} \right] \delta A_{iK} \frac{D A_{jM}}{t} +$$

$$\left[ \int _{{V_{0e}}_2} \varphi _{i} \sigma \varphi _{j,K} dV_e \right] \delta A_{iK} \frac{Dv_j}{Dt} +$$

$$\left. \left[ \int _{{A_{0e}}_{12}} e_{KLM} \varphi _i \varphi _{j,L} n_{2K} dA_e \right] \delta A_{jM} P_i \right \} =$$

$$-\sum_e \sum_{j=1}^N \left[ \int _{{A_{0e}}_{12}} e_{KLM} \varphi _j T_{0L} n_{2K} dA_e \right] \delta A_{jM}$$ (775)

$$\sum _e \sum _{i=1}^N \sum _{j=1}^N \left\{ \left[ \int _{{V_{0e}}_2} \varphi _{i,K} \sigma \varphi _i dV_e \right] \delta v_j \frac{DA_{iK}}{Dt} \right.$$

$$\left. \left[ \int _{{V_{0e}}_2} \varphi _{i,K} \sigma \varphi _{j,K} dV_e \right] \delta v_i \frac{Dv_j}{Dt} \right \} = 0$$ (776)

$$\sum _e \sum _{i=1}^N \sum _{j=1}^N \left\{ \left[ \int _{{V_{0e}}_1} \varphi _{i,K} \mu \varphi _{j,K} dV_e \right] \delta P_i P_j \right.$$ (777)

$$\left. \left[ \int _{{A_{0e}}_{12}} \varphi _i e_{KLM} \varphi _{j,L} n_{2K} dA_e \right] \delta P_i A_{jM} \right \}$$ (778)

$$= - \sum_e \sum_i \left[ \int _{{V_{0e}}_1} \mu \varphi _{i,K} T_{0K} dV_e \right] \delta P_i.$$ (779)

Notice that the first two equations apply to domain 2, the last one applies to domain 1. In domain 3 only the first equation applies, in which the time dependent terms are dropped.

This leads to the following matrices:

$$[K_{AA}]_{e(iK)(jM)}=\int_{{(V_{0e})}_{\Omega _2}} \frac{1}{\mu }... ...L} \delta_{KM} - \varphi_{i,M} \varphi_{j,K} + \varphi_{i,K} \varphi_{j,M}]dV_e$$ (780)

$$[K_{AP}]_{e(jM)(i)}=\int_{{(A_{0e})}_{\Gamma _{12}}} e_{KLM} \varphi_i \varphi_{j,L} n_{1K} dA_e$$ (781)

$$[K_{PA}]_{e(i)(jM)}=\int_{{(A_{0e})}_{\Gamma _{12}}} e_{KLM} \varphi_i \varphi_{j,L} n_{1K} dA_e$$ (782)

$$[K_{PP}]_{e(i)(j)}=- \int_{{(V_{0e})}_{\Omega _1}} \mu \varphi_{i,K} \varphi_{j,K} dV_e$$ (783)

$$[M_{AA}]_{e(iK)(jM)}=\int_{{(V_{0e})}_{\Omega _2}} \varphi_i \sigma \varphi_j \delta_{KM} dV_e$$ (784)

$$[M_{Av}]_{e(iK)(j)}=\int_{{(V_{0e})}_{\Omega _2}} \varphi_i \sigma \varphi_{j,K} \delta_{KM} dV_e$$ (785)

$$[M_{vA}]_{e(j)(iK)}=\int_{{(V_{0e})}_{\Omega _2}} \varphi_{j,K} \sigma \varphi_i \delta_{KM} dV_e$$ (786)

$$[M_{vv}]_{e(i)(j)}=\int_{{(V_{0e})}_{\Omega _2}} \varphi_{i,K} \sigma \varphi_{i,K} \delta_{KM} dV_e$$ (787)

$$\{F_{A}\}_{e(jM)}=- \int_{{(A_{0e})}_{\Gamma _{12}}} e_{KLM} \varphi_j T_{0L} n_{1K} dA_e$$ (788)

$$\{F_{P}\}_{e(i)}=- \int_{{(V_{0e})}_{\Omega _1}} \mu \varphi_{i,K} T_{0K} dV_e$$ (789)

Figure 153: Nonzero parts of the governing matrices

Repeated indices imply implicit summation. The $[K]$ matrices are analogous to the conductivity matrix in heat transfer analyses, the $[M]$ matrices are the counterpart of the capacity matrix. $\{F\}$ represents the force. The resulting system consists of first order ordinary differential equations in time:

$$[M] \{\dot {U} \} + [K] \{U\} = \{F\},$$ (790)

where the corresponding matrices look like in Figure 153. Solution of this system, which is virtually identical to the heat equation and can also be solved as such ([24], Chapter 7), yields the solution for $\boldsymbol{A}, v$ and $P$ (collected in the vector $\{U\}$ in the above equation) from which the magnetic field $\boldsymbol{B}$ and the electric field $\boldsymbol{E}$ can be determined using Equations (741,742).

The internal electromagnetic forces amount to:

$$\{F_A\}_{(iK)} = [K_{AA}]_{e(iK)(jM)} \cdot \{A\}_{(jM)} + [K_{AP}]_{e(iK)(j)} \cdot \{P\}_{(j)}$$

$$= \frac{1}{\mu } \int _{{(V_{0e})}_{\Omega _2}} [\varphi_{i,L} A_{K,L} - \varphi_{i,M} A_{M,K} + \varphi_{i,K} A_{M,M}]dV_e$$

$$= - \int_{{(A_{0e})}_{\Gamma _{12}}} e_{KLM} \varphi_{i,L} P n_{1M} dA_e$$ (791)

and

$$\{F\}_{(i)} = [K_{PA}]_{e(i)(jM)} \cdot \{A\}_{(jM)} + [K_{PP}]_{e(i)(j)} \cdot \{P\}_{(j)}$$

$$= \int_{{(A_{0e})}_{\Gamma _{12}}} \varphi_i e_{KLM} A_{M,L} n_{1K} dA_e$$

$$= - \int_{{(V_{0e})}_{\Omega _1}} \mu \varphi_{i,K} P_{,K} dV_e.$$ (792)

They have to be in equilibrium with the external forces.

What does the above theory imply for the practical modeling? The conductor containing the driving current is supposed to be modeled using shell elements. The thickness of the shell elements can vary. The current usually flows near the surface (skin effect), so the modeling with shell elements is not really a restriction. The current and potential in the conductor is calculated using the heat transfer analogy. This means that potential boundary conditions have to be defined as temperature, current boundary conditions as heat flow conditions. The driving current containing conductor is completely separate from the mesh used to calculate the magnetic and electric fields. Notice that the current in the driving electromagnetic coils is not supposed to be changed by the electromagnetic field it generates.

The volumetric domains of interest are $\Omega_1 \Omega_2$ and $\Omega_3$. These three domains represent the air, the conducting workpiece and that part of the air which, if filled with workpiece material, makes the workpiece simply connected. These three domains are to be meshed with volumetric elements. The meshes should not be connected, i.e., one can mesh these domains in a completely independent way. This also implies that one can choose the appropriate mesh density for each domain separately.

Based on the driving current the field $\boldsymbol{T_0}$ is determined in domain 1 with the Biot-Savart law. This part of the code is parallellized, since the Biot-Savart integration is calculationally quite expensive. Because of Equation (788) $\boldsymbol{T_0}$ is also determined on the external faces of domain 2 and 3 which are in contact with domain 1.

The following boundary conditions are imposed (through MPC's):

On the external faces of domain 1 which are in contact with domain 2 or domain 3:

• Calculation of $\boldsymbol{A}$ based on the external facial values in domain 2 and 3 (cf. the area integrals in $[K]$)

On the external faces of domain 2 and 3 which are in contact with domain 1:

• Calculation of $P$ based on the external facial values in domain 1.
• Imposition of $\boldsymbol{A} \cdot \boldsymbol{n}=0$.

On the faces between domain 2 and 3:

• Continuity of $\boldsymbol{A}$

These MPC's are generated automatically within CalculiX and have not to be taken care of by the user. Finally, the value of V has to be fixed in at least one node of domain 2. This has to be done by the user with a *BOUNDARY condition on degree of freedom 8.

The material data to be defined include:

• the electrical conductivity in the driving coils
• the magnetic permeability in the air
• the density, the thermal conductivity, the specific heat, the electrical conductivity and the magnetic permeability in the workpiece.

To this end the cards *DENSITY, *CONDUCTIVITY, *SPECIFIC HEAT, *ELECTRICAL CONDUCTIVITY, MAGNETIC PERMEABILITY can be used. In the presence of thermal radiation the *PHYSICAL CONSTANTS card is also needed.

The procedure card is *ELECTROMAGNETICS. For magnetostatic calculations the parameter MAGNETOSTATICS is to be used, for athermal electromagnetic calculations the parameter NO HEAT TRANSFER. Default is an electromagnetic calculation with heat transfer.

Available output variables are POT (the electric potential in the driving current coil) and U (vector magnetic potential $\boldsymbol{A}$ in domain 2 and 3, integral of the scalar electric potential $V$ in domain 2 and scalar magnetic potential in domain 1) on the *NODE FILE card and ECD (electric current density in the driving current coil), EMFE (electric field in the workpiece) and EMFB (magnetic field in the air and the workpiece) on the *EL FILE card. In the dat-file the heating power in any element set belonging to domain 2 can be stored (key EBHE).

For Alternating Current (AC) the steady state answer can be calculated by use of the parameter FREQUENCY and by specifying the frequency of the current (in 1/[T], where [T] is the unit of time) by means of the parameter OMEGA on the &sstarf#star;ELECTROMAGNETICS card. In that case the current takes the form

$$\boldsymbol{j}= \boldsymbol{j_0} e^{i \omega t},$$ (793)

and qa solution for the vector $\{U\}$ is looked for with the same time dependency. Substitution in the governing equation and separating real and imaginary part leads to

$$\begin{bmatrix}K & -\omega M \\ \omega M & K \end{bmatrix} \begin{Bmatrix}U_R \\ U_I \end{Bmatrix} = \begin{Bmatrix}F \\ 0 \end{Bmatrix}.$$ (794)

This is a nonsymmetric system of equations the solution of which constitutes the real and imaginary part of the electromagnetic potentials. The fields $\boldsymbol{E}$ and $\boldsymbol{B}$ now satisfy:

$$\boldsymbol{E} = - i \omega [\boldsymbol{A} - {\nabla} v ] e^{i \omega t},$$ (795)

and

$$\boldsymbol{B}= ({\nabla } \times \boldsymbol{A} ) e^{i \omega t},$$ (796)

where $\boldsymbol{E}$, $\boldsymbol{B}$, $\boldsymbol{A}$ and $v$ are now complex quantities. The heating power density amounts to:

$$\overline{\boldsymbol{j}} \cdot \boldsymbol{E} = \sigma \overline... ...{E} = \sigma ( \Vert \boldsymbol{E}_R \Vert^2 + \Vert\boldsymbol{E}_I \Vert^2),$$ (797)

where an overline denotes the complex conjugate. Alternating Current is known to induce the skin effect in structures [29], i.e. the electric field is concentrated in a thin layer underneath the surface. To order to capture this effect a sufficiently fine mesh is needed in that area.

Experience has shown that the matrices inthe above systems of equations are not always well conditioned, especially if the temperature is included and a coupled thermal-electromagnetic system is solved. To improve the conditiono a scaling of the diagonal terms can be introduced. For a diagonal term $a_{ii}$ the scaling term is defined as $a_{ii}^s = \vert a_{ii}\vert$ if $\vert a_{ii}\vert \ge 10^{-30}$ and $a_{ii}^s =1$ else. The original system of equations

$$a_{ij} x_j = b_i$$ (798)

can be written as (no implicit summation)

$$\sum_{j} \frac{a_{ij}}{\sqrt{\vert a_{ii}^s\vert \vert a_{jj}^s\v... ... } ( x_j \sqrt{\vert a_{jj}^s\vert}) = \frac{b_i}{\sqrt{\vert a_{ii}^s\vert} },$$ (799)

or

$$A_{ij} X_j = B_i,$$ (800)

where

$$A_{ij}=\frac{a_{\underline{i}\underline{j}}}{\sqrt{\vert a_{\underline{i}\underline{i}}^s\vert \vert a_{\underline{j}\underline{j}}^s\vert} },$$ (801)

$$X_j=x_{\underline{j}} \sqrt{\vert a_{\underline{j} \underline{j}}^s\vert}$$ (802)

and

$$B_i=\frac{b_{\underline{i}}}{\sqrt{\vert a_{\underline{i}\underline{i}}^s\vert} }.$$ (803)

The latter system is better behaved.

Examples are induction.inp, induction2.inp and induction3.inp.