Numerical experiment based on non-linear micropolar finite element to study micro-rotations generated by the non-symmetric Maxwell stress tensor

The main aim of the present work is to investigate the role of the Maxwell stress tensor in the study of active materials. Despite the importance of this tensor in modeling mechatronic devices used in sophisticated applications, its non–symmetry still generates controversies in the literature, probably because classical continuum mechanics assumes a symmetric Cauchy stress, although the sum of Cauchy and Maxwell stresses is non–symmetric. In the framework of generalised continuum mechanics–a more advanced formalism than the classical one–, each material point has an associated microstructure so that the micro–rotations of the electric/magnetic dipoles present in real active materials may be simulated. To this end, a modified total stress formulation, including an angular momentum balance, is developed and implemented into a finite element research code using a complex–step formulation. It is concluded that generalised mechanics allows for incorporating both symmetric and non–symmetric contributions of the Maxwell tensor. Consequently, the rotations generated by the electromagnetic field may be analysed. The influence of the complete Maxwell tensor in a magnetostrictive actuator is studied by several magneto–mechanical numerical experiments of a Terfenol–D rod encircled by air, and several conclusions are highlighted.


Introduction
Active materials, such as piezoelectrics and piezomagnetics, are used in mechatronic devices due to their inherent ability to combine mechanical and electromagnetic fields by two separate couplings: the first can be modelled by constitutive equations-obtained from thermodynamic potentials-, and the second by electromagnetic body forces, commonly incorporated into the linear momentum balance by the Maxwell Stress Tensor (MST).
Physically, the origin of these couplings comes from two separate interactions, see [27,31]: -Short-range forces (local effects) due to interatomic interactions among material points inside the body that generate both forces and couples.-Long-range forces (non-local effects) that emerge from the interactions between material points and external fields.
The short-range forces, see [37,39], are closely related with the piezoelectric/piezomagnetic constitutive equations, some reported in [15].As mentioned, the contribution of these forces to the linear momentum balance is often obtained from thermodynamic potentials, see [10,13].
The long-range forces are commonly expressed as the divergence of the MST for an amenable implementation in numerical approaches such as the Finite Element (FE) method.However, the MST is non-symmetric in transversely isotropic materials such as piezoelectrics and piezomagnetics [25,34,37].This lack of symmetry has generated a compre-hensive discussion in the literature since Noether's theorem [21] states that the angular momentum must be conserved, and non-symmetric stress tensors do not automatically fulfil this condition.For this reason, [31] reported that the nonsymmetric MST does not have physical meaning since it is caused by the transformation of the electromagnetic longrange body forces into the short-range divergence of the MST.
Despite the importance of the MST in modelling mechatronic devices (for instance, for vibrations and noise control), it should be considered in many research works, for example, in [29,32].This is due to two main difficulties: the non-linearity the possible non-symmetry The first difficulty is nowadays easily overcome by nonlinear algorithms and sophisticated computer codes.But the second remains unsolved since the proper solution is not clear from a theoretical point of view, even if it is essential since the MST skew-symmetric part can generate electromagnetic torques as argued in [6,14].
On this ground, two main approaches based on Classical Continuum Mechanics (CCM) may be found in the literature to model active materials considering the MST: (i) Analytically, [1,2,19,37] introduce a symmetric total stress tensor as the sum of two non-symmetric tensors: a Cauchy-like stress tensor and the MST itself.FE models using this total stress are reported in [11,22,28,30,38].(ii) The numerical approach from [25] allows us to calculate electromagnetic torques and forces with a nodal force method.The torques are expanded into non-local (longrange) and local (short-range) terms, the latter due to the MST non-symmetric part.
Although the controversy on the MST non-symmetry was partially clarified through physical arguments in [19,37], the main novelty of the present work is to establish the CCM limitations (including the use of the MST itself) to model active materials.For this purpose, a Generalised Continuum Mechanics (GCM) approach based on the Micropolar Mechanics (MM) theory is developed; see [17] for a general description of MM.Roughly, this theory extends CCM assuming that each material point in the continuum is associated with a micro-rigid body and, consequently, three micro-rotations are incorporated to the FE degrees of freedom (dofs).
Based on MM, in the current work both linear and angular momentum balances are stated, concluding that the MST holds a skew-symmetric part.Furthermore, the electromagnetic torque generated by this non-symmetry may only be studied with MM since, in the framework of CCM, each , mechanical traction t i , couple stress vector m i , and other body forces f i material point lays in the position of its centre of mass [24].Then the point only transmits linear momentum to neighbours, not micro-rotations.As demonstrated in [19,37], the absence of rotation implies that the MST must be symmetric for basic CCM.
Finally, and for completeness, a numerical study based on the FE method is developed to obtain some conclusions on the influence of the MST using MM.In particular, a mini-magnetostrictive actuator made out of the transversely isotropic Terfenol-D material is simulated.Four numerical experiments are developed, demonstrating that the microrotations can be of the macro order for magnetic fields of normal levels.Also, regardless of intensity, not all magnetic fields can produce rotations.
Dynamic effects are not considered in the present work for simplicity, and for the FE formulation, only the magnetomechanical coupling is considered.
For clarity, index notation is used throughout the present work.Non-symmetric tensors are denoted by (•) i j , with their symmetric and skew-symmetric parts denoted by (•) (i j) and (•) [i j] , respectively.In addition, the alternating third-order and identity second-order tensors are denoted by e i jk and δ i j .Finally, the Einstein summation convention is used throughout the document.∈ Ω E M , mechanical forces represented by traction t i and couple stress m i vectors so that {t i , m i } ∈ Ω M .Other body forces f i ∈ Ω O are included for completeness due to long-range effects such as gravitation.

Balance equations
As introduced before, MM considers that each material point has an associated microstructure mathematically represented by a trihedron at each material point C. Consequently, the microstructure is modelled as a rigid body adding three micro-rotations θ i to the three classical macrodisplacements u i .
For the most general case, the electromagnetic forces f E M i are composed of: -Long-range forces from electromagnetic source: free electric charge density ρ f q and free electric current j f i .-Short-range forces due to polarisation P i and magnetisation M i vectors related to the material media.
The vectors t i and m i exclusively represent short-range forces.The former is the classical traction related to a generic stress tensor σ i j by the Cauchy relation of ( 1) top.The latter is due to the microstructure, and it is related to the couple stress tensor τ i j of (1) bottom, also by a Cauchy-like relation: ( Notice that both σ i j and τ i j are a priori non-symmetric tensors.

Electromagnetic forces
Continuum Electrodynamics is founded on a set of four empirical laws called Maxwell equations [28] expressed as: where E i , D i , H i , and B i denote electric field, electric displacement, magnetic field, and magnetic induction, respectively.In addition, t denotes time and ∂ partial derivative.
The electromagnetic constitutive equations that relate P i and M i vectors are defined as: where 0 and μ 0 are the vacuum permittivity and permeability, respectively.The balance of linear electromagnetic momentum is obtained from Poynting's theorem [28] and in quasi-static regime is expressed as: where σ E M ji is the mathematical symbol for the MST.There exist at least four expressions [9] for f E M i or σ E M ji , but the present work uses: According to [39] the previous expression satisfies the electromagnetic boundary conditions and may be obtained from a thermodynamic formulation as in [10] or from energymomentum considerations as in [8].
Since any tensor may uniquely be decomposed into symmetric and skew-symmetric parts, the MST is split into both by introducing (3) into (5) to give: where the second contribution represents the tensor form of the classical electromagnetic torque.
According to [23], applying the divergence theorem to the total MST allows to transform the long-range forces f E M i into short-range ones by: where t E M i denotes the Maxwell traction vector.

Total linear momentum balance
From Noether's theorem, the linear momentum must be conserved since the Lagrangian (function of movement) must be independent of the origin.Considering both mechanical and electromagnetic forces, the balance of linear momentum establishes the relationship between long-range and shortrange forces by: Applying the divergence theorem to the previous first integral, taking into account (1) top and (7), the local form of (8) becomes:

Total angular momentum balance
Also, from Noether's theorem, the angular momentum must be explicitly conserved since the Lagrangian function must be independent of the measurement angle.Depending on the material constitution, the MST skew-symmetric part may be: = 0, then the angular momentum balance is automatic stated.σ E M [i j] = 0, then the angular momentum balance must be explicitly guaranteed in the formulation.With x i denoting the position vector that locates each material point, for the second case and applying the definition of angular momentum, the second balance reads: Again, the divergence theorem may be applied to the first term on the left-hand side of (10).Taking into account (1) bottom and ( 9), the balance of angular momentum in local form reads:

On the symmetry of the Maxwell tensor
Based on the equations of the previous section, the main objectives are now: i) to analyse the CCM limitations to study the MST skew-symmetric part, ii) to formulate a new approach based on MM to include the MST skew-symmetric part.To address i), a revision of the different approaches for the MST is conducted.

Classical continuum mechanics approaches
In the CCM framework, the couple stress tensor and the micro-rotations are absent: τ ji = θ i = 0.In this context, at least three approaches exist to study the MST in the literature.The most important approaches are listed to highlight the absence of rotations.

Rinaldi
This approach was developed in [31] and is grounded on the lack of physical meaning in the mixture of long-and shortrange forces.Specifically, even though the long-range force f E M i (called maxwellian in the reference) may be mathematically expressed as a short-range force t E M i related with σ E M ji , as in (7), the reference argues that ( 7) is not unique but "arbitrary to within an additive divergence-less tensor".In this context, [20] also discusses the wrong "conceptual replacement" since t E M i has no physical sense "unless it is integrated over a closed surface".
Rinaldi concludes that the short-range forces must exclusively come from the classical and symmetric Cauchy stress tensor σ C (i j) .Consequently, the linear (9) and angular (11) balances simplify to: Finally, Rinaldi remarks that f E M i can only be replaced by σ E M ji, j if the electromagnetic torque is added to (12) bottom.

Eringen
This approach reported in [4] modifies the balance of linear electromagnetic momentum (4) by splitting f E M i as the sum of long-range Lorentz f L i and of short-range ponderomotive f P M i forces: Also, Eringen argues that f P M i only generates angular momentum, and therefore, it is not included in the linear momentum (9):

Total stress
This approach is the most published: [1,2,19,37].It is grounded in the definition of a symmetric total stress σ T as the sum of a non-symmetric and Cauchy-like tensor σ i j plus the MST tensor: As argued in [1], this σ T (i j) includes both long-range (electromagnetic) and short-range (mechanical) forces and satisfies boundary conditions.Along that line, the reference criticised Rinaldi's approach for not considering proper boundary conditions.Still within CCM, the linear (9) and angular (11) balances become: Fig. 2 Polarisable body composed of electric dipoles p e i defined by two electric charges q + , q − represented by rigid bodies (zoom): stretching and rotation of p e i under electric field This approach results in the same set of equations as (13) of Eringen from the previous subsection.

Micropolar mechanics
As commented, CCM represents each material point only by its centre of mass.A more sophisticated formalism such as MM can better describe the electric and magnetic dipole rotations at the micro-continuum scale.
Consider the polarisable medium shown in Fig. 2 (the same arguments hold for magnetisable media) composed of electric dipoles q i , defined by two electric charges q + , q − .As observed in the zoom and argued in [7], q i will stretch and rotate an angle θ i upon application of an electric field.CCM captures the stretch, but a GCM approach such as MM must be used to study the rotation.Under this approach, the balances of linear and angular momentum read: where σ C [ ji] may be considered null in a first approximation, see Sect. 4. Notice that with this nullity and according to [18], [ jk] is closely related to a parameter called coupling number, which depends on the microstructure (intrinsic lattice size scale) of each material.
In conclusion, the current MM approach has several advantages, allowing to: -sum the symmetric Cauchy σ C (i j) to the non-symmetric Maxwell σ E M ji , since the angular momentum balance is fulfilled by the additional (17) bottom.
-use the classical constitutive equations of piezoelectrics and/or piezomagnetics (see Sect. 4) since they are expressed using only the symmetric part σ C (i j) .-calculate the rotation of the electric/magnetic dipoles, and from it, the generated electromagnetic torque τ ji , relevant for mechatronic devices.
-recover the standard total stress method if the microstructure is not considered and τ ji = 0.

Magnetostrictive governing equations
This section particularises the MM governing equations for the magnetostrictive problem used in the numerical experiments of Sect.6.Seven dofs are considered: three macro-displacements u i , three micro-rotations θ i , and the scalar magnetic potential ϕ.
The governing equations comprise balance (linear and angular) momenta, compatibility equations, boundary conditions, and material constitution.For the sake of clarity, this section summarises all equations.

Balance equations
The balance equations for the MM approach given in (17) are rewritten here without the separation into symmetric and skew-symmetric contributions.In addition, the magnetic field is balanced by the Gauss law, the second of (2) to obtain:

Assuming σ C
[ ji] = 0, the tensor σ ji without the electric field and the new pseudo-vector σ × i read: where the latter may be expressed in matrix form as:

Compatibility equations
Again, two compatibility equations must be considered due to the existence of mechanical and magnetic fields.For the former, the strain measures-energetically conjugated to σ ji and τ ji from the first and second (18)-may be obtained from the additive decomposition of motion shown in Fig. 3.Then, the total deformation process may be mathematically represented by the gradient of macro-displacements u j,i composed of: Fig. 3 Micropolar mechanics deformation process, combining mechanical deformations and dipole micro-rotations, as in Fig. 2 I A mechanical deformation given by γ ji .II A pure rotation of the microstructure given by the skewsymmetric tensor Θ ji = −e jik θ k .
Unlike the CCM model, for which there is only one strain measure, in the MM framework and as in some beam theories, there are two deformation measures: where χ ji is the gradient of micro-rotations.As observed, the CCM strain measure is recovered without micro-rotations.The magnetic compatibility equation may be obtained by applying Helmholtz's theorem to the last Maxwell law (2) without electric terms to get:

Material constitution
The piezomagnetic constitutive equations used in the numerical experiments of Sect.6 are presented as a set of four equations: where C i jkl , e ϕ i jk , μ i j and C * i jkl are the elastic fourth-order tensor, piezomagnetic third-order tensor, magnetic permeability second-order tensor and micro-elastic fourth-order tensor, respectively.Furthermore, a new fourth-order tensor denoted by Ci jkl depends on the coupling number, see [18], assumed to be zero in the present work for simplicity.
Assuming the nine-component Voigt notation of Table 1 and considering that the material is magnetised along the x 3 direction as shown in Fig. 5, these constitutive tensors may be expressed in matrix forms as in the Appendix 1.For C * i jkl , the two micromechanical material coefficients l t and l b are the twisting and bending characteristic lengths; see [17].

Boundary conditions
Denoting by ūi and φ the prescribed displacements and scalar magnetic potential; and by tC the prescribed Cauchy traction, both boundary conditions may be expressed as: Notice that there are no natural boundary conditions for the MST since they emerge from the body forces; see (4).

Finite element formulation
This section briefly describes the Finite Element (FE) formulation conducted to implement the theoretical MM formulation reported in Sect. 4. In particular, a three-dimensional, non-linear, and displacement-based FE equation set is developed.There are no geometrical non-linearities since small strain is considered, but the MST introduces material nonlinearities that quadratically depend on the magnetic field, as shown in (20) top.For this reason, the FE formulation has two peculiarities: It is based on three residuals: linear and angular momentum balances plus magnetic Gauss law, to be explicitly defined in (28).A complex-step formulation with a perturbation parameter h n is conducted, as explained in [12,16,35].

Weak forms
To obtain the weak forms, the balance equations ( 18) are multiplied by test functions δu i , δθ i , δϕ of the dofs and integrated over the whole domain Ω: Applying the divergence theorem to the first term on the left-hand side of each (25), the final weak forms become:

Discretisations and residuals
In the present FE formulation, eight-noded elements with standard Lagrange functions N a at node a are used for the discretisation.In particular, the dofs and their test functions may be expressed by an isoparametric formulation as: where I denotes a complex number.The residuals at each node a are obtained by introducing the discretisations of (27) in the weak forms of (26) to give: where the uppercase subscripts refer to Voigt's notation of Table 1 and the matrices B i J and B i at each local node a are:

Tangent matrices
The complex-step method easily allows for computing the multi-coupled tangent matrices.These matrices are the derivatives of the residuals R a (u b ) for node a in terms of a generic nodal unknown u b .Then, writing a Taylor series expansion about u b in terms of an imaginary change of each dependent variable at node n, one may write: which upon collecting terms, yields: For cases where the higher derivatives are well behaved and using small values h n = 10 −40 , all residuals and tangent matrices are computed to full numerical accuracy from: This formulation is implemented into the research code FEAP [36] by using one of its dummy elements; the present numerical code without micro-rotations is tested against solutions of the total MST formulation developed in [26] obtaining perfect agreement.Notice that the multi-coupled tangent matrix is non-symmetric due to the inclusion of the MST and, consequently, the FEAP command UTAN must be used to invert the tangent matrix.6 Numerical experiments This section presents four numerical experiments to study the influence of the skew-symmetric part of the MST on the magneto-mechanical response in a magnetostrictive material under magnetic fields.For this purpose and as in [26], a magnetostrictive rod of length 6 [mm] and diameter 2 [mm] is simulated.The rod is made out of the alloy Terfenol-D, with properties listed in Table 2; the order of magnitude for l b , l t is taken from [3] for the crystal grains of many magnetic materials.This alloy presents transversely isotropic planes x 1 -x 2 since the magnetisation M 3 is along the x 3 axis, see the Fig. 4; therefore, the planes containing x 3 are anisotropic.
The skew-symmetric MST plays an essential role in activating the rotations as long as two conditions are met: -The magnetic field H i interacts with an anisotropic plane, not with the transversely isotropic plane.-The field is oblique to this anisotropic plane.
The idea is to prescribe an H i that activates the skewsymmetric term σ × i of (20), generating couple stresses τ ji and therefore micro-rotations θ i .These rotations are closely related to the existence of the microstructure, introduced in the FE model by the l b and l t scale parameters in the third tensor of the Appendix 1.
For all numerical experiments, the rod geometry is discretised by a FE mesh composed of 3,840 eight-noded elements, including 6,300 nodal points with seven dofs, as shown in Fig. 4 left.The mechanical dofs u i and θ i are fixed at the rod bottom, the rest are free.Since the present is a displacementbased FE formulation, no tractions are prescribed due to the absence of mechanical forces or pressures.
To prescribe nonuniform magnetic fields for the first three experiments with the FE method, a discretised air domain with dimensions 3×3×6 [mm] wraps the Terfenol-D rod: the additional mesh of 2,560 elements and 4,200 nodes is represented in Fig. 4 right.The vacuum permeability of air is μ 0 = 1.256×10 −6 [N/A 2 ], and its mechanical dofs from the coupled FE of Sect. 5 are deactivated so that only one FE type is used though the analyses for both materials.
In Table 3, the magnetic boundary conditions for the cases in the next subsections are listed; if no prescription is applied, the word "free" is indicated.

Magnetic field in transversely isotropic planes
The quadrupole magnetic field reported in [5] is chosen for the first three experiments.In the laboratory, the quadrupole is created by four magnets arranged so that H varies with the radial distance in planes x 1 -x 2 , but the field is constant along the x 3 axis.The magnetic potential boundary conditions are sketched in Fig. 5 left and listed in Table 3, producing a field that curves from one vertical side to the contiguous one (centre figure): then, H 1 ≈ H 2 and H 3 = 0.The field is almost axisymmetric, as seen in the figure of the right for H (smoothed through the whole mesh by the SPSVERBa1 software).But the critical domain is the Terfenol-D rod, where the field is uniform.
The maximum generated H is 22 [kA/m], a value that could be obtained by a commercial magneto cell under an electric intensity of only 5 [A].  Figure 6 shows the resulting contour plots for the three micro-rotations θ i generated inside the rod.As observed in the figure and expected, the micro-rotations are strictly zero in the air (green level) since its mechanical rigidity is null.The order of magnitude inside the Terfenol-D is 10 −13 [rad], not zero but almost numerical noise.This result is logical since under the quadrupole the magnetic field is contained in the transversely isotropic planes x 1 -x 2 , and no rotations are activated as it will be demonstrated in (36).
The main shortcoming of the present MM formulation is the necessary CPU time: Table 4 shows the calculation times of three models for the current numerical experiment.The first row is for CCM without Maxwell stress tensor (CCMw/o), the second for CCM with total Maxwell stress (CCMw), and the third for the present GCM.
All simulations are executed in a dual-core Intel Core i5-2500T running at 2.3 [GHz].For the first two models, the FE element developed in [26] was used with the mesh of Fig. 4.
As observed, the CPU time is substantially higher for MM, mainly because the three micro-rotations are extra nodal unknowns.The dofs are the same for the CCM-w/o and CCM-w models; however, the CPU time differs due to the non-linearity of the MST, which requires one more Newton-Raphson iteration.In this context, SPSVERBa1 calculates for all models the value: The residuals' superscripts k + 1 and 0 refer to the current and the first iterations, respectively.The solution is reached when val is lower than the computer-dependent tolerance tol = 10 −16 .

Symmetric magnetic field in anisotropic planes
With the previous mesh, the micro-rotations are now studied with the boundary conditions of Fig. 7 left and of Table 3.Now H varies in the vertical planes but not in the horizontal ones.In this way, the magnetic field curves from ϕ = 0 [A] to all planes but primarily to the vertical ones crossing the anisotropic planes as in the centre figure.
In Fig. 8, the resulting rotations are plotted.The microrotation θ 3 is again only numerical noise (almost zero) since the only magnetic field to produce this component would be in the isotropic plane x 1 -x 2 (of the Sect.6.1 type), but as mentioned and demonstrated in Sect.6.4, these fields do not create micro-rotations.The other two components are six or seven orders of magnitude higher than those of Sect.6.1, thanks to the prevalence of the H 3 component; H 1 is also high but only at the edges of x 1 constant.The value of θ 2 is almost three times larger than that of θ 1 due to the magnetic field concentration on some of the edges.
In any case, the micro-rotations are still much smaller than expected for an H partially applied in an anisotropic plane.There are two reasons for these low values: i) the perpendicularity of H to the planes x 1 -x 2 in some zones of the rod, and ii) the symmetry of this field that partially cancels θ 1 and θ 2 since the field concentrations at edges of the plane x 1 = 1.5 equilibrates the effect of the other concentrations at x 1 = −1.5, see the Fig. 7 right.3, the magnetic field curves from the ϕ = −10 [A] face towards all others in the whole domain.

Nonsymmetric magnetic field in anisotropic planes
In Fig. 10, the three rotations are represented: θ 1 is still not activated, but θ 2 is already a measurable 0.01 • , and again θ 3 is practically zero.The reason for θ 2 > θ 1 is that most of the magnetic field curves in the x 2 -constant vertical plane to the external face.In addition, the micro-rotation concentrates close to the face under ϕ =−10 [A] for the same reason, and any other rotation does not cancel it.
In the central figure of the highest component θ 2 , only the Terfenol-D part of the mesh has been represented to appreciate that the micro-rotation of the solid is uniform in horizontal planes and non-linear to fulfil the zero boundary condition at the rod base plus the free movement at the top.Although appreciable, these micro-rotations are still small since the magnetic field is primarily parallel (not oblique) to the vertical anisotropic planes in the central part of the rod.
Figure 11 presents two distributions: the left one relates the intensity of the magnetic potential up to a value of -50 [A] (achievable with commercial magnets) against the maximum θ 2 of Fig. 10 centre in absolute value.The trend is approximately quadratic, producing a substantial 25 increase in rotations for only five in magnetic potential.
As commented before, the micro-rotations depend on the micro-scale parameters l b = l t .For this reason and to study their influence on θ 2 , the right Fig. 11 plots the order of magnitude of the maximum micro-rotation versus the order of magnitude of l b .These rotations range from a non-physical 10 5 for very low l b to an almost zero 10 −10 [rad] for very high l b ; the real values given in [3] range from 10 −4 to 10 −6 .
To understand the linearity of the result, consider the equations that govern the macro-mechanical bending of a beam: with M b , m b , E, I , and φ denoting bending moment, moment source, Young modulus, second moment of area, and crosssection rotation in any plane, respectively.The integration of φ does not conserve constants due to the "cantilever"-type boundary conditions of Fig. 4.
For the current MM analogy, m is closely related to the skew-symmetric part of the MST and I to the bending scale length l b , which can be interpreted as the inertia of the micro-continuum.Consequently, the maximum of the three rotations from Fig. 7 (right) inversely depends on l b as in (34), that is, the lower its value, the higher the micro-rotations and the relationship is linear.

Uniform oblique magnetic field in anisotropic planes
The objective of the last numerical experiment is to analyse the influence of a completely oblique magnetic field on the generation of substantial micro-rotations.The importance of this study is that, for real materials, the initial magnetisation direction M 3 of the grains can rotate a different amount under the presence of a variable H, provoking the domain As observed in the left Fig. 12, the numerical experiment considers the Terfenol-D rod magnetised along x 3 under the action of an oblique and constant magnetic field function of the variable angle α.A vertical H 3 = 22 [kA/m] is first considered; to study the effect of different directions of the magnetic field in the plane x 2 -x 3 , a rotation with a variable angle α is calculated: Although the present FE is not mixed and, therefore, first derivatives cannot directly be prescribed, thanks to the fully nonlinear formulation, the calculated H1 , H2 , H3 are introduced in the SPSVERBa1 input as Neumann boundary conditions.Then, they can be automatically assigned to the nodal values of the magnetic field during the first iteration.
Figure 12 right plots the calculated maximum microrotation inside the rod in absolute value versus α.As observed, the micro-rotations are null for both α = 0 • and α = 90 • since the applied magnetic field is parallel (almost as in Sect.6.2) and perpendicular (as in Sect.6.1) to M 3 .Increasing α from zero, the micro-rotations increase to a maximum of 0.0665 [rad] for α = 45 • .
The explanation of this distribution can be found in the second (18), which can be expanded to: where the diagonal entries τ ii can be regarded as MM torsional moments and the off-diagonal τ i j as MM bending moments.Consider first that B i is approximately proportional to H i (u j,k is very small); then the second and third entries of the second vector at the left-hand side of (36) are zero (in this experiment ∀α → H1 = 0 and therefore B 1 ≈ 0).Since from  The right-hand side previous equation is maximum at α = 45 • and zero at α = 0 • or 90 • (vertical or horizontal H respectively), values that explain the maximum and minimums of Fig. 12 right.These two values also explain the nullity of the micro-rotations of Fig. 6 and, in part, the low values of Fig. 9, since for these angles, the second vector of the (36) right-hand side is nil.In these two experiments, the three equations give the trivial solution τ ji = 0, and without source, all micro-rotations are zero.
The discussion of ( 36) is qualitative but not quantitative: the difficulty lies in finding the distribution of the rotations inside the Terfenol-D rod, for which the present numerical method is necessary.Figure 13 shows FE contour plots (deformed mesh with zoom ×50,000) generated with the boundary conditions of Fig. 12.In the left figure, a nonlinear distribution is observed due to the clamped fixation of the rod at the bottom end.The same non-linearity is observed for the τ 31 component of the couple stress in the central figure, with a non-negligible reaction of 9.1 [N•m] at the fixed end; the deformed mesh shows a global bending of the rod, a shape that cannot be produced with, for instance, the CCM model.
To quantify the influence of the MST skew-symmetric part, the Fig. 13 right shows a ratio between the Frobenius norm of σ E M [ ji] and that of σ E M ji : The figure shows that the ratio is approximately constant since the applied magnetic field is also constant in the bar.
The MST skew-symmetric part is significant given that it generates about 6% of the total.One interesting outcome of the present study is that the mentioned "domain switching" phenomena could be explained with the micro-rotations and the influence of the magnetic dipoles' rotation on the magnetic material properties.These investigations will be studied in future work.

Concluding remarks
The present article has presented a theoretical and numerical formulation based on the finite element method to investigate the influence of the non-symmetry of the Maxwell stress tensor.This tensor contains a skew-symmetric contribution necessary for the study of rotations in non-isotropic materials; physically, these local rotations arise from the conservation of angular momentum.
In the framework of classical continuum mechanics, the total stress tensor formulation allows for calculating displacements of the center of mass at each material point but does not calculate rotations.On the contrary, extended formulations such as the current generalised continuum mechanics take into account the skew-symmetric part of the Maxwell tensor, which generates local rotations and couple stresses.
Numerically, a three-dimensional finite element formulation with seven degrees of freedom (three displacements, three micro-rotations, and the magnetic scalar potential) using a complex-step approach for the tangent matrices has been formulated and implemented in the research finite element code FEAP.
Four numerical experiments have been studied to extract the following consequences: -The skew-symmetric part of the Maxwell stress tensor generates both micro-rotations and couple stresses that depend on i) scale parameters (microsize-dependency) closely related with the internal structure of the magnetic materials and ii) orientation of the prescribed magnetic field.
-In the classical continuum mechanics framework, the scale parameters of i) are null and, therefore, there are no micro-rotations.-The dependency of ii) could explain the highly nonlinear domain switching mechanisms in magnetostrictive/electrostrictive materials since the rotations of the internal grains change the magnetisation/polarisation of materials.-For a magnetic field of 22 [kA/m] and for Terfenol-D material, the order of magnitude of these micro-rotations is up to 0.066 radians and its corresponding couple stress −9 [N•m].
In summary, the skew-symmetric part of the Maxwell tensor is a relevant 6% of the total.Therefore this part should be considered in modern and sophisticated electro-magnetomechanical devices.
Funding Funding for open access publishing: Universidad de Granada/CBUA.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A constitutive matrices
This appendix reports the matrix forms of the constitutive equations used in the present work.

Fig. 1
Fig. 1 Domain Ω = Ω E M ∪ Ω M ∪ Ω O with boundary Γ and outward normal n i .The domain is under electromagnetic f E M i Consider a domain Ω = Ω E M ∪Ω M ∪Ω O and its boundary Γ , containing electromagnetic Ω E M , mechanical Ω M and other Ω O matter, as shown in Fig. 1.The present section states the balances of linear and angular momentum considering electromagnetic forces f E M i

Fig. 4 to 3 Table 3
Fig.4 Finite element mesh of the Terfenol-D rod (left) for all numerical experiments; finite element mesh of air (right) for numerical experiments 1 to 3

Fig. 5
Fig. 5 Quadrupole magnetic field.Boundary conditions of magnetic scalar potential ϕ in [A] applied to external vertical sides (left); direction of the magnetic field lines (centre); contour plot of the resulting magnetic field module in [A/m] (right).External dimensions 3×3×6 [mm]

Fig. 6 Fig. 7 Fig. 8
Fig. 6 Terfenol-D rod wrapped by air under the nonuniform magnetic field of Fig. 5. Contour plots of micro-rotations in radians generated by the skew-symmetric part of the Maxwell stress tensor

Fig. 9 Fig. 10
Fig.9 Modified quadrupole magnetic field.Boundary conditions of magnetic scalar potential ϕ in [A] applied to external sides (left); direction of the magnetic field lines (centre); contour plot of the resulting magnetic module in [A/m] (right).External dimensions 3×3×6[mm]

Fig. 11 Fig. 12
Fig.11For experiment of Fig.9: left: absolute value of maximum micro-rotation in plane x 1 -x 3 versus prescribed magnetic potential; right: exponent of 10 for micro-rotation versus same exponent for microstructure scale parameter

Table 1
Tensorial notation with double indices (•) i j and equivalent Voigt notation with single indices (•) I

Table 2
[33]enol-D material properties at |H| = 100 [kA/m] from[33]: C i j stiffness, e ϕ i j piezo-magnetic coupling, μ i j magnetic permeability and l b = l t micromechanical characteristic lengths

Table 4
Calculation time (with respect to CCM-w/o) to solve the numerical experiment of Fig. 5 with three models.Number of Newton-Raphson iterations and converged residual norm for each experiment