1 Introduction
Micromorphic theory envisions a material body as a continuous collection of deformable particles; each possesses finite size and inner structure. On the other hand, classical continuum mechanics envisions a material body as a continuous collection of material points, each with infinitesimal size and no inner structure. The purpose of going beyond the classical continuum mechanics is to take into account the microstructure of the material body in question while still keeping the advantages of continuum theory intact. Because micromorphic theory can incorporate more complex substances, it can be utilized to treat blood, bubbly fluids, liquid crystals [1–5], etc. Incorporating the coupling between thermomechanics and electromagnetics into micromorphic theory can further enhance the understanding and the prediction of large classes of physical phenomena and provide innumerable technological applications. To cite a few, piezoelectric crystals, magnetic memory devices, energy generation by means of controlled plasma, and the study of many natural phenomena, such as earth magnetism and stellar dynamics. Phenomenologically important cross effects such as electostriction, magnetostriction, Peltier, Seebeck, Hall, Ettingshausen, Righi-Leduc, and Nernst effects can now be studied theoretically and numerically [6–8].
Microcontinuum field theories constitute extensions of the classical field theories concerned with deformations, motions, and electromagnetic interactions of material media as continua. It is emphasized that in the classical continuum theory, a point is represented by a geometrical point, infinitesimal in size and without inner structure. Then the question becomes How can one represent the intrinsic deformation of a point particle in a microcontinuum? Eringen settled this question by replacing the deformable point particle with a geometric point P and some vectors attached to P, which denote the orientations and intrinsic deformations of the material points in the deformable point particle [5,7,9]. Geometrically, a particle P is identified with its position vector X in a Lagrangian coordinate system and vectors Ξα(α=1,2,3,....,N) attached to P, representing the inner structure of P. Here N is the number of discrete material points in the particle. Now the motions may be expressed as
x=x(X,t)(1)
ξα=ξα(X,Ξα,t),α∈(1,2,3,...,N)(2)
where t is the time, ξα is the Eulerian coordinates of α-th material point in the particle. A medium with such general motions was named by Eringen as the microcontinuum of grade N. Obviously, this is too complicated. Therefore, in a so-called two-level continuum model, first let the position vector of a material point be decomposed as the sum of the position vector of the centroid of the particle and the position vector of a material point relative to the centroid:
x′=x+ξ,X′=X+Ξ(3)
Then let the motions be reduced to
x=x(X,t),ξ=ξ(X,Ξ,t)(4)
If the micromotions ξ=ξ(X,Ξ,t) is further reduced to an affine motion, i.e.,
ξ=χK(X,t)ΞKorξk=χkK(X,t)ΞK(5)
Then we arrive at the doorstep of the micromorphic theory.
The ultimate goal of this work is to formulate the finite element equations for a micromorphic electromagnetic continuum. To begin with, notice that the balance laws of micromorphic electromagnetic continuum consist of two parts: the thermomechanical (TM) part and the electromagnetic (EM) part. In Section 2, 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 [10], as well as in its special case, non-relativistic electromagnetic theory. In Section 3, the Eulerian description of the basic laws of micromorphic theory, including Conservation of Mass, Conservation of Microinertia, Balance of Linear Momentum, Balance of Moments of Momentum, Conservation of Energy, and Entropy Principle, is introduced. Also included in the basic laws of micromorphic theory are the electromagnetic (EM) part of the body force, body moment, and energy source [6,7,11–13]. It is emphasized that this inclusion establishes the link between the thermomechanical part (TM) and the electromagnetic part (EM) of the micromorphic electromagnetic theory. The constitutive theory for thermo-visco-elastic-plastic-electromagnetic (TVEP-EM) materials is formulated in Section 4. 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 the independent constitutive variables and, of course, one needs to supply a set of governing equations for the newly added internal variables [5,14]. The materially linear constitutive equations are formulated in Section 5. 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, scalar potential ϕ and vector potential A are introduced in Section 6. The scalar potential ϕ and the vector potential A are related to the electric field E and magnetic flux B [10] as
B≡∇×A(6)
E≡−∇ϕ−1c∂A∂t(7)
Eventually, the governing equations for ϕ and A, not E and B, are obtained in Section 6. In Section 7, for finite element formulation, we link the displacement U, micromotion χ, 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β, χkKβ, Tβ, ϕβ, and Akβ are obtained. It is noticed that the governing equations for displacements uand micromotions χ are mainly acoustic wave equations; the governing equations for temperature T are mainly diffusion equations; The governing equations for scalar potential φ and vector potential A are mainly wave equations with wave speed in the order of the speed of light. The speed of acoustic waves and the speed of light are not in the same order of magnitude. To solve the five sets of differential equations, Eqs. (172)–(176), simultaneously is numerically impossible and physically meaningless. In Section 8, we propose to solve the problem in two stages. To the electromagnetic waves in terms of the scalar and vector potentials, everything else seems to stand still. Therefore, in StageI, we solve Eqs. (172)–(174) dynamically and seek the static or nearly static solutions. In StageII, based on the static or nearly static solutions obtained in StageI, namely u,T, and χ, we calculate the forcing terms for the differential equations for ϕ and A. Discussions and conclusions are given in Section 9.
2 Maxwell’s Equations
The balance laws of micromorphic electromagnetic 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 [7,10,12] as
∇⋅D=qorDk,k=q(8)
∇×E+1c∂B∂t=0 or εijkEk,j+1c∂Bi∂t=0(9)
∇⋅B=0 or Bk,k=0(10)
∇×H−1c∂D∂t=1cJ or εijkHk,j−1c∂Di∂t=1cJi(11)
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(12)
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 [10]
E∗=γ(E+β×B)−γ2γ+1β(β⋅E)(13)
B∗=γ(B−β×E)−γ2γ+1β(β⋅B)(14)
where β≡vc,β=β⋅β,γ≡11−β2, and v is the velocity of co-moving frame RC relative to the fixed laboratory frame RG, and c is the speed of light in vacuum. Jackson [10] 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∗)(15)
B=γ(B∗+β×E∗)−γ2γ+1β(β⋅B∗)(16)
The proof of Jackson’s formula, Eqs. (13)–(16), and further discussion of the Lorentz Transformation are included in Appendix A.
Following the same pattern of Eqs. (13)–(16), one may write and prove the validity of
D∗=γ(D+β×H)−γ2γ+1β(β⋅D),D=γ(D∗−β×H∗)−γ2γ+1β(β⋅D∗)(17)
H∗=γ(H−β×D)−γ2γ+1β(β⋅H),H=γ(H∗+β×D∗)−γ2γ+1β(β⋅H∗)(18)
M∗=γ(M+β×P)−γ2γ+1β(β⋅M),M=γ(M∗−β×P∗)−γ2γ+1β(β⋅M∗)(19)
P∗=γ(P−β×M)−γ2γ+1β(β⋅P),P=γ(P∗+β×M∗)−γ2γ+1β(β⋅P∗)(20)
Recall that γ≡(1−β2)−12, the Taylor series expansion of γ about β=0 is obtained as
γ≈1+12β2+34β4+....(21)
Let’s keep the Taylor series at most to β2. Then, for non-relativistic electromagnetic theory, we have
E∗=E+β×B,E=E∗−β×B∗(22)
B∗=B−β×E,B=B∗+β×E∗(23)
D∗=D+β×H,D=D∗−β×H∗(24)
H∗=H−β×D,H=H∗+β×D∗(25)
P∗=P−β×M,P=P∗+β×M∗(26)
M∗=M+β×P,M=M∗−β×P∗(27)
Also, it is noticed that Eringen [6] gave a detailed discussion and derivation and showed that
q∗=q,J∗=J−qv(28)
3 Balance Laws of Micromorphic Theory
To begin with, let’s briefly derive the balance laws of micromorphic theory [7]. There are three parts:
(I) The concepts of mass and inertia require a finite volume of the material body. Before taking a limiting process of volume densities, consider a particle P having volume element ΔV in the reference state and its image p in the spatial frame at time t. The total mass of these particles is the sum of the masses of microelements, i.e.,
ρoΔV=∫ΔVρo′dV′,ρΔv=∫Δvρ′dv′(29)
where the primed quantities refer to microelements of P and p. The relative position vector Ξ and ξ of the microelements are taken with respect to the centroids of ΔV and Δv, so that
∫ΔVρo′ΞdV′=0,∫Δvρ′ξdv′=0(30)
However, the second moments of ρo′dV′ and ρ′dv′ do not vanish and they are given by
ρoIKLΔV=∫ΔVρo′ΞKΞLdV′,ρiklΔv=∫Δvρ′ξkξldv′(31)
We assume the mass of the microelement is conserved during the motion
ρo′dV′=ρ′dv′(32)
This implies that, in the limit as ΔV→0 and Δv→0,
ρodV=ρdv(33)
ikl=IKLχkKχlL,IKL=iklχ¯Kkχ¯Ll(34)
These two balance laws can be re-written as
ddt∫vρdv=0(35)
ddt∫vρiklχ¯Kkχ¯Lldv=0(36)
(II) In classical continuum mechanics the energy associated with stresses can be written as ∫Δatklvldak. Now in micromorphic theory, follow Eqs. (4) and (5), we write
∫Δatkl′(vl+ξ˙l)dak′=(tklvl+mklmυlm)Δak(37)
where the microgyration tensor υkl, moment stress tensor mklm, and body moment ρlkl are defined through
ξ˙k=υklξl(38)
mklmΔak=∫Δatkl′ξmdak′(39)
ρllmΔv=∫Δvρ′fl′ξmdv′(40)
The kinetic energy per unit mass can now be expressed as
K=12{vkvk+iklυmkυml}(41)
Therefore, for micromorphic theory, the conservation of energy can be expressed as
ddt∫vρ(e+K)dv=∮s(tklvl+mklmυlm−qk)dak+∫vρ(fkvk+lklυkl+h)dv(42)
Consider a time-dependent rigid motion of the reference frame
x~k=Qkl(t)xl+bk(t)(43)
where Qkl is an orthogonal tensor, i.e.,
QQT=QTQ=I,det(Q)=1(44)
Then the velocity, angular velocity, microinertia tensor, and its material time rate can be expressed as
v~k=Qklvl+Q˙klxl+b˙k(45)
Ωkl=Q˙kmQlm(46)
i~kl=QkmimpQlp(47)
di~kldt=QkmdimpdtQlp+2Ωkmiml(48)
Suppose that at time t, the body is brought back to the original orientation (i.e., Q=I, having only translational and angular velocities, i.e., Ω and b˙ are constants). Then we have
v~k=vk+Ωklxl+b˙k(49)
Ωkl=Q˙kl(50)
i~kl=ikl(51)
di~kldt=dikldt+2Ωkmiml(52)
and also, it is noticed that
υ~kl=Ωkl+υkl(53)
Ωkl=−Ωlk(54)
Under these transformations, the mass density ρ, the internal energy density e, Cauchy stress tkl, moment stress mklm, and heat flux qk are not affected since the motion is rigid. However, the body force fk and body moment lkl must be accommodated by corresponding linear and micro-accelerations (spin inertia), i.e.,
f~k−v~k=fk−vk,l~kl−σ~kl=lkl−σkl(55)
Also, notice that
K~=12{v~kv~k+i~klυ~mkυ~ml}(56)
Then subtracting the energy balance law in the x−frame from that in the x~−frame gives
ddt∫vρ(K~−K)dv=∮s{tkl(v~l−vl)+mklm(υ~lm−υlm}dak+∫vρ{f~kv~k−fkvk+l~klυ~kl−lklυkl}dv(57)
After lengthy but straightforward derivation one obtains
∫v{[Ωlmxm+b˙l]f¯l+Ωlml¯lm}dv−∫v{ρ^(K~−K)+i~kl(υmkΩml+12ΩmkΩml)}=0(58)
where
f¯l≡tkl,k+ρ(fl−v˙l),l¯lm≡mklm,k+tml+ρ(llm−σlm)(59)
For arbitrary and independent variations of b˙l and Ωkl=−Ωlk, the coefficients these quantities must vanish. Hence
f¯l=0⇒tkl,k+ρ(fl−v˙l)=0(60)
l¯lm=l¯ml⇒mklm,k+tml−slm+ρ(llm−σlm)=0(61)
where the introduction of a symmetric sml=slm is due to the antisymmetric Ωkl=−Ωlk, and it is named microstress.
The law of conservation of energy can simply be derived from
ddt∫vρ(e+K)dv=∮s(tklvl+mklmυlm−qk)dak+∫vρ(fkvk+lklυkl+h)dv(62)
where e is the internal energy density, qk is the heat flux and h is the energy source density. The conservation of energy in local form can be obtained as
ρe˙=mklmυlm,k+tkl(vl,k−υlk)+sklυlk−qk,k+ρh(63)
(III) The second law of thermodynamics can be derived as
ddt∫vρηdv≥−∮s1θqkdak+∫vρhθdv(64)
where η is the entropy density per unit mass, θ is the absolute temperature. The local form can now be written as
ρη˙+(qiθ),i−ρhθ≥0(65)
Now the Eulerian description of the basic laws of micromorphic thermomechanical (TM)-electromagnetic (EM) theory can be expressed as [7]:
Conservation of Mass
ρ˙+ρvk,k=0(66)
Conservation of Microinertia
dikldt=ikmυlm+ilmυkm(67)
Balance of Linear Momentum
tji,j+ρ(fi−v˙i)+Fi=0(68)
Balance of Moments of Momentum
mklm,k+tml−sml+ρ(llm−σlm)+Llm=0(69)
Conservation of Energy
ρe˙=mklmυlm,k+tkl(vl,k−υlk)+sklυkl−qk,k+ρh+W(70)
Entropy Principle (Second Law of Thermodynamics)
ρη˙+(qiθ),i−ρhθ≥0(71)
where ρ,ikl=ilk,tkl≠tlk,skl=slk,mklm,e,qk,η,θ,vi,υlm,fi,Fi,llm,Llm,h, and W are the mass density, microinertia, Cauchy stress, microstress, moment stress, internal energy density, hear flux, entropy density, absolute temperature, velocity, microgyration, body force (TM), body force (EM), body moment (TM), body moment (EM), energy source (TM), and energy source (EM), respectively; the spin inertia σ is defined as
σkl≡iml(υ˙km+υknυnm)(72)
The EM part of the body force, body moment, and energy source are given as [6,7,11–13]
F=qE+(∇E)⋅P+(∇B)⋅M+1c{J×B+∇⋅(vP×B)+∂∂t(P×B)}(73)
L=E∗⊗P+B⊗M∗(74)
W=E∗⋅(P˙+P∇⋅v)−M∗⋅B+J∗⋅E∗(75)
Or they can be expressed as
Fk=qEk+El,kPl+Bl,kMl+1cεklm{JlBm+(vnPlBm),n+∂∂t(PlBm)}(76)
Lkl=Ek∗Pl+BkMl∗(77)
W=Ek∗(P˙k+Pkvl,l)−Mk∗B˙k+Jk∗Ek∗(78)
Notice that, if A is a function of Eulerian coordinate xi and time t, then A˙=dAdt=∂A∂t+A,ivi.
Define the Helmholtz free energy density as
ψ≡e−θη−ρ−1Ek∗Pk(79)
Notice that the Helmholtz free energy density defined this way is quite different from its usual way. Then one obtains the Lagrangian description of the law of conservation of energy and the Clausius-Duhem (CD) inequality as (for detailed derivation, see Appendix B)
ρoe˙=MKLMγ˙KLM+TKLmα˙KL+SKLmβ˙KL−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗+ρoh(80)
−ρo(ψ˙+ηθ˙)+MKLMγ˙KLM+TKLmα˙KL+SKLmβ˙KL−1θQKθ,K−PKE˙K∗−MK∗B˙K+JK∗EK∗≥0(81)
where the mechanical part of the Cauchy stress and the microstress are defined as
tklm≡tkl+PkEl∗+Mk∗Bl,TKLm≡JtklmXK,kχlL(82)
sklm≡skl+12(PkEl∗+PlEk∗+Mk∗Bl+Ml∗Bk),SKLm≡12Jsklmχ¯Kkχ¯Ll(83)
In other words, the electromagnetic parts are defined as
tklem≡−PkEl∗−Mk∗Blsklem≡−12(PkEl∗+PlEk∗+Mk∗Bl+Ml∗Bk)(84)
If there is no coupling with electromagnetics, i.e., P=E∗=M∗=B=0, then one has
tklm=tkl,tklem=0sklm=skl,sklem=0(85)
4 Constitutive Equations
In this section we are going to formulate the micromorphic 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 [14]. For micromorphic thermo-visco-elastic-plastic continuum, a set of internal variables is introduced as [5]
W={αP,βP,γP,R}(86)
where αp,βp,γp are respectively the plastic strains corresponding to the generalized Lagrangian strains α,β,γ defined as
αKL≡xk,Kχ¯Lk−δKLβKL≡χkKχkL−δKL=βLKγKLM≡χ¯KkχkL,M(87)
and R, named as the hardening parameters, is a generalized vector of internal variables.
Remark 1:
In this work, the set of internal variables W includes the Lagrangian strains in micromorphic theory and generalized vector of hardening parameters R. So far, it might be considered as the most general form for internal variables.
Now let the formulation of the micromorphic constitutive theory for the TVEP-EM materials begin with
Z=Z(Y)(88)
where
Z={Tm,Sm,M,Q,ψ,η,P,M∗,J∗}={TKLm,SKLm,MKLM,QK,ψ,η,PK,MK∗,JK∗}(89)
Y={α,β,γ,α˙,β˙,γ˙,θ,∇θ,E∗,B,W}={α,β,γ,α˙,β˙,γ˙,θ,∇θ,E∗,B,αp,βp,γp,R}={αKL,βKL,γKLM,α˙KL,β˙KL,γ˙KLM,θ,θ,K,EK∗,BK,αKLp,βKLp,γKLMp,R}(90)
The definitions of EK∗,BK,PK,MK∗, and JK∗ are given in Appendix B.
To separate the material behavior into two distinct parts: thermo-visco-elastic-electromagnetic (TVE-EM) part and thermo-visco-elastic-plastic-electromagnetic (TVEP-EM) part, a scalar-valued yield function is introduced as
f=f(Y)(91)
For a set of fixed values of W, a hyper surface, named yield surface, is determined by
f(Y)=0(92)
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∗,B, i.e.,
λ=∂f∂αKLα˙KL+∂f∂βKLβ˙KL+∂f∂γKLMγ˙KLM+∂f∂θθ˙+∂f∂α˙KLα¨KL+∂f∂β˙KLβ¨KL+∂f∂γ˙KLMγ¨KLM+∂f∂θ,Kθ˙,K+∂f∂EK∗E˙K∗+∂f∂BKB˙K(93)
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. Following the axiom of equipresence, the governing equations for the internal variables of plasticity W are proposed to be
W˙=λ∗Φ(Y)(94)
where
λ∗={0iff<0orλ<0λiff=0andλ≥0(95)
Now we substitute Eqs. (88) and (94) into the Clausius-Duhem inequality, Eq. (81). It leads to
−ρo{∂ψ∂αKLα˙KL+∂ψ∂βKLβ˙KL+∂ψ∂γKLMγ˙KLM+∂ψ∂α˙KLα¨KL+∂ψ∂β˙KLβ¨KL+∂ψ∂γ˙KLMγ¨KLM+∂ψ∂θθ˙+∂ψ∂θ,Kθ˙,K+∂ψ∂EK∗E˙K∗+∂ψ∂BKB˙K+ηθ˙+∂ψ∂W⋅λ∗Φ}+MKLMγ˙KLM+TKLmα˙KL+SKLmβ˙KL−1θQKθ,K−PKE˙K∗−MK∗B˙K+JK∗EK∗≥0(96)
Since this inequality is linear in α¨KL,β¨KL,γ¨KLM,E˙K∗,B˙K,θ˙,θ˙,K, it leads to
ψ=ψ(αKL,βKL,γKLM,θ,EK∗,BK,αKLp,βKLp,γKLMp,R)(97)
η=−∂ψ∂θ(98)
TKLme≡ρo∂ψ∂αKL,TKLm=TKLme+TKLmd(Y)(99)
SKLme≡ρo∂ψ∂βKL,SKLm=SKLme+SKLmd(Y)(100)
MKLMe≡ρo∂ψ∂γKLM,MKLM=MKLMe+MKLMd(Y)(101)
PK=−ρo∂ψ∂EK∗(102)
MK∗=−ρo∂ψ∂BK(103)
MKLMdγ˙KLM+TKLmdα˙KL+SKLmdβ˙KL−1θQKθ,K+JK∗EK∗−ρo∂ψ∂W⋅λ∗Φ≥0(104)
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. (104), must be satisfied for any value of the loading rate λ, it implies
∂ψ∂W⋅Φ≤0(105)
MKLMdγ˙KLM+TKLmdα˙KL+SKLmdβ˙KL−1θQKθ,K+JK∗EK∗≥0(106)
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 [15]. Therefore, in the case of loading, we have f˙=0, i.e.,
f˙=λ∗+∂f∂W⋅W˙=λ∗+λ∗∂f∂W⋅Φ=λ∗(1+∂f∂W⋅Φ)=0(107)
This gives another constitutive constraint to the plasticity
∂f∂W⋅Φ=−1(108)
5 Linear Constitutive Equations
For micromorphic TVEP-EM solid, define the potential energy density function Σ as
Σ≡ρoψ=Σ(α,β,γ,θ,E∗,B,αp,βp,γp,R)(109)
Following the idea of Green and Naghdi [16], and assuming the hardening parameters R do not appear in the potential energy density function, Eq. (109) can now be written as
Σ=Σ(α−αp,β−βp,γ−γp,T,E∗,B)≡Σ(αe,βe,γe,T,E∗,B)=Σo−ρoηoT−12ρoγT2To−A¯KLαKLeT−B¯KLβKLeT−C¯KLMγKLMeT+aKTEK∗+bKTBK+{12AKLMNαKLeαMNe+12BKLMNβKLeβMNe+12CKLMNPQγKLMeγNPQe+DKLMNαKLeβMNe+EKLMNPαKLeγMNPe+FKLMNPβKLeγMNPe}+FKLM1αKLeEM∗+FKLM2αKLeBM+GKLM1βKLeEM∗+GKLM2βKLeBM−HKLMN1γKLMeEN∗−HKLMN2γKLMeBN−12D¯KL1EK∗EL∗−12D¯KL2BKBL−D¯KL3EK∗BL(110)
where To is the reference temperature and T is the temperature variation, and
T≡θ−To,To>0,|T|<<To(111)
Remark 2:
Green and Naghdi’s idea is to replace α,β,γ,αp,βp,γp by αe,βe, and γe, which are defined as αe≡α−αp, βe≡β−βp,γe≡γ−γp. This simplifies a lot. Also, it shows the physical meaning clearly.
From now on we assume that the material is isotropic. Therefore, all the odd order material property tensors are vanishing. Therefore, Eq. (110) can be re-written as
Σ=Σ(α−αp,β−βp,γ−γp,T,E∗,B)≡Σ(αe,βe,γe,T,E∗,B)=Σo−ρoηoT−12ρoγT2To−A¯KLαKLeT−B¯KLβKLeT+{12AKLMNαKLeαMNe+12BKLMNβKLeβMNe+DKLMNαKLeβMNe+12CKLMNPQγKLMeγNPQe}−HKLMN1γKLMeEN∗−HKLMN2γKLMeBN−12D¯KL1EK∗EL∗−12D¯KL2BKBL−D¯KL3EK∗BL(112)
Then it is straightforward to obtain
η=−∂ψ∂θ=−∂ψ∂T=ηo+γTTo+{A¯KLαKLe+B¯KLβKLe}/ρo(113)
TKLme=ρo∂ψ∂αKLe=∂Σ∂αKLe=−A¯KLT+AKLMNαMNe+DKLMNβMNe(114)
SKLme=ρo∂ψ∂βKLe=∂Σ∂βKLe=−B¯KLT+BKLMNβMNe+DMNKLαMNe(115)
MKLM=ρo∂ψ∂γKLMe=CKLMNPQγNPQe−HKLMN1EN∗−HKLMN2BN(116)
PK=−ρo∂ψ∂EK∗=HLMNK1γLMNe+D¯KL1EL∗+D¯KL3BL≡PKm+D¯KL1EL∗+D¯KL3BL(117)
MK∗=−ρo∂ψ∂BK=HLMNK2γLMNe+D¯KL2BL+D¯LK3EL∗≡MK∗m+D¯KL2BL+D¯LK3EL∗(118)
Moreover, the 2nd order, 4th order, and 6th order material property tensors are made of a few material constants, for example,
A¯KL=A¯δKL,B¯KL=B¯δKLAKLMN=λ1δKLδMN+μ11δKMδLN+μ12δKNδLMBKLMN=λ2δKLδMN+μ2(δKMδLN+δKNδLM)DKLMN=λ3δKLδMN+μ3(δKMδLN+δKNδLM)CKLMNPQ=C1δKLδMNδPQ+C2δKLδMPδNQ+C3δKLδMQδNP+C4δKNδLMδPQ+C5δKMδNLδPQ+C6δKMδLPδNQ+C7δKNδLPδMQ+C8δKPδLQδMN+C9δKNδLQδMP+C10δKPδNLδMQ+C11δKQδLPδMN+C12δKQδLMδNP+C13δKMδLQδNP+C14δKPδLMδNQ+C15δKQδNLδMP(119)
Onsager’s Postulate
Now we follow Onsager’s Postulate [17,18] and construct a potential of dissipation quadratic as
D=12dKLMN1α˙KLα˙MN+dKLMN2α˙KLβ˙MN+12dKLMN3β˙KLβ˙MN+12fIJKLMN1γ˙IJKγ˙LMN+fIJKM2γ˙IJKθ,M+fIJKM3γ˙IJKEM∗+12gMN1θ,Mθ,N+gMN2θ,MEN∗+12gMN3EM∗EN∗(120)
where α˙,β˙,γ˙,∇θ, and E∗ are the thermodynamic forces. Correspondingly, the thermodynamic fluxes Tmd,Smd,Md,Q, and J∗ can be obtained as
TKLmd=∂D∂α˙KL=dKLMN1α˙MN+dKLMN2β˙MN(121)
SKLmd=∂D∂β˙KL=dMNKL2α˙MN+dKLMN3β˙MN(122)
MIJKd=∂D∂γ˙IJK=fIJKLMN1γ˙LMN+fIJKM2θ,M+fIJKM3EM∗(123)
−1θQM=∂D∂θ,M=fIJKM2γ˙IJK+gMN1θ,N+gMN2EN∗(124)
JM∗=∂D∂EM∗=fIJKM3γ˙IJK+gNM2θ,N+gMN3EN∗(125)
One may verify that the CD inequality Md⋮γ˙+Tmd:α˙+Smd:β˙−1θQ⋅∇θ+J∗⋅E∗≥0 is now reduced to
D≥0(126)
which is precisely what the Onsager’s Postulate means.
Because the material is assumed to be isotropic, similar to Eq. (119), one could have
gKL1=g¯1δKL,gKL2=g¯2δKL,gKL3=g¯3δKLdKLMN1=λ¯1δKLδMN+μ¯11δKMδLN+μ¯12δKNδLMdKLMN2=λ¯2δKLδMN+μ¯2(δKMδLN+δKNδLM)dKLMN3=λ¯3δKLδMN+μ¯3(δKMδLN+δKNδLM)fKLMN2=κ2δKLδMN+τ21δKMδLN+τ22δKNδLMfKLMN3=κ3δKLδMN+τ31δKMδLN+τ32δKNδLMfKLMNPQ1=f1δKLδMNδPQ+f2δKLδMPδNQ+f3δKLδMQδNP+f4δKNδLMδPQ+f5δKMδNLδPQ+f6δKMδLPδNQ+f7δKNδLPδMQ+f8δKPδLQδMN+f9δKNδLQδMP+f10δKPδNLδMQ+f11δKQδLPδMN+f12δKQδLMδNP+f13δKMδLQδNP+f14δKPδLMδNQ+f15δKQδNLδMP(127)
Return Mapping Algorithm for Micromorphic Plasticity
Define the generalized total strain tensor, elastic strain tensor, and plastic strain tensor as
E≡|αβγ|,Ee≡|αeβeγe|,Ep≡|αpβpγp|(128)
Notice that
αe≡α−αp,βe≡β−βp,γe≡γ−γp(129)
Therefore
E=Ee+Ep(130)
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 [5,19]
f=‖E−Ep−G‖−(EY+κHg)=‖Ee−G‖−(EY+κHg)(131)
where EY,H, and κ are material constants.
Remark 3:
This is a special form. One may want to improve it or to generalize it.
Let the evolution equations for the plastic strain tensor and the hardening parameters be
E˙p=Ωξ‖ξ‖,G˙=ΩC1(1−κ)Hξ‖ξ‖,g˙=ΩC2(132)
where
ξ≡E−Ep−G=Ee−G,‖ξ‖≡{ξ⋅ξ}12(133)
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)(134)
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)(135)
It is worthwhile to mention that in this case (f~n+1≤0,t=tn+1), 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‖(136)
Gn+1=Gn+ΔΩC1(1−κ)Hξ~n+1/‖ξ~n+1‖(137)
gn+1=gn+ΔΩC2(138)
where
ΔΩ=f~n+11+C1H−(C1−C2)κH(139)
Now one can prove that yield function, based on the updated En+1p,Gn+1, and gn+1, returns to zero [5] i.e.,
fn+1=‖E−En+1p−Gn+1‖−(EY+κHgn+1)=0(140)
Remark 4:
The yield function based on the updated values returning to zero means the return mapping algorithm works. No need to worry theoretically.
From Eqs. (114)–(118), it is noticed that the elastic parts of the stresses, TKLme,SKLme,MKLMe, and PK,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 well-known that the electric field E and the magnetic flux B are not independent of each other, nor are the dielectric displacement D and the magnetic field H. Mathematically, this is noticed from Maxwell’s equations, especially, Eqs. (9) and (11). To resolve this problem, especially in the finite element formulation presented in the next section, introduce the scalar potential ϕ and the vector potential A, which are related to the electric field E and magnetic flux B as
B≡∇×AorBi=eijkAk,j(141)
E≡−∇ϕ−1c∂A∂torEi≡−ϕ,i−1c∂Ai∂t(142)
Now we follow the same pattern as Eqs. (117) and (118) and write the polarization Pk and magnetization Mk as
Pk=Pkm+D1Ek+D3BkMk=Mmm+D2Bk+D3Ek(143)
where Pkm and Mkm are the polarization and magnetization induced by strains, i.e.,
Pkm≡Hlmnk1γlmneMkm≡Hlmnk2γlmne(144)
Admittedly the introduction of Eqs. (143) and (144) involves approximations. Without making these approximations, it is difficult, if not impossible, to derive the governing equations for scalar potential ϕ and vector potential A.
Now the dielectric displacement vector D and the magnetic field vector H can be expressed as
Dk=Ek+Pk =Ek+Pkm+D1Ek+D3Bk =(1+D1)Ek+D3Bk+Pkm ≡εEk+D3Bk+Pkm(145)
Hk=Bk−Mkm =Bk−(Mkm+D2Bk+D3Ek) =(1−D2)Bk−D3Ek−Mkm ≡μ−1Bk−D3Ek−Mkm(146)
Substituting Eqs. (141) and (142) into the Maxwell’s equations leads to
Dk,k=εEk,k+D3Bk,k+Pk,km=−ε{ϕ,kk+1cA˙k,k}+Pk,km=q(147)
∇×E+1c∂B∂t=∇×{−∇ϕ−1cA˙}+1c{∇×A˙}=0(148)
∇⋅B=∇⋅(∇×A)=0(149)
∇×H−1c∂D∂t=∇×[μ−1B−D3E−M∗m]−1c{εE˙+D3B˙+P˙m}=μ−1∇×B−∇×M∗m−εcE˙−1cP˙m=μ−1∇×∇×A+εc[∇ϕ˙+1cA¨]−1cP˙m−∇×Mm=1cJ(150)
It is seen that two of the four Maxwell’s equations, Eqs. (148) and (149), are identically satisfied; the other two equations, Eqs. (147) and (150) can be re-written as
ϕ,kk+1cA˙k,k={Pk,km−q}/ε(151)
μ−1[Aj,ij−Ai,jj]+εcϕ˙,i+εc2A¨i=1c[Ji+P˙im]+eijkMk,jm(152)
Gauge Transformation
The famous gauge transformation is represented by
ϕ′=ϕ−1cψ˙A′=A+∇ψ(153)
Substituting Eq. (153) into Eqs. (151) and (152) results
∇2ϕ′+1c∇⋅A˙′=ϕ,ii−1cψ˙,ii+1c{∇⋅A˙+∇2ψ˙}=ϕ,ii+1c∇⋅A˙(154)
−μ−1Ai,jj′+μ−1Aj,ji′+ε1cϕ˙′,i+1c2εA¨′i=μ−1{−Ai,jj−ψ,ijj+Aj,ji+ψ,jji}+ε{1c[ϕ˙,i−1cψ¨,i]+1c2[A¨i+ψ¨,i]}=−μ−1{Ai,jj−Aj,ji}+ε{1cϕ˙,i+1c2A¨i}=1c[Ji+P˙im]+εijkMk,jm(155)
In other words, the function ψ is arbitrary. We utilize this freedom to arrive at the Lorentz Condition.
Lorentz Condition
ε1cϕ˙+μ−1∇⋅A=0(156)
Substituting the Lorentz condition into Eqs. (147) and (150) results
−ε{∇2ϕ+1c∇⋅A˙}=−ε{∇2ϕ−εμ1c2ϕ¨}=q−∇⋅Pm(157)
−μ−1∇2Ai+1c2εA¨i=Jic+εijkMk,jm+1c∂Pim∂t(158)
In other words, the Maxwell’s equations are reduced to
εμ1c2ϕ¨−∇2ϕ=a(159)
εμ1c2A¨i−∇2Ai=Ci(160)
where
a≡1ε(q−Pk,km)(161)
Ck≡μ{Jkc+εkijMj,im+1cP˙km}(162)
This means the Maxwell’s equations become two wave equations with forcing terms and the speed of wave propagation is equal to
v=c1εμ(163)
Since both the material constants ε and η are positive and greater than one, therefore Eq. (163) says the speed of EM wave propagation in material is smaller than that in vacuum or in air.
7 Finite Element Formulation
Since the Lagrangian stresses and the Eulerian stresses are related as
TKL=JtklXK,kχlL,tkl=J−1TKLxk,Kχ¯LlSKL=12Jsklχ¯Kkχ¯Ll,skl=2J−1SKLχkKχlLMKLM=JmmklXM,mχkKχ¯Ll,mmkl=J−1MKLMxm,Mχ¯KkχlL(164)
One may verify that the balance law of linear momentum and balance law of moment of momentum, can be written as
(TKLχ¯Ll),K+ρo(fl−v˙l)+JFl=0(165)
(MLMKχ¯LlχmM),K+TMLxm,Mχ¯Ll−2SLMχlLχmM+ρo(llm−σlm)+JLlm=0(166)
Recall that the conservation of energy can be expressed as
ρoe˙=MKLMγ˙KLM+TKLmα˙KL+S¯KLmβ˙KL−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗+ρoh(167)
Based on the constitutive equations, Eqs. (113)–(118) and (120)–(134), Eq. (167) can be further derived to be [5]
ρoγT+ToToT˙+(T+To)(A¯KLα˙KLe+B¯KLβ˙KLe)=TKLmeα˙KLp+SKLmeβ˙KLp+MKLMeγ˙KLMp+TKLdα˙KL+SKLdβ˙KL+MKLMdγ˙KLM−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗+ρoh(168)
Recall the governing equations of the scalar potential ϕ and the vector potential A as
εμc2∂2ϕ∂t2−ϕ,ll=a(169)
εμc2∂2Ak∂t2−Ak,ll=Ck(170)
We now link the displacement, micromotion, 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αχkK=NαχkKαT=NαTαϕ=NαϕαAk=NαAkα(171)
where Nα[α=1,2,3,....,8] are the 8 shape functions for a standard 8-noe 3D element. The finite element equations are obtained as (for detailed derivation, see Appendix C)
MαβU¨Mβ=fMα1+fMα2+fMα3(172)
IMαβχ¨kKβ=ϕkKα1+ϕkKα2+ϕkKα3+ϕkKα4+ϕkKα5+ϕkKα6(173)
γMαβT˙β=qα1+qα2+qα3+qα4+qα5+qα6+qα7(174)
mαβϕ¨β+kαβϕβ=ξ^α+a~α(175)
mαβA¨kβ+kαβAkβ=A^kα+C^kα(176)
where
Mαβ≡∫VρoNαNβdV=Mβα(177)
IMαβ≡∫VρoINαNβdV=IMβα(178)
γMαβ≡∫VρoγNαNβdV=γMβα(179)
mαβ≡∫v1c¯2NαNβdv(180)
kαβ≡∫vNα,iNβ,idv(181)
fMα1≡∫SFF^lδlMNαdS≡∫SFF^MNαdSfMα2≡−∫VTKLχ¯LlδlMNα,KdVfMα3≡∫V(ρofl+JFl)δlMNαdV(182)
ϕkKα1≡∫SMM^kKNαdSϕkKα2≡−∫VMLMPχ¯LkχmMχ¯Km,PNαdVϕkKα3≡−∫VMLKMχ¯LkNα,MdVϕkKα4≡∫VTMLxm,Mχ¯Lkχ¯KmNαdVϕkKα5≡∫V−2SLKχkLNαdVϕkKα6≡∫V(ρolkm+JLkm)χ¯KmNαdV(183)
qα1≡∫V−(T+To)[A¯KLα˙KLe+B¯KLβ˙KLe]NαdVqα2≡∫V{TKLmeα˙KLp+SKLmeβ˙KLp+MKLMeγ˙KLMp}NαdVqα3≡∫V{TKLdα˙KL+SKLdβ˙KL+MKLMdγ˙KLM}NαdVqα4≡∫V{−PKE˙K∗−MK∗B˙K+JK∗EK∗}NαdVqα5≡−∫SqQ^NαdSqα6≡∫VQKNα,KdVqα7≡∫VρohNαdV(184)
ξ^α=∫sξϕ,iniNαds≡∫sξξNαdsa^α≡∫vaNαdv(185)
A^kα=∫sςAk,iniNαds≡∫sςA~kNαdsC^kα≡∫vCkNαdv(186)
In this work, for simplicity, we assume that the microinertia is isotropic, i.e.,
IKL=IδKL(187)
which implies
ikl=IKLχkKχlL=IδKLχkKχlL=IχkKχlK(188)
Also, it is worthwhile to note that on the boundary the surface force, the surface moment and the heat flow are specified respectively as
F^l=TKLχ¯LlnK,M^kK=MLKMχ¯LknM,Q^=QKnK(189)
where nKis the unit outward normal to the boundary surface.
8 Two-Stage Solutions
It is noticed that Eqs. (165)–(170) are governing equations for u,χ,T,ϕ, and A, which are coupled together—that is why it is coined with the name micromorphic thermomechanical–electromagnetic coupling. To solve the five sets of differential equations, Eqs. (165)–(170), simultaneously, is numerically impossible and physically meaningless. Also, relative to the EM wave, everything else is standing still. The strategy of a two-stage solution, in general, can be described as follows.
Consider that there are two sets of variables, Set A and Set B. Based upon their general governing equations, the variables in Set A move faster and the variables in Set B move slower. Then, in Stage I, we seek the static solutions of the variables in Set B. Once the solutions in Set B are obtained, one may find the forcing terms for those variables in Set A. Then, for Stage II, we solve for the dynamic solutions for variables in Set A.
In this work, in StageI, we are seeking the static solutions of a hollow cylinder subjected to twisting. The finite element results (T,χavg,Pavg,Tavg) of StageI are plotted on the deformed shape (cf. Fig. 1). Since in micromorphic theory the Cauchy stress tensor and the corresponding 2nd order Piola-Kirchhoff stress are not symmetric, the L2norm of the 2nd order Piola-Kirchhoff stress is defined as
Tavg≡‖T‖=T112+T222+T332+T232+T322+T312+T132+T122+T212(190)
Figure 1: (a) Temperature variation, (b) L2−norm of the micromotion, (c) L2−norm of the plastic strain, (d) L2−norm of the 2nd order Piola-Kirchhoff stress. The color bars indicate the magnitudes of those variables.
The L2−norm of the plastic strain αp is defined as
Pavg≡{αKLpαKLp}1/2(191)
Similarly, one may define the L2−norms of βp and γp as
Pavg2≡{βKLpβKLp}1/2,Pavg3≡{γKLMpγKLMp}1/2(192)
The L2−norm of micromotion is defined as
χavg≡‖χ−δ‖=(χ11−δ11)2+(χ22−δ22)2+(χ33−δ33)2+χ232+χ322+χ312+χ132+χ122+χ212(193)
This means if there is no micromotion, then χkK=δkK and δkK is the shifter. Therefore in Eq. (193) we have three terms like χ11−δ11,χ22−δ22, and χ33−δ33 and the other six terms like χ23,χ32,χ31,χ13,χ12,χ21.
In this work, we are using the MKS system, i.e., meter, kilogram, second are the units for length, mass, and time, respectively. Also, let the units of an electric charge, electric field, and temperature be q, e, and oK, respectively. Some typical numerical values of the material constants common to both stages are specified as
A¯=0.002K/M/s2/oK,B¯=0.002K/M/s2/oKλ1=300K/M/s2,μ11=225K/M/s2,μ12=75K/M/s2λ2=3K/M/s2,μ2=3K/M/s2λ3=0.75K/M/s2,μ3=0.75K/M/s2Ci=10KM/s2,i∈[1,15]λ¯1=0.06KM/s,μ¯11=0.045KM/s,μ¯12=0.015KM/sλ¯2=0.0006KM/s,μ¯2=0.0004KM/sλ¯3=0.00015KM/s,μ¯3=0.0001KM/sg1=2000KM/s3/(oK)2g2=10K/(oK)/s3/eg3=10K/M/s3/e2EY=1,κ=0.2,H=5,C1=0.6,C2=0.4γ=200000M2/s2/oK,To=300oKρo=1000K/M3,I=0.002M2ε=2,μ=2(194)
These numerical values are our trial values, to the most are our educated guessing. Because, to the best of our knowledge, we don’t know whether there exits a testing machine which can measure microstresses, moment stresses, and micromotions.
In StageI, we solve Eqs. (165), (166) and (168) dynamically and seek the static or nearly static solutions with E∗=B=P=M∗=0. The assumption E∗=B=P=M∗=0 implies that F=L=W=0 (cf. Eqs. (36)–(38)) and also tm=t,sm=s,Tm=T,Sm=S (cf. Eqs. (42)–(44)). In StageII, based on the static or nearly static solutions obtained in StageI, namely U,T, and χ, we calculate the forcing terms, a and C, and solve the differential equations, Eqs. (169) and (170), for ϕ and A. Of course, one may generalize the procedure a little bit: assume E,B,P,M are constants, not necessarily vanishing. This mean F,L, and W are not vanishing either. They should be calculated according to Eqs. (36)–(38), respectively.
In a simpler case with the assumption that there are no free charge and no current, the forcing terms are reduced to
a=−1εPk,km(195)
Ck=μ{εkijMj,im+1cP˙km}(196)
Also, because Pkm and Mkm are the polarization and magnetization induced by strains (cf. Eq. (144)) and from StageI the strains are static or nearly static solutions, Eq. (196) is further reduced to
Ck=μεkijMj,im(197)
What if we would like to interchange the order of treatment for Set A and Set B. Of course, the variables in Set A move fast, but if we seek static solutions for those variables in Set A and seek dynamic solutions for those variables in Set B, then the order is reversed. For example, now in Stage I, we seek static solutions of scalar potential and vector potentials—they are mutually independent of each other. Then, in Stage II, we seek the dynamic solutions of temperature, displacements, and micromotions. Of course, they are moving faster than the static scalar and vector potential—their velocities are zero now.
Numerical Results of the Torsion Problem
Let the specimen be a hollow cylinder occupying
V={R,θ,Z|1.0≤R≤2.0,0≤θ≤2π,0≤Z≤30}(198)
The finite element mesh has 2520 nodes and 1920 8-node solid element. In StageI, each node has 13 unknowns (three for U, nine for χ, and one for T), therefore this sample problem has 32,760 degrees of freedom. Here central difference method is employed to solve Eqs. (165)–(168) which are two wave equations and one diffusion equation. To ensure the numerical stability, the time step Δt cannot be too large. For this sample problem, in StageI, we use Δt=0.001sec and have verified that such Δt not only provides stability but also accuracy.
The boundary conditions are specified as: the surface force, surface moment, and heat flow are vanishing (cf. Eq. (189)), i.e., F^M=M^kK=Q^=0, except
(I) At the bottom surface, Z=0, the displacements and temperature are specified as
UX=UY=UZ=T=0(199)
(II) At the top surface, Z=30, a torsion and an elongation are applied, and temperature is specified as
Ux=rcos(θ+Δθ)−X,Uy=rsin(θ+Δθ)−Y,Uz=3.0,T=5.0(200)
where Δθ=0.25π and
r≡X2+Y2,θ=tan−1(Y/X)(201)
The corresponding system of equations, Eqs. (165)–(168), statically can be expressed as
fMα1+fMα2+fMα3=0orf(U,χ,T)=0(202)
ϕkKα1+ϕkKα2+ϕkKα3+ϕkKα4+ϕkKα5+ϕkKα6=0orφ(U,χ,T)=0(203)
qα1+qα2+qα3+qα4+qα5+qα6+qα7=0orq(U,χ,T)=0(204)
Usually we would solve Eqs. (202)–(204) by using Newton-Raphson method. Since, in the finite element analysis of micromorphic theory each node has 13 degrees of freedom, this simple model has 32,760 degrees of freedom. In this case, Newton-Raphson method involves a stiffness matrix of size 32,760 × 32,760. This means it is practically impossible to solve it by Newton-Raphson method. Another way is to solve Eqs. (165)–(168) dynamically with proper damping coefficients so that, with reasonable number of time steps, one may obtain nearly static solutions.
The solutions for StageI are presented in Fig. 1.
Fig. 1 are the nearly static solutions of temperature variation, L2-norms of micromotion, plastic strain, and 2nd order Piola-Kirchhoff stress. One may observe the axis symmetry, the twist of the hallow cylinder through the finite element mesh, and the magnitude through the color bars.
In StageII, we investigate four cases. In general, we have
ϕ=Ax=Ay=Az=0(205)
Except for the following four cases:
Case 1: εμc2A¨1−∇2A1=C1(206)
Initial Conditions: A1=A˙1=0
Boundary Conditions: at Z=0 A1={−sin(ωt)if0≤ωt≤4π0ifωt≥4π
At Z=30, A1=0
Case2:εμc2A¨2−∇2A2=C2(207)
Initial Conditions: A2=A˙2=0
Boundary Conditions: at Z=0 A2={−sin(ωt)if0≤ωt≤4π0ifωt≥4π
At Z=30, A2=0
Case 3:εμc2A¨3−∇2A3=C3(208)
Initial Conditions: A3=A˙3=0
Boundary Conditions: at Z=0 A3={−sin(ωt)if0≤ωt≤4π0ifωt≥4π
At Z=30, A3=0
Case 4:εμc2ϕ¨−∇2ϕ=a(209)
Initial Conditions: ϕ=ϕ˙=0
Boundary Conditions: at Z=0 ϕ={−sin(ωt)if0≤ωt≤4π0ifωt≥4π
At Z=30, ϕ=0
Notice that, in Eqs. (206)–(209), ω=4π/(400Δt).
It is worthwhile to mention that the time step in StageII is Δt=0.25×10−9s, comparing with Δt=0.001 s, in StageI. This huge difference justifies the treatment of the coupling problem by a two-stage solution procedure. Also, the boundary conditions, associated with Eqs. (206)–(209), show that there is a half sine wave started at the bottom surface (Z=0) and then bounced back and forth between the top and bottom surfaces (cf. Figs. 2–5)—this is indicated through boundary conditions at Z=0 after ωt≥4π and at Z=30.
Figure 2: Numerical results of Ax for Case 1: The first, second, third, fourth, and fifth column are referred to time at 400Δt,800Δt,1200Δt,1600Δt,and2000Δt, respectively. The color bars indicate the magnitudes of the variables.
Figure 3: Numerical results of Ay for Case 2: The first, second, third, fourth, and fifth column are referred to time at 400Δt,800Δt,1200Δt,1600Δt,and2000Δt, respectively. The color bars indicate the magnitudes of the variables.
Figure 4: Numerical results of Az for Case 3: The first, second, third, fourth, and fifth column are referred to time at 400Δt,800Δt,1200Δt,1600Δt and 2000Δt, respectively. The color bars indicate the magnitudes of the variables.
Figure 5: Numerical results of ϕ for Case 4: The first, second, third, fourth, and fifth column are referred to time at 400Δt,800Δt,1200Δt,1600Δt, and 2000Δt, respectively. The color bars indicate the magnitudes of the variables.
To understand the StageII solutions about electromagnetic wave propagation in deformed micromorphic medium, we present the following discussion. For example, in Case 1, additional to the forcing term (cf. Eqs. (161) and (162)), and boundary condition, Eq. (206), one will obtain ϕ=Ay=Az=0 for all the time and a wave propagation of Ax which will induce a wave propagation of E and B according to B=∇×A,E=−∇ϕ−A˙/c. In other words, we obtain
B=|ijk∂∂x∂∂y∂∂zAx00|=∂Ax∂zj−∂Ax∂ykEx=−1cdAxdt(210)
Similarly, for case 2, Case 3, and Case 4, one obtains the followings, respectively.
B=−∂Ay∂zi+∂Ay∂xkEy=−1cdAydt(211)
B=∂Az∂yi−∂Az∂xjEz=−1cdAzdt(212)
B=0Ex=−∂ϕ∂x,Ey=−∂ϕ∂y,Ez=−∂ϕ∂z(213)
The numerical results for the four cases as functions of time are shown in Figs. 2–5.
9 Conclusion
First, we notice that the balance laws of micromorphic electromagnetic continuum consist of two parts: the thermomechanical (TM) part and the electromagnetic (EM) part. The TM part of the balance laws are the Conservation of Mass, Conservation of Microinertia, Balance of Linear Momentum, Balance of Moments of Momentum, Conservation of Energy, and Entropy Principle in Eringen’s micromorphic theory. The EM part of the balance laws is the Maxwell’s equations. The two parts are linked together by the electromagnetic body force, body moment, and energy source. We formulated a constitutive theory for thermo-visco-elastic-plastic-electromagnetic (TVEP-TM) materials. 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 well-known that the electric field E and the magnetic flux B are not independent of each other. To proceed with finite element analysis, one has to utilize the scalar potential ϕ and the vector potential A from which the electric field E and magnetic flux B can be derived.
Finally, five sets of differential equations are obtained for displacements, micromotions, temperature, scalar potential, and vector potential (cf. Eqs. (172)–(176)). The governing equations for displacements uand micromotionsχ are mainly acoustic wave equations; the governing equations for temperature T are mainly diffusion equations; The governing equations for scalar potential φ and vector potential A are mainly wave equations with a wave speed in the order of the speed of light. The huge difference between the acoustic wave speed and the speed of light dictates that solving these differential equations simultaneously is numerically impossible and physically meaningless. We solved this problem in two stages. To the electromagnetic waves in terms of the scalar and vector potentials, everything else seems to stand still. Therefore, in StageI, we solve the differential equations for temperature, displacements, and micromotions dynamically and seek the static or nearly static solutions. In StageII, based on the static or nearly static solutions obtained in StageI, we calculate the forcing terms for the differential equations of ϕ and A. Once ϕ or A is solved, the induced electric field E and magnetic flux B can be derived and obtained. For future study, it is worthwhile to recommend some papers related to micromorphic theory [20–24].
Acknowledgement: Not applicable.
Funding Statement: The authors received no specific funding for this study.
Author Contributions: The authors confirm contribution to the paper as follows: study conception and design, data collection, analysis and interpretation of results, draft manuscript preparation: Jiaoyan Li and James D. Lee. All authors reviewed and approved the final version of the manuscript.
Availability of Data and Materials: The authors confirm that the data supporting the findings of this study are available within the article.
Ethics Approval: Not applicable.
Conflicts of Interest: The authors declare no conflicts of interest.
Nomenclature
| A | A vector Ak, vector potential |
| Ak | In finite element analysis, Ak=NβAkβ |
| B | A vector Bk, magnetic flux, B=∇×A |
| C | Speed of light |
| D | A vector Dk, dielectric displacement vector |
| e | Internal energy density |
| E | A vector Ek, electric field, E=−∇ϕ−1c∂A∂t |
| E,Ee,Ep | Total strains, elastic strains, and plastic strains, E=Ee+Ep |
| f(Y) | A scalar-valued yield function of independent constitutive variables |
| f(Y)=0 | A hyper surface, named yield surface which separates the TVE-EM part and the TVEP-EM part of the material |
| f | A vector fi, body force (Thermomechanical) |
| F | A vector Fi, body force (Electromagnetic) |
| H | A vector Hk, magnetic field |
| h | Energy source (Thermomechanical) |
| ikl | A symmetric microinertia tensor, ikl=ilk |
| J | A vector Jk, total current vector |
| llm | A 2nd order tensor, body moment (Thermomechanical) |
| Llm | A 2nd order tensor, body moment (Electromagnetic) |
| M | Magnetization vector Mk, M≡B−H |
| mklm | A third order moment stress tensor |
| P | Polarization vector Pk, P≡D−E |
| q | Free charge density |
| qk | A vector, heat flux |
| RG | Fixed laboratory frame to where q,J,D,E,P,B,H,M are referred |
| RC | Co-moving reference frame to where q,J,D,E,P,B,H,M are referred |
| R | A generalized vector of internal variables, named as the hardening parameters i.e., R={G,g}. G has the same dimension as the strain,g is a scalar function. |
| skl | A symmetric microstress tensor, skl=slk, skl=sklm+sklem |
| sklm | The mechanical part of the microstress tensor, sklm=skl−12(PkEl∗+PlEk∗+Mk∗Bl+Ml∗Bk) |
| sklem | The electromagnetic part of the microstress tensor, sklem=−12(PkEl∗+PlEk∗+Mk∗Bl+Ml∗Bk) |
| tkl | A non-symmetric Cauchy stress tensor, tkl≠tlk, tkl=tklm+tklem |
| tklm | The mechanical part of the Cauchy stress tensor, tklm=tkl+PkEl∗+Mk∗Bl |
| tklem | The electromagnetic part of the Cauchy stress tensor, tklem=−PkEl∗−Mk∗Bl |
| T | Temperature |
| Tβ | In finite element analysis, T=NβTβ |
| U | A vector UK, displacement |
| UKα | In finite element analysis, UK=NβUKβ |
| v | Velocity vector vk |
| W | Energy source (Electromagnetic) |
| W | The set of internal variables, W={αp,βp,γp,R} |
| X | Eulerian coordinate of a particle P |
| X | Lagrangian coordinate of a particle P |
| Y | Set of independent objective constitutive variables |
| Z | Set if dependent objective constitutive variables |
| β | In finite element analysis, β means β−th shape function |
| α | First strain tensor, αKL≡xk,Kχ¯Lk−δKL |
| β | Second strain tensor (symmetric), βKL=χkKχkL−δKL=βLK |
| γ | Third strain tensor (third order), γKLM≡χ¯KkχkL,M |
| E≡|αβγ| | Total strains |
| Ee≡|αeβeγe| | Elastic strains |
| Ep≡|αpβpγp| | Plastic strains |
| αp,βp,γp | The corresponding plastic strains of α,β,γ |
| ρ | Mass density |
| ϕ | Scalar potential |
| ϕβ | In finite element analysis, ϕ=Nβϕβ |
| η | Entropy density |
| ψ | Helmholtz free energy density, ψ≡e−θη−ρ−1Ek∗Pk |
| θ | Absolute temperature |
| υlm | Microgyration tensor |
| σkl | Spin inertia tensor, σkl≡iml(υ˙km+υknυnm) |
| χkK | The affine micromotion which links ξk and ΞK, i.e., ξk=χkK(X,t)ΞK |
| χ¯Kk | The inverse of χkK, i.e., χkKχ¯Kl=δkl,χ¯KkχkL=δKL |
| χkKβ | In finite element analysis, χkK=NβχkKβ |
| ξα | The Eulerian coordinates of α−th material point in the particle |
| Ξα | Vector of α−th material point that attached to P, representing the inner structure of P |
Appendix A
It is noticed that the quantities D,E,P,B,H,M are all referred to a fixed laboratory frame RG. On the other hand, D∗,E∗,P∗,B∗,H∗,M∗ are referred to a co-moving frame RC. In relativistic electromagnetic theory, Jackson [10] indicated and explained that the Lorentz Transformation of E and B fields can be expressed as
E∗=γ(E+β×B)−γ2γ+1β(β⋅E),E=γ(E∗−β×B∗)−γ2γ+1β(β⋅E∗)(A1)
B∗=γ(B−β×E)−γ2γ+1β(β⋅B),B=γ(B∗+β×E∗)−γ2γ+1β(β⋅B∗)(A2)
where
β≡vc,β=β⋅β,γ≡11−β⋅β=11−β2(A3)
This means the inverse can be found by interchanging the star and the no-star quantities and by changing β to −β.
To prove Jackson’s formula, Eqs. (A1) and (A2), first notice that
γ2(1−β2)=1,−2γ3γ+1+γ2+γ4(γ+1)2β2=0(A4)
β×(β×E)=0β×(β×B)=0β×(β×E)=β(β⋅E)−β2Eβ×(β×B)=β(β⋅B)−β2B(A5)
Then substitute E into E* of Eq. (A1) and it results
E=γ(E∗−β×B∗)−γ2γ+1β(β⋅E∗)=γ{[γ(E+β×B)−γ2γ+1β(β⋅E)]−β×[γ(B−β×E)−γ2γ+1β(β⋅B)]}−γ2γ+1β{β⋅[γ(E+β×B)−γ2γ+1β(β⋅E)]}=γ{γE−γ2γ+1β(β⋅E)+γ[β(β⋅E)−β2E]}−γ2γ+1β{β⋅[γ(E−γ2γ+1β(β⋅E))]}+γ2β×B−γβ×[γB−γ2γ+1β(β⋅B)]=−γ2γ+1β{β⋅[−γ2γ+1β(β⋅E)]}=E(γ2−γ2β2)+β(β⋅E){−γ3γ+1+γ2−γ3γ+1+γ4(γ+1)2β2}=E(A6)
Now substitute B into B* of Eq. (A2) and it results
B=γ(B∗+β×E∗)−γ2γ+1β(β⋅B∗)=γ{[γ(B−β×E)−γ2γ+1β(β⋅B)]+β×[γ(E+β×B)−γ2γ+1β(β⋅E)]}−γ2γ+1β{β⋅[γ(B−β×E)−γ2γ+1β(β⋅B)]}=γ{−γβ×E+β×[γ(E)−γ2γ+1β(β⋅E)]}+γ2γ+1β{β⋅[γ(β×E)]}+γ{[γB−γ2γ+1β(β⋅B)]+β×[γ(β×B)]}−γ2γ+1β{β⋅[γ(B)−γ2γ+1β(β⋅B)]}=γ{[γB−γ2γ+1β(β⋅B)]+γ[β(β⋅B)−β2B]]}−γ2γ+1β{β⋅[γ(B)−γ2γ+1β(β⋅B)]}=B{γ2−γ2β2}+β(β⋅B){−γ3γ+1+γ2−γ3γ+1+γ4β2(γ+1)2}=B(A7)
This proves that Eqs. (A1) and (A2) are valid formula in Lorentz Transformation.
Following the same pattern of Eqs. (A1) and (A2), we write
D∗=γ(D+β×H)−γ2γ+1β(β⋅D)(A8)
H∗=γ(H−β×D)−γ2γ+1β(β⋅H)(A9)
M∗=γ(M+β×P)−γ2γ+1β(β⋅M)(A10)
P∗=γ(P−β×M)−γ2γ+1β(β⋅P)(A11)
Similar to the derivation of Eqs. (A6) and (A7), one may prove that
D=γ(D∗−β×H∗)−γ2γ+1β(β⋅D∗)(A12)
H=γ(H∗+β×D∗)−γ2γ+1β(β⋅H∗)(A13)
M=γ(M∗−β×P∗)−γ2γ+1β(β⋅M∗)(A14)
P=γ(P∗+β×M∗)−γ2γ+1β(β⋅P∗)(A15)
Since polarization and magnetization are defined as
P=D−E,P∗=D∗−E∗M=B−H,M∗=B∗−H∗(A16)
One may verify that
P=D−E=γ(D∗−β×H∗)−γ2γ+1β(β⋅D∗)−γ(E∗−β×B∗)+γ2γ+1β(β⋅E∗)=γ(P∗+β×M∗)−γ2γ+1β(β⋅P∗)=P(A17)
M=B−H=γ(B∗+β×E∗)−γ2γ+1β(β⋅B∗)−γ(H∗+β×D∗)+γ2γ+1β(β⋅H∗)=γ(M∗−β×P∗)−γ2γ+1β(β⋅M∗)=M(A18)
Recall that γ≡(1−β2)−12, the Taylor series expansion of γ about β=0 is obtained as
γ≈1+12β2+34β4+....(A19)
For the time being, let’s keep the Taylor series only to β2. Then we have
γ2γ+1≈1+β22+12β2≈12+38β2(A20)
Therefore, for non-relativistic electromagnetic theory we have
E∗=E+β×B,E=E∗−β×B∗(A21)
B∗=B−β×E,B=B∗+β×E∗(A22)
D∗=D+β×H,D=D∗−β×H∗(A23)
H∗=H−β×D,H=H∗+β×D∗(A24)
P∗=P−β×M,P=P∗+β×M∗(A25)
M∗=M+β×P,M=M∗−β×P∗(A26)
For illustrative purpose, let’s substitute E into E* and substitute B* into B of Eqs. (A21) and (A22) and it results
E=E′−β×B′=E+β×B−β×{B−β×E}=E+β×(β×E)=E+β(β⋅E)−β2E≈E(A27)
B=B′+β×E′=B−β×E+β×(E+β×B)=B+β×(β×B)=B+β(β⋅B)−β2B≈B(A28)
Notice that the approximations in Eqs. (A27) and (A28) are made due to the assumption that non-relativistic electromagnetic theory is valid up to small β≡vc≈0, β2 is not included because it is even smaller.
Appendix B
Recall the following equations
ρe˙=mklmυlm,k+tkl(vl,k−υlk)+sklυkl−qk,k+ρh+W(70)
ρη˙+∇⋅(q/θ)−ρh/θ≥0(71)
W=E∗⋅(P˙+P∇⋅v)−M∗⋅B˙+J∗⋅E∗(75)
ψ≡e−θη−ρ−1Ek∗Pk(79)
One may verify that Eq. (33) can be re-written as
ρoe˙=MKLMγ˙KLM+TKLα˙KL+SKLβ˙KL−QK,K+ρoh+JW(A29)
where the description of the Lagrangian stresses and heat flux, MKLM,TKL,SKL,QK, and the description of the Eulerian stresses and heat flux, mklm,tkl,skl,qk. are related as [5]
TKL=JtklXK,kχlL,tkl=J−1TKLxk,Kχ¯LlSKL=12Jsklχ¯Kkχ¯Ll,skl=2J−1SKLχkKχlLMKLM=JmmklXM,mχkKχ¯Ll,mmkl=J−1MKLMxm,Mχ¯KkχlLQK=JqkXK,k,qk=J−1QKxk,K(A30)
and the corresponding strain rates are defined as
α˙KL≡aklxk,Kχ¯Ll≡(vl,k−υlk)xk,Kχ¯Llβ˙KL≡2bklχkKχlL≡(υkl+υlk)χkKχlLγ˙KLM≡cklmχ¯KkχlLxm,M≡υkl,mχ¯KkχlLxm,M(A31)
One may readily verify that Eqs. (34) and (39) can be re-written as
ρoθη˙−1θQKθ,K+QK,K−ρoh≥0(A32)
ρoψ≡ρo(e−θη)−JEk∗Pk(A33)
Define EK∗,BK,PK,MK∗, and JK∗ as
EK∗≡Ek∗xk,KBK≡Bkxk,KPK≡JPkXK,kMK∗≡JMk∗XK,kJK∗≡JJk∗XK,k(A34)
One may readily verify that
Ek∗≡EK∗XK,kBk≡BKXK,kPk≡J−1PKxk,KMk∗≡J−1MK∗xk,KJk∗≡J−1JK∗xk,K(A35)
E˙K∗=(E˙k∗+El∗vl,k)xk,KB˙K=(B˙k+Blvl,k)xk,KEK∗PK=JEk∗Pk(A36)
Theorem 1: The entropy principle can be expressed as
−ρo(ψ˙+ηθ˙)−1θQKθ,K−EK∗P˙K−E˙K∗PK+MKLMγ˙KLM+TKLα˙KL+SKLβ˙KL+JW≥0(A37)
Proof:
Rewrite Eq. (A33) as
ρoψ˙=ρoe˙−ρoθ˙η−ρoθη˙−E˙K∗PK−EK∗P˙K(A38)
Following Eq. (A33), Eq. (81) can be re-written as
ρoθη˙−1θQKθ,K+QK,K−ρoh=−ρo(ψ˙+ηθ˙)+ρoe˙−EK∗P˙K−E˙K∗PK−1θQKθ,K+QK,K−ρoh(A39)
Substituting Eq. (A29) into Eq. (A39) results
−ρo(ψ˙+ηθ˙)−EK∗P˙K−E˙K∗PK+MKLMγ˙KLM+TKLα˙KL+SKLβ˙KL−1θQKθ,K+JW≥0(A40)
Thus, Theorem 1 is proved. □
Theorem 2: The entropy principle can also be expressed as
−ρo(ψ˙+ηθ˙)+MKLMγ˙KLM+TKLα˙KL+SKLβ˙KL−1θQKθ,K−PKE˙K∗−MK∗B˙K+JK∗EK∗+Jvl,k(PkEl∗+Mk∗Bl)≥0(A41)
Proof:
Now we rewrite Eq. (38) as
JW=J{Ek∗(P˙k+Pk∇⋅v)−Mk∗B˙k+Jk∗Ek∗}≡α1+α2+α3(A42)
and derive the followings
α1≡J{Ek∗Pk∇⋅v+Ek∗ddt[J−1PKxk,K]}=J{Ek∗Pk∇⋅v+Ek∗[−J−1∇⋅vPKxk,K+J−1P˙Kxk,K+J−1PKvk,lxl,K]}=J{Ek∗[J−1P˙Kxk,K+J−1PKvk,lxl,K]}=EK∗P˙K+JEl∗Pkvl,k(A43)
α2≡−JMk∗B˙k=−J(J−1MK∗xk,K)ddt[BLXL,k]=−MK∗xk,K[B˙LXL,k−BLvl,kXL,l]=−MK∗B˙K+Jvl,kMk∗Bl(A44)
α3≡JJk∗Ek∗=J(EK∗XK,k)(J−1JL∗xk,L)=JK∗EK∗(A45)
Therefore, we arrive at
JW=EK∗P˙K−MK∗B˙K+JK∗EK∗+Jvl,k(PkEl∗+Mk∗Bl)(A46)
Since SKL=SLK and skl=slk, one may readily verify that
TKLα˙KL=Jtkl(vl,k−υlk),SKLβ˙KL=12Jskl(υkl+υlk)(A47)
Now one may further derive TKLα˙KL+SKLβ˙KL+Jvl,k(PkEl∗+Mk∗Bl) as
TKLα˙KL+SKLβ˙KL+Jvl,k(PkEl∗+Mk∗Bl)=J{tkl(vl,k−υlk)+sklυlk+vl,k(PkEl∗+Mk∗Bl)}=J{vl,k(tkl+PkEl∗+Mk∗Bl)+(skl−tkl)υlk}(A48)
Define the mechanical part of the Cauchy stresses and microstresses, tklm and sklm, and the electromagnetic part of the Cauchy stresses and microstresses, tklem and sklem, as
tkl=tklm+tklem,skl=sklm+sklem(A49)
tklem≡−PkEl∗−Mk∗Bl,tklm≡tkl+PkEl∗+Mk∗Blsklem≡−PkEl∗−Mk∗Bl,sklm≡skl+PkEl∗+Mk∗Bl(A50)
Now one may verify that Eq. (A48) as
TKLα˙KL+SKLβ˙KL+Jvl,k(PkEl∗+Mk∗Bl)=J{tklm(vl,k−υlk)+sklmυlk}(A51)
Define
SKLm≡12Jsklmχ¯Kkχ¯Ll(A52)
Then verify that
sklm=2J−1SMNmχkMχlN(A53)
Notice that SKLm≠SLKm,sklm≠slkm, and β˙KL=β˙LK. Now we further derive
SKLmβ˙KL=S¯KLmβ˙KL=12Jχ¯Kkχ¯Lls¯klmχmKχnL(υmn+υnm)=12J(υkl+υlk)s¯klm=Jυlks¯klm(A54)
Notice that S¯KLm is the average of SKLm and SLKm, i.e., S¯KLm≡12(SKLm+SLKm)=S¯LKm. Now we follow the pattern of Eqs. (A52) and (A53) and obtain the following
S¯KLmβ˙KL=Jυlks¯klm(A55)
where s¯klm is the average of sklm and slkm, i.e., s¯klm≡12(s¯klm+s¯lkm)=s¯lkm. From Eq. (A50), we have
s¯klm=skl+12(PkEl∗+PlEk∗+Mk∗Bl+Ml∗Bk)(A56)
Now Eq. (A51) can be re-written as
TKLα˙KL+SKLβ˙KL+Jvl,k(PkEl∗+Mk∗Bl)=J{tklm(vl,k−υlk)+sklmυlk}=TKLmα˙KL+S¯KLmβ˙KL(A57)
Now the conservation of energy and the Clausius-Duhem inequality can be expressed as
ρoe˙=MKLMγ˙KLM+TKLmα˙KL+S¯KLmβ˙KL−QK,K−PKE˙K∗−MK∗B˙K+JK∗EK∗+ρoh(A58)
−ρo(ψ˙+ηθ˙)+MKLMγ˙KLM+TKLmα˙KL+S¯KLmβ˙KL−1θQKθ,K−PKE˙K∗−MK∗B˙K+JK∗EK∗≥0(A59)
From now on if there is no ambiguity, we drop the bar on top of SKLm and sklm. □
Appendix C
First, we link the displacement, micromotion, 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αχkK=NαχkKαT=NαTαϕ=NαϕαAk=NαAkα(A60)
The weak form of Eq. (165) can be written as
∫V{(TKLχ¯Ll),K+ρo(fl−v˙l)+JFl}δuldV=0(A61)
Term by term, we may obtain the followings:
∫V(TKLχ¯Ll),KδuldV=∮STKLχ¯LlnKδuldS−∫VTKLχ¯Llδul,KdV=δUMα{∮STKLχ¯LlnKδlMNαdS−∫VTKLχ¯LlδlMNα,KdV}=δUMα{∫SσF^MNαdS−∫VTKLχ¯LlδlMNα,KdV}(A62)
where on the boundary the surface force is specified as
F^M=TKLχ¯LlnKδlM(A63)
Also, we define
fMα1≡∫SFF^MNαdSfMα2≡−∫VTKLχ¯LlδlMNα,KdV(A64)
The other two terms are derived to be
∫V(ρofl+JFl)δuldV=δUMα∫V(ρofM+JFM)NαdV≡δUMαfMα3(A65)
∫Vρov˙lδuldV=U¨MβδUMα∫VρoNαNβdV≡U¨MβδUMαMαβ(A66)
Now Eq. (A61) can be re-written as
δUMα{fMα1+fMα2+fMα3−MαβU¨Mβ}=0(A67)
Because Eq. (A67) has to be valid for arbitrary virtual displacement δUMα, the finite element equation for the balance of linear momentum becomes
MαβU¨Mβ=fMα1+fMα2+fMα3(A68)
For the finite element formulation of the Balance of Moment of Momentum, first we recall the microgyration and spin inertia as [7]
υkl=χ˙kKχ¯Kl(A69)
σkl≡iml(υ˙km+υknυnm)≡imlCkm(A70)
From Eq. (A69), one may obtain
χ˙kK=υklχlK(A71)
Now we derive the followings
χ¨kK=dχ˙kKdt=ddt(υklχlK)=υ˙klχlK+υklχ˙lK=υ˙klχlK+υklυlnχnK=(υ˙kn+υklυln)χnK=CknχnK(A72)
Define
δAkl≡χ¯KlδχkK(A73)
Then it is seen that
σklδAkl=imlCkmδAkl=IMLχmMχlLCkmδAkl=IMLχmMχlLCkmχ¯KlδχkK=IMLχmMCkmδχkL=IMLχ¨kMδχkL(A74)
Now the weak form of Eq. (166) can be written as
∫V{(MLMKχ¯LlχmM),K+TMLxm,Mχ¯Ll−2SLMχlLχmM+ρo(llm−σlm)+JLlm}δAlmdV=0(A75)
The inertia term can be derived as
∫VρoσlmδAlmdV=∫VρoIKLχ¨kKδχkLdV=χ¨kLβδχkKα∫VρoIKLNαNβdV≡δχkKαMKLαβχ¨kLβ(A76)
If we assume the material is spin isotropic, i.e., IKL=IδKL, then Eq. (A76) is reduced to
∫VρoσlmδAlmdV=δχkKαIMαβχ¨kKβ(A77)
The first term in the weak form becomes [5]
∫V(MLMKχ¯LlχmM),KδAlmdV=∫V(MLMKχ¯LkχmM),KδAkmdV=∫V(MLMKχ¯LkχmM),KδAkmdV=∫V(MLMPχ¯LkχmM),Pχ¯KmδχkKdV=δχkKα{∫V(MLMPχ¯LkχmM),Pχ¯KmNαdV−∫VMLMPχ¯LkχmM(χ¯KmNα),PdV}=δχkKα{∫SMM^kKNαdS−∫VMLMPχ¯LkχmMχ¯Km,PNαdV−∫VMLKMχ¯LkNα,MdV}(A78)
where on the boundary the surface moment is specified as
M^kK≡MLKMχ¯LknM(A79)
The second, third, and fourth terms become
∫VTMLxm,Mχ¯LlδAlmdV=∫VTMLxm,Mχ¯LkδAkmdV=∫VTMLxm,Mχ¯Lkχ¯KmδχkKdV=δχkKα∫VTMLxm,Mχ¯LkNαdV(A80)
∫V−2SLMχlLχmMδAlmdV=−∫V2SLMχkLχmMδAkmdV=−δχkKα∫V2SLMχkLχmMχ¯KmNαdV=−δχkKα∫V2SLKχkLNαdV(A81)
Define
ϕkKα1≡∫SMM^kKNαdSϕkKα2≡−∫VMLMPχ¯LkχmMχ¯Km,PNαdVϕkKα3≡−∫VMLKMχ¯LkNα,MdVϕkKα4≡∫VTMLxm,Mχ¯Lkχ¯KmNαdVϕkKα5≡∫V−2SLKχkLNαdVϕkKα6≡∫V{ρolkm+JLkm}χ¯KmNαdV(A82)
With the assumption of spin isotropy, i.e., IKL=IδKL, Eq. (A75) can now be re-written as
δχkKα{ϕkKα1+ϕkKα2+ϕkKα3+ϕkKα4+ϕkKα5+ϕkKα6−IMαβχ¨kKβ}=0(A83)
Because Eq. (A83) has to be valid for arbitrary virtual micromotion, δχkKα the finite element equation for the balance law of moment of momentum becomes
IMαβχ¨kKβ=ϕkKα1+ϕkKα2+ϕkKα3+ϕkKα4+ϕkKα5+ϕkKα6(A84)
Based on the energy equation, the weak form of the conservation law of energy can be written as (cf. Eq. (168)) [5]
∫V{ρoγ(T+To)ToT˙+(T+To)[A¯KLα˙KLe+B¯KLβ˙KLe]−TKLmeα˙KLp−SKLmeβ˙KLp−MKLMeγ˙KLMp−TKLmdα˙KL−SKLmdβ˙KL−MKLMdγ˙KLM+PKE˙K∗+MK∗B˙K−JK∗EK∗+QK,K−ρoh}δTdV=0(A85)
The first term in the weak form becomes
∫Vρoγ(T+To)ToT˙δTdV=T˙βδTα∫Vρoγ(NξTξ+To)ToNαNβdV≡T˙βδTαGαβ(A86)
The second term yields
∫V−(T+To)[A¯KLα˙KLe+B¯KLβ˙KLe]δTdV=δTα∫V−(T+To)[A¯KLα˙KLe+B¯KLβ˙KLe]NαdV≡δTαqα1(A87)
The third term results
∫V{TKLmeα˙KLp+SKLmeβ˙KLp+MKLMeγ˙KLMp}δTdV=δTα∫V{TKLmeα˙KLp+SKLmeβ˙KLp+MKLMeγ˙KLMp}NαdV≡δTαqα2(A88)
The fourth term is identified and derived as
∫V{TKLdα˙KL+SKLdβ˙KL+MKLMdγ˙KLM}δTdV=δTα∫V{TKLdα˙KL+SKLdβ˙KL+MKLMdγ˙KLM}NαdV≡δTαqα3(A89)
The fifth term is identified and derived as
∫V{−PKE˙K∗−MK∗B˙K+JK∗EK∗}δTdV=δTα∫V{−PKE˙K∗−MK∗B˙K+JK∗EK∗}NαdV≡δTαqα4(A90)
The sixth term can be derived as
∫V−QK,KδTdV=−∮SQKnKδTdS+∫VQKδT,KdV=δTα{−∫SqQ^NαdS+∫VQKNα,KdV}≡δTα{qα5+qα6}(A91)
The last term becomes
∫VρohδTdV=δTα∫VρohNαdV≡δTαqα7(A92)
Now Eq. (A85) can be re-written as
δTα{GαβT˙β−[qα1+qα2+qα4+qα4+qα5+qα6+qα7]}=0(A93)
Therefore, the finite element equation for conservation law of energy becomes
GαβT˙β=qα1+qα2+qα3+qα4+qα5+qα6+qα7(A94)
where
qα1≡∫V−(T+To)[A¯KLα˙KLe+B¯KLβ˙KLe]NαdV(A95)
qα2≡∫V{TKLmeα˙KLp+SKLmeβ˙KLp+MKLMeγ˙KLMp}NαdV(A96)
qα3≡∫V{TKLdα˙KL+SKLdβ˙KL+MKLMdγ˙KLM}NαdV(A97)
qα4≡∫V{−PKE˙K∗−MK∗B˙K+JK∗EK∗}NαdV(A98)
qα5≡−∫SqQ^NαdS(A99)
qα6≡∫VQKNα,KdV(A100)
qα7≡∫VρohNαdV(A101)
For scalar potential, the weak form of Eq. (159) is expressed as
∫v{1c¯2ϕ¨−(ϕ,i),i−a}δϕdv=0(A102)
Term by term one obtains
∫v1c¯2ϕ¨δϕdv=δϕαϕ¨β∫v1c¯2NβNβdv=δϕαmαβϕ¨β(A103)
∫v−(ϕ,i),iδϕdv=∫v{−(ϕ,iδϕ),i+ϕ,iδϕ,i}dv=δϕα∮s−ϕ,iniNαds+δϕαϕβ∫vNα,iNβ,idv=δϕα∫sξ−ϕ,iniNαds+δϕαϕβ∫vNα,iNβ,idv=−δϕαξ^α+δϕαkαβϕβ(A104)
∫v−aδϕdv=δϕα∫v−aNαdv=−δϕαa~α(A105)
Because the weak form has to be valid for any arbitrary δϕα, one obtains
mαβϕ¨β+kαβϕβ=ξ^α+a~α(A106)
where sξ is part of the enclosing surface of the volume v on which ϕ,ini is specified as ξ; and
mαβ≡∫v1c¯2NαNβdvkαβ≡∫vNα,iNβ,idvξ^α≡∫sξξNαds≡∫sξϕ,iniNαdsa^α=∫vaNαdv(A107)
For vector potential, the weak form of Eq. (159) is expressed as
∫v{1c¯2∂2Ak∂t2−Ak,ii−Ck}δAkdv=0(A108)
Term by term one obtains
∫v1c¯2∂2Ak∂t2δAkdv=δAkαA¨kβ∫v1c¯2NαNβdv=δAkαmαβA¨kβ(A109)
∫v−Ak,iiδAkdv=∫v{−[Ak,iδAk],i+Ak,iδAk,i}dv=δAkα∮s−Ak,iniNαds+δAkαAkβ∫vNα,iNβ,idv=δAkα∫sς−A~kNαds+δAkαAkβ∫vNα,iNβ,idv≡δAkα{−A^kα+kαβAkβ}(A110)
∫v−CkδAkdv=−δAkα∫vCkNαdv≡−δAkαC^kα(A111)
Because the weak form has to be valid for any arbitrary δAkα, one obtains
mαβA¨kβ+kαβAkβ=A^kα+C^kα(A112)
where sς is part of the enclosing surface of the volume v on which Ak,ini is specified as A~k; and
A^kα≡∫sςAk,iniNαds≡∫sςA~kNαdsC^kα≡∫vCkNαdv(A113)
References
1. Lee JD, Eringen AC. Alignment of nematic liquid crystals. J Chem Phys. 1971;55(9):4504–8. doi:10.1063/1.1676781. [Google Scholar] [CrossRef]
2. Lee JD, Eringen AC. Boundary effects of orientation of nematic liquid crystals. J Chem Phys. 1971;55(9):4509–12. doi:10.1063/1.1676782. [Google Scholar] [CrossRef]
3. Lee JD, Eringen AC. Wave propagation in nematic liquid crystals. J Chem Phys. 1971;54(12):5027–34. doi:10.1063/1.1674793. [Google Scholar] [CrossRef]
4. Lee JD, Eringen AC. Continuum theory of smectic liquid crystals. J Chem Phys. 1973;58(10):4203–11. doi:10.1063/1.1678976. [Google Scholar] [CrossRef]
5. Lee JD, Li J. Advanced continuum theories and finite element analyses. Singapore: World Scientific; 2019. doi:10.1142/11312. [Google Scholar] [CrossRef]
6. Eringen AC. Mechanics of continua. Huntington, NY, USA: Robert E. Krieger Publishing Co.; 1980. [Google Scholar]
7. Eringen AC. Microcontinuum field theories. New York, NY, USA: Springer; 1999. doi:10.1007/978-1-4612-0555-5. [Google Scholar] [CrossRef]
8. Eringen AC. Microcontinuum field theories II: fluent media. New York, NY, USA: Springer; 2001. [Google Scholar]
9. Eringen AC. Simple microfluids. Int J Eng Sci. 1964;2(2):205–17. doi:10.1016/0020-7225(64)90005-9. [Google Scholar] [CrossRef]
10. Jackson JD. Classical electrodynamics. Hoboken, NJ, USA: John Wiley & Sons, Inc.; 1999. [Google Scholar]
11. de Groot SR, Suttorp IG. Foundation of electrodynamics. Amsterdam, The Netherlands: North-Holland Publishing Company; 1972. [Google Scholar]
12. Eringen AC, Maugin GA. Electrodynamics of continua I: foundations and solid media. New York, NY, USA: Springer; 1990. doi:10.1007/978-1-4612-3226-1. [Google Scholar] [CrossRef]
13. Eringen AC, Maugin GA. Electrodynamics of continua II: fluids and complex media. New York, NY, USA: Springer; 1990. doi:10.1007/978-1-4612-3236-0. [Google Scholar] [CrossRef]
14. Coleman BD, Noll W. The thermodynamics of elastic materials with heat conduction and viscosity. Arch Ration Mech Anal. 1963;13(1):167–78. doi:10.1007/BF01262690. [Google Scholar] [CrossRef]
15. Casey J. On elastic-thermo-plastic materials at finite deformations. Int J Plast. 1998;14(1–3):173–91. doi:10.1016/S0749-6419(97)00047-8. [Google Scholar] [CrossRef]
16. Green AE, Naghdi PM. A general theory of an elastic-plastic continuum. Arch Ration Mech Anal. 1965;18(4):251–81. doi:10.1007/BF00251666. [Google Scholar] [CrossRef]
17. Onsager L. Reciprocal relations in irreversible processes. I. Phys Rev. 1931;37(4):405–26. doi:10.1103/physrev.37.405. [Google Scholar] [CrossRef]
18. Onsager L. Reciprocal relations in irreversible processes. II. Phys Rev. 1931;38(12):2265–79. doi:10.1103/physrev.38.2265. [Google Scholar] [CrossRef]
19. Li J, Thiesen J, Robert KP, Yang X, Lee JD. Finite element analysis of fracture mechanics in micromorphic theory. J Micromech Mol Phys. 2022;7(2):143–56. doi:10.1142/s2424913021420170. [Google Scholar] [CrossRef]
20. Chen Z, Chu X, Duan Q. A micromorphic peridynamic model and the fracture simulations of quasi-brittle material. Eng Fract Mech. 2022;271:108631. doi:10.1016/j.engfracmech.2022.108631. [Google Scholar] [CrossRef]
21. Reges PDN, Pitangueira RLS, Silva LL. Elastic degradation models for the micromorphic continuum. Int J Non Linear Mech. 2023;154(5):104450. doi:10.1016/j.ijnonlinmec.2023.104450. [Google Scholar] [CrossRef]
22. Hütter G. Interpretation of micromorphic constitutive relations for porous materials at the microscale via harmonic decomposition. J Mech Phys Solids. 2023;171(1):105135. doi:10.1016/j.jmps.2022.105135. [Google Scholar] [CrossRef]
23. Malik A, Hütter G, Abendroth M, Kiefer B. Micromorphic FE2 simulation of plastic deformations of foam structures. Int J Mech Sci. 2024;282(3):109551. doi:10.1016/j.ijmecsci.2024.109551. [Google Scholar] [CrossRef]
24. Enakoutsa K, Burson R, Li Y. Hybrid modeling of porous metal plasticity: integrating micromorphic theories with data-driven methods. Next Mater. 2025;9(1):101145. doi:10.1016/j.nxmate.2025.101145. [Google Scholar] [CrossRef]
Login
Register
Submit a Paper
Propose a Special lssue
Open Access
View Full Text
Download PDF
Copyright © 2026 The Author(s). Published by Tech Science Press.
Downloads
Citation Tools
