Finite Element Analysis of the Electromagnetics of Continuum

1  Introduction

Electromagnetics of Continuum is defined as a branch of the physical sciences concerned with the interaction of electromagnetic fields with a deformable body. The ultimate goal of this work is to formulate the finite element equations for an electromagnetic continuum. To begin with, notice that the balance laws of electromagnetics continuum consist of two parts: the thermomechanical (TM) part and the electromagnetic (EM) part. The work can be applied to more complex substances, so it can be utilized to treat electrodynamic-related composite materials. Incorporating the coupling between thermomechanics and electromagnetics into a unified theory can further enhance understanding and prediction of large classes of physical phenomena and provide many technological applications. Phenomenologically important cross effects, such as Peltier, Seebeck, Hall, Ettingshausen, Righi-Leduc, and Nernst effects, can now be studied theoretically and numerically.

The electromagnetic (EM) balance laws, i.e., the famous Maxwell’s equations, are introduced. We also introduce the Lorentz Transformation of E and B fields in relativistic electromagnetic theory [1], as well as in its special case, non-relativistic electromagnetic theory.

Then the Eulerian description of the basic laws of continuum theory, i.e., Conservation of Mass, Balance of Linear Momentum, Balance of Angular Momentum, Conservation of Energy, and Entropy Principle, is introduced. Also included in these basic laws are the electromagnetic (EM) part of the body force, body moment, and energy source [2–6]. It is emphasized that this inclusion establishes the link between the thermomechanical part (TM) and the electromagnetic part (EM) of the electrodynamics of continuum.

Then we formulated the constitutive theory for thermo-visco-elastic-plastic-electromagnetic (TVEP-TM) materials. The formulation of constitutive theory, including plasticity, is unique in the sense that one needs to add a set of internal variables to the list of independent constitutive variables and, of course, one needs to supply a set of governing equations for the newly added internal variables [7]. Following that, we specialized the materially linear constitutive equations. It is emphasized that this work is based on geometrically nonlinear, or say large-strain, approaches. Onsager’s postulate is utilized for the derivation of viscosity. Return-Mapping-Algorithm is invoked for plasticity.

It is a well-known physical fact that the electric field E and the magnetic flux B are not independent from each other, nor are the dielectric displacement D and the magnetic field H. To resolve this problem, especially for the finite element formulation, the scalar potential ϕ and the vector potential A are introduced [1], which are related to the electric field E and magnetic flux B as

B≡∇×A,E≡−∇ϕ−1c∂A∂t

Eventually, the governing equations for ϕ and A, not E and B, are obtained in Section 6.

For the finite element formulation, we link the displacement U, temperature T, scalar potential φ, and vector potential A at a generic point in a finite element to the corresponding nodal values associated with that finite element. Finally, the finite element equations for UKβ, Tβ, ϕβ, and Akβ are obtained. It is proposed to solve a torsion problem of a cylinder in two approaches, each of which has two stages. The results will show the coupling between electromagnetics and thermomechanics.

2  Maxwell’s Equations

The balance laws of electromagnetics of continuum consist of two parts: the thermomechanical (TM) part and the electromagnetic (EM) part. The EM balance laws are the well-known Maxwell’s equations that can be expressed in the Heaviside-Lorentz system [1,5,6] as

∇⋅D=q or Dk,k=q(1)

∇×E+1c∂B∂t=0 or eijkEk,j+1c∂Bi∂t=0(2)

∇⋅B=0 or Bk,k=0(3)

∇×H−1c∂D∂t=1cJ or eijkHk,j−1c∂Di∂t=1cJi(4)

where D is the dielectric displacement vector, B the magnetic induction vector, E the electric field vector, H the magnetic field vector, J the total current vector, and q the free charge density.

Define the polarization vector P and the magnetization vector M as

P=D−E or Pk=Dk−EkM=B−H or Mk=Bk−Hk(5)

It is noticed that the quantities q,J,D,E,P,B,H,M are all referred to a fixed laboratory frame RG. On the other hand, q∗,J∗,D∗,E∗,P∗,B∗,H∗,M∗ are referred to a co-moving reference frame RC.

In relativistic electromagnetic theory, the Lorentz Transformation of E and B fields can be expressed as [1]

E∗=γ(E+β×B)−γ2γ+1β(β⋅E)(6)

B∗=γ(B−β×E)−γ2γ+1β(β⋅B)(7)

where β≡vc,β=β⋅β,γ≡11−β2, and v is the velocity of the co-moving frame RC relative to the fixed laboratory frame RG, and c is the speed of light in vacuum. Jackson also explained that the inverse can be found by interchanging the star and the no-star quantities and by changing β to −β, i.e.,

E=γ(E∗−β×B∗)−γ2γ+1β(β⋅E∗)(8)

B=γ(B∗+β×E∗)−γ2γ+1β(β⋅B∗)(9)

Recall that γ≡(1−β2)−12, the Taylor series expansion of γ about β=0 is obtained as

γ≈1+12β2+34β4+....(10)

Let’s keep the Taylor series at most to β2. Then, for non-relativistic electromagnetic theory, we have

E∗=E+β×B,E=E∗−β×B∗(11)

B∗=B−β×E,B=B∗+β×E∗(12)

D∗=D+β×H,D=D∗−β×H∗(13)

H∗=H−β×D,H=H∗+β×D∗(14)

P∗=P−β×M,P=P∗+β×M∗(15)

M∗=M+β×P,M=M∗−β×P∗(16)

Also, it is noticed that Eringen [3] gave a detailed discussion and derivation, and showed that

q∗=q,J∗=J−qv(17)

3  Balance Laws of Continuum Theory

The Eulerian description of the basic laws of continuum theory can be expressed as [4]:

Conservation of Mass

ρ˙+ρvk,k=0(18)

Balance of Linear Momentum

tji,j+ρ(fi−v˙i)+Fi=0(19)

Balance of Angular Momentum

εkijtij+Ck=0(20)

Conservation of Energy

ρe˙−tklvl,k+qk,k−ρh−W=0(21)

Entropy Principle (Second Law of Thermodynamics)

ρη˙+(qiθ),i−ρhθ≥0(22)

where ρ,tkl≠tlk,e,qk,η,θ,vi,fi,Fi,Ck,h, and W are the mass density, Cauchy stress, internal energy density, heat flux, entropy density, absolute temperature, velocity, body force (TM), body force (EM), body moment (EM), energy source (TM), and energy source (EM), respectively.

The EM part of the body force, body moment, and energy source are given as [2–6]:

Fk=qEk+El,kPl+Bl,kMl+1ceklm{JlBm+(vnPlBm),n+∂∂t(PlBm)}(23)

Ck=εijk{PiEj∗+Mi∗Bj}(24)

W=Ek∗(P˙k+Pkvl,l)−Mk∗B˙k+Jk∗Ek∗(25)

Notice that, if A is a function of the Eulerian coordinate xi and time t, then A˙=∂A∂t+A,ivi.

Let the Helmholtz free energy density in this work be defined as

ψ≡e−θη−ρ−1Ek∗Pk(26)

Notice that this expression of Helmholtz free energy density is different from that in classical continuum mechanics, i.e., the inclusion of −ρ−1Ek∗Pk. However, the expression of Helmholtz free energy density is the same as that in micromorphic theory [8]. In other words, the expression of the Helmholtz free energy density is hinged on whether the work has thermomechanical and electromagnetic coupling. Classical continuum mechanics and microcontinuum theory [4] (including micromorphic theory) are two completely different theories, and yet they have many essentials and fundamental principles in common. Therefore, in order to reduce overlap, interested readers are recommended to read [8].

Then one obtains the Lagrangian description of the law of conservation of energy and the Clausius-Duhem (CD) inequality as [7]

ρoe˙=TKLmE˙KL−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗−ρoh(27)

−ρo(ψ˙+ηθ˙)+TKLmE˙KL−1θQKθ,K−PKE˙K∗−MK∗B˙K+JK∗EK∗≥0(28)

where the mechanical part of the Cauchy stress, tklm is defined as

tklm≡tkl−tklem,TKLm≡JtklmXK,kXL,l(29)

And the electromagnetic part of the Cauchy stress, tklem is defined as

tklem≡−PkEl∗−Mk∗Bl(30)

Now one may verify that

εikltklm=εikl{tkl−tklem}=εikl{tkl+PkEl∗+Mk∗Bl}=εikl{tkl+Ci}⇐ Eq.(30)=0⇐ Eq.(20)(31)

which implies the symmetry of the mechanical part of the Cauchy stress tensor, i.e.,

tklm=tlkm,TKLm=TLKm(32)

Notice that tkl≠tlk,tklem≠tlkem. If there is no coupling with electromagnetics, i.e., P=E∗=0, and M∗=B=0, then one has

tklm=tkl,tklem=0(33)

4  Constitutive Equations

In this section, we are going to formulate the constitutive theory for thermo-visco-elastic-plastic-electromagnetic (TVEP-TM) materials. The formulation of constitutive theory, including plasticity, is unique in the sense that one needs to add a set of internal variables to the list of independent constitutive variables and, of course, one needs to supply a set of governing equations for the newly added internal variables [9]. For the electrodynamics of a continuum, a set of internal variables is introduced as [7]

W={Ep,R}(34)

where Ep is the plastic strain corresponding to the Lagrangian strain E defined as

EKL≡12(xk,Kxk,L−δKL)(35)

And R, named as the hardening parameters, is a generalized vector of internal variables.

Now, let the formulation of the constitutive theory for the TVEP-EM materials begin with

Z=Z(Y)(36)

where

Z={Tm,Q,ψ,η,P,M∗,J∗}={TKLm,QK,ψ,η,PK,MK∗,JK∗}(37)

Y={E,E˙,θ,∇θ,E∗,B,Ep,R}(38)

The definitions of TKLm,TKL,QK,EK∗,BK,PK,MK∗, and JK∗ are given as

TKLm≡JtklmXK,kXL,l,TKL≡JtklXK,kXL,l,QK≡JqkXK,k

PK≡JPkXK,k,MK∗≡JMk∗XK,k,JK∗≡JJk∗XK,k(39)

EK∗≡Ek∗xk,K,BK≡Bkxk,K

Following the same procedures in [8], a scalar-valued yield function is defined as

f=f(Y)(40)

For a set of fixed values of W, a hyper surface, named the yield surface, is determined by

f(Y)=0(41)

Define the loading rate λ as the scalar product between the outward normal to the yield surface and the tangent vector to the trajectory of E,E˙,θ,∇θ,E∗,B, i.e.,

λ=∂f∂EKLE˙KL+∂f∂θθ˙+∂f∂E˙KLE¨KL+∂f∂θ,Kθ˙,K+∂f∂EK∗E˙K∗+∂f∂BKB˙K(42)

Three distinct cases, unloading, neutral loading, and loading, are defined as: (a) f<0, (b) f=λ=0, and (c) f=0,λ>0, respectively. The internal variables of plasticity W will remain unchanged in the cases of unloading and neutral loading. The governing equations for the internal variables of plasticity W are proposed to be

W˙=λ∗Φ(Y)(43)

where

λ∗={0iff<0orλ<0λiff=0andλ≥0(44)

Now we substitute Eqs. (36) and (42) into the Clausius-Duhem inequality, Eq. (28). It leads to

−ρo{∂ψ∂EKLE˙KL+∂ψ∂E˙KLE¨KL+∂ψ∂θθ˙+∂ψ∂θ,Kθ˙,K+∂ψ∂EK∗E˙K∗+∂ψ∂BKB˙K+ηθ˙+∂ψ∂W⋅λ∗Φ}+TKLmE˙KL−1θQKθ,K−PKE˙K∗−MK∗B˙K+JK∗EK∗≥0(45)

Since this inequality is linear in E¨KL,E˙K∗,B˙K,θ˙,θ˙,K, it leads to

ψ=ψ(EKL,θ,EK∗,BK,EKLp,R)(46)

η=−∂ψ∂θ(47)

TKLm≡ρo∂ψ∂EKL,TKLm=TKLme+TKLmd(Y)(48)

PK=−ρo∂ψ∂EK∗(49)

MK∗=−ρo∂ψ∂BK(50)

TKLmdE˙KL−1θQKθ,K+JK∗EK∗−ρo∂ψ∂W⋅λ∗Φ≥0(51)

where the superscript d refers to the dissipative part and superscript e refers to the elastic part, or say, the reversible part. Since the inequality, Eq. (51), must be satisfied for any value of the loading rate λ, it implies

∂ψ∂W⋅Φ≤0(52)

TKLmdE˙KL−1θQKθ,K+JK∗EK∗≥0(53)

Also, it should be emphasized that, in the case of loading, a TVEP-EM state leads to another TVEP-EM state. In other words, the consistency condition of plasticity requires that f=0 and λ>0 lead to another state with f=0 [10]. Therefore, in the case of loading, we have f˙=0, i.e.,

f˙=λ∗+∂f∂W⋅W˙=λ∗+λ∗∂f∂W⋅Φ=λ∗(1+∂f∂W⋅Φ)=0(54)

This gives another constitutive constraint to the plasticity

∂f∂W⋅Φ=−1(55)

5  Linear Constitutive Equations

For TVEP-EM solid, define the potential energy density function Σ as

Σ≡ρoψ=Σ(E,θ,E∗,B,Ep,R)(56)

Following the idea of Green and Naghdi [11], and assuming the hardening parameters R do not appear in the potential energy density function, Eq. (56) can now be written as

Σ=Σ(E−Ep,T,E∗,B)≡Σ(Ee,T,E∗,B)=Σo−ρoηoT−12ρoγT2To−A¯KLEKLeT+aKTEK∗+bKTBK+12AKLMNEKLeEMNe+FKLM1EKLeEM∗+FKLM2EKLeBM−12DKL1EK∗EL∗−12DKL2BKBL−DKL3EK∗BL(57)

where To is the reference temperature and T is the temperature variation, and

T≡θ−To,To>0,|T|<<To(58)

From now on, we assume that the material is isotropic. Therefore, all the odd-order material property tensors are vanishing. Then it is straightforward to obtain

η=−∂ψ∂θ=−∂ψ∂T=ηo+γTTo+A¯KLEKLeρo(59)

TKLme=ρo∂ψ∂EKL=∂Σ∂EKL=−A¯KLT+AKLMNEMNe(60)

PK=−ρo∂ψ∂EK∗=DKL1EL∗+DKL3BL(61)

MK∗=−ρo∂ψ∂BK=DKL2BL+DLK3EL∗(62)

Moreover, the 2nd order and 4th order material property tensors are made of a few material constants, i.e.,

A¯KL=A¯δKL,DKL1=D1δKL,DKL2=D2δKL,DKL3=D3δKLAKLMN=λ1δKLδMN+μ11δKMδLN+μ12δKNδLM(63)

Onsager’s Postulate

Now we follow Onsager’s Postulate and construct a potential of dissipation quadratic as [12,13]

D=12dKLMNE˙KLE˙MN+12gMN1θ,Mθ,N+gMN2θ,MEN∗+12gMN3EM∗EN∗(64)

where E˙,∇θ, and E∗ are the thermodynamic forces. Correspondingly, the thermodynamic fluxes Tmd,Q, and J∗ can be obtained as

TKLmd=∂D∂E˙KL=dKLMNE˙MN→d1E˙MMδKL+2d2E˙KL(65)

−1θQM=∂D∂θ,M=gMN1θ,N+gMN2EN∗→g1θ,M+g2EM∗(66)

JM∗=∂D∂EM∗=gNM2θ,N+gMN3EN∗→g2θ,M+g3EM∗(67)

One may verify that the CD inequality Tmd:E¨−1θQ⋅∇θ+J∗⋅E∗≥0 is now reduced to D≥0, which is precisely what the Onsager’s Postulate means.

Return Mapping Algorithm for Plasticity

Notice that E=Ee+Ep. Now, let the hardening parameters R be decomposed into two parts: the vector part G, which has the same dimension as E, and a scalar part g. We propose a yield function as [7]

f=‖E−Ep−G‖−(EY+κHg)=‖Ee−G‖−(EY+κHg)(68)

where EY,H, and κ are material constants. Following the same procedures in [8], let the evolution equations for the plastic strain tensor and the hardening parameters be

E˙p=Ωξ‖ξ‖,G˙=ΩC1(1−κ)Hξ‖ξ‖,g˙=ΩC2(69)

where

ξ≡E−Ep−G=Ee−G,‖ξ‖≡{ξ⋅ξ}12(70)

The trial value of the yield function at t=tn+1, based on elastic prediction, is calculated as

f~n+1=‖En+1−Enp−Gn‖−(EY+κHgn)≡‖ξ~n+1‖−(EY+κHgn)(71)

If f~n+1≤0, which means the elastic prediction is correct, then the plastic strain and hardening parameters at t=tn+1 maintain their values at t=tn, i.e.,

fn+1=‖En+1−Enp−Gn‖−(EY+κHgn)(72)

It is worthwhile to mention that, in this case, i.e., f~n+1≤0,t=tn+1, one has En+1p=Enp, Gn+1=Gn, and gn+1=gn.

If f~n+1>0, which means the elastic prediction is incorrect, then the values of plastic strain and hardening parameters at t=tn+1 should be updated as follows:

En+1p=Enp+ΔΩξ~n+1/‖ξ~n+1‖(73)

Gn+1=Gn+ΔΩC1(1−κ)Hξ~n+1/‖ξ~n+1‖(74)

gn+1=gn+ΔΩC2(75)

where

ΔΩ=f~n+11+C1H−(C1−C2)κH(76)

Now, one can prove that the yield function is based on the updated En+1p,Gn+1,gn+1 returns to zero [7] i.e.,

fn+1=‖E−En+1p−Gn+1‖−(EY+κHgn+1)=0(77)

From Eqs. (60)–(62), it is noticed that the elastic parts of the stresses, TKLme, and PK, and MK∗ depend on the elastic strains, not the total strains. Using the Return Mapping Algorithm, one may find the plastic strains exactly, and once the total strains are found, usually through finite element analysis (to be discussed later), then immediately elastic strains and elastic parts of the stresses as well are obtained. Also, it is seen that the Return Mapping Algorithm is strain-based, not stress-based. Because in this formulation the stresses have two parts: elastic parts and dissipative parts, which depend on strain rates. A strain-based Return Mapping Algorithm works independently of strain rates.

6  Scalar and Vector Potentials

It is a well-known physical fact that the electric field E and the magnetic flux B are not independent from each other, nor are the dielectric displacement D and the magnetic field H. Mathematically, this is noticed from Maxwell’s equations, especially, Eqs. (2) and (4). To resolve this problem, especially in the finite element formulation presented in the next section, let’s introduce the scalar potential ϕ and the vector potential A, which are related to the electric field E and magnetic flux B as

B≡∇×A or Bi=eijkAk,j(78)

E≡−∇ϕ−1c∂A∂t or Ei≡−ϕ,i−1c∂Ai∂t(79)

Now we follow the same pattern as Eqs. (61) and (62), and write the polarization Pk, magnetization Mk, and current Jk as

Pk=D1Ek∗+D3BkMk=D2Bk+D3Ek∗Jk=g2θ,k+g3Ek∗(80)

Now the dielectric displacement vector D and the magnetic field vector H can be expressed as

Dk=Ek+Pk=Ek+D1Ek∗+D3Bk=(1+D1)Ek+D3Bk+D1(Ek∗−Ek)≡εEk+D3Bk+D1(Ek∗−Ek)≈εEk+D3Bk(81)

Hk=Bk−Mk=Bk−(D2Bk+D3Ek∗)=(1−D2)Bk−D3Ek+D3(Ek∗−Ek)≡μ−1Bk−D3Ek+D3(Ek∗−Ek)≈μ−1Bk−D3Ek(82)

Admittedly, the introduction of Eqs. (80)–(82) involve approximations. Without making these approximations, it is difficult, if not impossible, to derive the governing equations for the scalar potential ϕ and vector potential A.

Substituting Eqs. (78) and (79) into Maxwell’s equations leads to

Dk,k=εEk,k+D3Bk,k=−ε{ϕ,kk+1cA˙k,k}=q(83)

∇×E+1c∂B∂t=∇×{−∇ϕ−1cA˙}+1c{∇×A˙}=0(84)

∇⋅B=∇⋅(∇×A)=0(85)

∇×H−1c∂D∂t=∇×[μ−1B−D3E]−1c{εE˙+D3B˙}=μ−1∇×B−εcE˙=μ−1∇×∇×A+εc[∇ϕ˙+1cA¨]=1cJ(86)

It is seen that two of the four Maxwell’s equations, Eqs. (84) and (85), are identically satisfied; the other two equations, Eqs. (83) and (86) can be rewritten as

ϕ,kk+1cA˙k,k=−q/ε(87)

μ−1[Aj,ij−Ai,jj]+εcϕ˙,i+εc2A¨i=1cJi(88)

Following the same procedures in [8], the Lorentz condition is obtained as

ε1cϕ˙+μ−1∇⋅A=0(89)

Substituting the Lorentz condition into Eqs. (83) and (86) results

−ε{∇2ϕ+1c∇⋅A˙}=−ε{∇2ϕ−εμ1c2ϕ¨}=q(90)

−μ−1∇2Ai+1c2εA¨i=Jic(91)

In other words, Maxwell’s equations are reduced to

εμ1c2ϕ¨−∇2ϕ=a≡qε(92)

εμ1c2A¨k−∇2Ak=Ck≡μcJk(93)

This means the Maxwell’s equations become two wave equations with forcing terms, and the speed of wave propagation is equal to

v=c/εμ(94)

Since both the material constants ε and η are positive and greater than one, therefore Eq. (94) says the speed of EM wave propagation in a material is lower than that in vacuum or in air.

7  Finite Element Formulation

Since the mechanical parts of the Lagrangian stress tensor and the Eulerian stress tensor are related as

TKL=JtklmXK,kXL,l,tkl=J−1TKLmxk,Kxl,L(95)

One may verify that the balance law of linear momentum, Eq. (19), can be rewritten as

(TKLmxl,L),K+ρo(fl−v˙l)+JFl=0(96)

The law of conservation of energy can be written as [7]

ρoe˙=TKLmE˙KL−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗−ρoh(97)

Based on the constitutive equations, Eqs. (59)–(62), (65)–(67) and (97) can be further derived to be [7]

ρoγT+ToToT˙+(T+To)A¯KLE˙KLe=TKLmeE˙KLp+TKLdE˙KL−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗+ρoh(98)

Recall the governing equations of the scalar potential ϕ and the vector potential A as

εμc2ϕ¨−∇2ϕ=qε≡a(99)

εμ1c2A¨k−∇2Ak=Ck≡μc{g2θ,k+g3Ek∗}(100)

We now link the displacement, temperature, scalar potential, and vector potential at a generic point in an element to the corresponding nodal values associated with this element as follows:

UK=NαUKαT=NαTαϕ=NαϕαAk=NαAkα(101)

where Nα[α=1,2,3,....,8] are the 8 shape functions for a standard 8-node 3D element. The finite element equations are obtained as (for detailed derivation, see Appendix A)

MαβU¨Mβ=fMα1+fMα2+fMα3(102)

γMαβT˙β=qα1+qα2+qα3+qα4+qα5+qα6+qα7(103)

mαβϕ¨β+kαβϕβ=ξ^α+a~α(104)

mαβA¨kβ+kαβAkβ=A^kα+C^kα(105)

where

Mαβ≡∫VρoNαNβdV=Mβα(106)

γMαβ≡∫VρoγNαNβdV=γMβα(107)

mαβ≡∫v1c¯2NαNβdv(108)

kαβ≡∫vNα,iNβ,idv(109)

and fmα1,fmα2,fmα3,qα1,qα2,qα3,qα4,qα5,qα6,qα7,ξ^α,a^α,A^kα, and C^kα are defined in Appendix A.

8  Torsion Problem

It is noticed that (a) Eq. (96) is mainly a governing equation for displacement u and it is a wave equation, (b) Eq. (98) is mainly a governing equation for temperature T and it is a diffusion equation, (c) Eq. (99) is a wave equation for scalar potential ϕ, (d) Eq. (100) is a wave equation for the vector potential A. Of course, all those differential equations have forcing terms that make u,T,ϕ, and A coupled together—that is why it is coined with the name thermomechanical–electromagnetic coupling.

However, we also noticed that the wave speed for the scalar potential ϕ and the vector potential A is equal to c/εμ which, for most materials, is at the same order as the speed of light; and yet the wave speed for the displacement is just the speed for an acoustic wave in the material; the speed involved in the diffusion equation is related to the thermal conductivity. Therefore, we consider solving the four differential equations, Eqs. (96) and (98)–(100), simultaneously, is numerically difficult (if not impossible) and physically meaningless. We propose to solve this problem in two stages.

In this work, we are using the MKS system, i.e., meter, kilogram, and second are the units for length, mass, and time, respectively. Also, let the units of temperature and electric charge beoK and q, respectively. In the Heavisine-Lorentz system, the dimensions of the electric field and the magnetic field are [E]=[B]=KM/s2/q. The numerical values of the material constants are specified as

A¯=0.002K/M/s2/oKλ1=6K/M/s2,μ11=4.5K/M/s2,μ12=1.5K/M/s2λ¯1=0.03K/M/s,μ¯11=0.0225K/M/s,μ¯12=0.0075K/M/sEY=10000,κ=1,H=5,C1=0.6,C2=0.4γ=400000M2/s2/oK,To=300oKg1=2000KM/s3/oK,g2=10M/s/oK,g3=10/sρo=1000K/M3,ε=2,μ=2(110)