Multiphysics and Thermodynamic Formulations for Equilibrium and Non-equilibrium Interactions: Non-linear Finite Elements Applied to Multi-coupled Active Materials

Combining several theories this paper presents a general multiphysics framework applied to the study of coupled and active materials, considering mechanical, electric, magnetic and thermal fields. The framework is based on thermodynamic equilibrium and non-equilibrium interactions, both linked by a two-temperature model. The multi-coupled governing equations are obtained from energy, momentum and entropy balances; the total energy is the sum of thermal, mechanical and electromagnetic parts. The momentum balance considers mechanical plus electromagnetic balances; for the latter the Abraham representation using the Maxwell stress tensor is formulated. This tensor is manipulated to automatically fulfill the angular momentum balance. The entropy balance is formulated using the classical Gibbs equation for equilibrium interactions and non-equilibrium thermodynamics. For the non-linear finite element formulations, this equation requires the transformation of thermoelectric coupling and conductivities into tensorial form. The two-way thermoelastic Biot term introduces damping: thermomechanical, pyromagnetic and pyroelectric converse electromagnetic dynamic interactions. Ponderomotrix and electromagnetic forces are also considered. The governing equations are converted into a variational formulation with the resulting four-field, multi-coupled formalism implemented and validated with two custom-made finite elements in the research code FEAP. Standard first-order isoparametric eight-node elements with seven degrees of freedom (dof) per node (three displacements, voltage and magnetic scalar potentials plus two temperatures) are used. Non-linearities and dynamics are solved with Newton-Raphson and Newmark-β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} algorithms, respectively. Results of thermoelectric, thermoelastic, thermomagnetic, piezoelectric, piezomagnetic, pyroelectric, pyromagnetic and galvanomagnetic interactions are presented, including non-linear dependency on temperature and some second-order interactions.


Introduction
Many modern technological devices, in particular electronics and sensors/actuators, extensively use multi-coupled materials characterized by the interaction of four fields: thermal, mechanical, electric and magnetic. These materials are increasingly subjected to sophisticated manufacturing processes that include miniaturization for microdevices. In addition, they are used in applications under very fast phenomena such as ultrasound waves or pulsed operation modes. Under these extreme conditions and given the relatively low-strength nature of practical active materials, stresses can reach limit values and produce total or partial failure.
Before exposing the objectives of the present work and in order to provide an idea of the importance of these materials, the following paragraphs list some interactions and their applications in modern technologies.
Piezoelectricity or electro-mechanical interaction, was discovered by brothers Pierre and Jackes Curie in 1880 (''direct'' effect) and by Lipmann in 1881 (''inverse'' or ''converse'' effect). Piezoelectric materials are used as sensors (e.g. for pressure), actuators (high precision positioning devices), in sonar applications (ultrasonics, hydrophones), in energy conversion (piezotransformers, energy harvesters); viz. [78] for more applications. Figure 1 shows a piezoelectric sensor fabricated with polarized ceramics. Currently, a new generation of piezo-elastomers is used in many applications. Piezomagnetism is the magnetic counterpart of piezoelectricity. The first experimental observation of piezomagnetism was made in 1960 by Borovik-Romanov, [19]. Piezomagnetic materials are also applied in energy conversion and medical applications [82].
Magnetostriction or magneto-mechanical interaction, was first reported by Joule in the early 1840's, observing the change in length of iron bars when they were magnetized. The first important application of these materials occurred during World War II to build transducers for sonars. Currently, they are used in active vibration control, control surface deployment and energy harvesting, viz. [4] for a full revision on applications. Figure 2 shows an actuator composed of the magnetostrictive Terfenol-D inside a solenoid that generates a magnetic field. Electrostriction is the electric counterpart of magnetostriction, also applied in sonars and sensors/actuators for small displacements, viz. [37] for a full review of applications. Note that both electrostriction and magnetostriction are second-order interactions related to the first-order ones described for piezoelectric and piezomagnetic interactions, respectively.
Pyroelectricity denotes the interaction between heat and electric fields. It is observed in certain materials that generate electric voltage when they are heated or cooled. According to [38], the first reference to pyroelectric interaction was reported by Theophrastus in 314 BC when he wrote that ''lyngourion'', probably the mineral tourmaline, had the property of attracting straws and bits of wood. Much later, the naturalist Linné related this pyroelectric property to electricity; he called tourmaline the ''electric stone''. Nowadays, a new generation of artificial pyroelectric materials in the form of thin films has been developed. For instance, lithium tantalite is used to create small-scale nuclear fusion or ''pyroelectric fusion'', viz. [95]. Pyromagnetism is the magnetic counterpart of pyroelectricity. There exist many sensors that use the pyromagnetic interaction, for example to detect radiation [144] or to measure temperature [107].
Thermoelectricity or interaction between heat and electric fluxes, was observed by Seebeck and Peltier even before the quantitative Ohm's law was formulated in 1855. Three different transport effects are present in thermoelectricity: Seebeck, Peltier and Thomson. In addition, the Ohm and Fourier laws inherent to electric and thermal fluxes are also present. Thermoelectric parts are used as coolers (cooling electronic devices, refrigeration, air conditioning), for power generation (energy harvesters, photovoltaic cells), as energy sensors (detection of water condensation and fluid flow, infrared thin film, cryogenic heat flux sensors), viz. [125] for more applications. A thermoelectric behavior has recently been observed in the Oriental hornet (Vespa orientalis) cuticular abdomen, [46]. Figure 3 shows several thermoelectric material devices and the hornet.
Several scientists related with the discovery of the interactions present in this work are shown in Fig. 4. Linné (1707-1778), Swedish botanist, physician and zoologist, discovered the pyroelectric interaction. Fourier (1768-1830) was a French physicist known for his works on heat transfer. As described before, Seebeck (1770-1831) discovered the thermoelectric effect in 1821. Peltier (1785-1845) was a French physicist who discovered the effect named after him in 1835. Ohm (1789-1854) was a German physicist known for his contribution to electric current fluxes. Joule (1818-1889), an English physicist author of the namesake effect. Fick (1829-1901), German physiologist, started to study mathematics and physics, but  Hall  was an American physicist who conducted thermoelectric research at Harvard University. The Curie brothers , French physicists, discovered the piezoelectric effect, as mentioned previously. Nernst (1864-1941), German physical chemist and physicist, is known for the formulation of the third law of thermodynamics.
The previous descriptions are just examples of interactions that can occur in a material subjected to up to four fields. Therefore, it is important to describe the problem by means of an introduction to multi-coupled formulations.

Multi-coupled Formulations
Every cause has an effect; in Physics, cause and effect are described by intensive and extensive variables, respectively (see Sect. 3.1). The interaction of both types of variables is expressed through phenomenological equations, generally non-linear and written in terms of certain coefficients denoted as material properties: observables that can be measured. Figure 5 shows a conceptualization of the phenomenological equations, mathematically representing interactions. Traditionally, these interactions have been classified into two groups: • Equilibrium interactions, EI • Non-equilibrium interactions, NEI Material properties are tensors that relate intensive and extensive variables, describing linear or non-linear interactions. These properties are classified into three groups, [101,141]: The interaction among variables is not straightforward, therefore the properties often depend on frequency and on the movement of domain-walls; from a macroscopic point of view the interaction is then non-reversible or dissipative. Then, the effect associated to a given cause not only depends on the cause itself, but also on the whole history of states that the material has undergone to reach the current state: important memory effects exist in the definition of the phenomenological equations. Polymeric photovoltaic cell. A photovoltaic-like behavior has been observed in Vespa orientalis hornet. Meso Peltier cooler. Hornet taken from [77], rest from manufacturers Fig. 4 Great scientists who observed principal and coupled interactions among elastic, electric, magnetic and thermal fields. These and other following pictures of scientists taken from [133] Multiphysics and Thermodynamic Formulations for Equilibrium and Non-equilibrium Interactions… 537 In Fig. 6 a schematic representation of the interaction between extensive and intensive variables is shown. Three regions can be differentiated: (a) Linear and reversible; first-order material properties suffice to describe the interaction. (b) Non-linear but still reversible; the region can be described with second-order material properties. (c) Non-linear and irreversible due to domain-wall properties, represented by a hysteretic branch.
The following subsections introduce equilibrium interactions (EI), non-equilibrium interactions (NEI) and their related material properties.

Equilibrium Interactions
The EI are formulated in thermodynamically reversible processes. These processes can be approximated by a succession of equilibrium states, therefore, EI is studied by equilibrium thermodynamics.
The present intensive variables in EI are stress, electric and magnetic fields and finally temperature; the respective extensive ones are strain, polarization and magnetization and finally entropy, see Table 5. Following the convention of Fig. 5, these variables are represented by rectangles (intensive) and circles (extensive) in the Heckmann diagram of field interactions shown in Fig. 7; the first-and second-order properties are represented by upward and downward vertex triangles. Only first-order properties of mechanical, thermal or electrical fields and their interactions are usually represented in the diagram, viz. [5,6,139,141]; in the present work the magnetic field is also included, see Sect. 1.1.3 for discussion. Therefore, a pair of conjugate variables (intensive and extensive) directly relates each of the four fields. According to [101,141], we count in total 28 interactions for these four fields: • 4 principal interactions relating conjugate variables. • 24 coupled interactions relating all intensive and extensive variables.
The properties' definitions are listed in Tables 1, 2 for the principal and coupled interactions, respectively. The four principal and for clarity only 12 (including the direct and converse) coupled interactions are represented in Fig. 7.
As mentioned previously, second-order properties represent reversible but non-linear interactions and in general their inclusion allows for an almost exact modelling of the phenomenological equations. These properties are weak correction terms of the first-order properties but in some particular cases may be dominant, for instance under high field strengths in electrostriction or magnetostriction.

Non-equilibrium Interactions
The NEI are studied by non-equilibrium thermodynamics or, in other words, thermodynamics of irreversible processes. In the framework of this type of thermodynamics, the intensive (voltage, concentration, temperature gradients) and extensive (electric, thermal, mass fluxes) variables are denoted as driving forces and fluxes, respectively.
As before, they are represented by rectangles and circles in Fig. 8. Also, as in Sect. 1.1, material properties are classified into first-order, second-order and domain-wall, and the    Multiphysics and Thermodynamic Formulations for Equilibrium and Non-equilibrium Interactions… 539 three regions shown in Fig. 6 remain valid. However, now all processes are irreversible (return arrows should be deleted) since the entropy of a system out of equilibrium always increases, see Sect. 3.1. In addition the interpretation of these domain-wall properties requires extended non-equilibrium thermodynamics (defined in Sect. 2.1) and, according to [57], still represents a challenge for the theoretical physicists. For first-order properties, three pairs of conjugate variables with a total of 15 interactions are considered: • Three principal conductivities between fluxes and gradients • Six coupled power properties between fluxes and gradients • Three interactions between fluxes • Three interactions between gradients The last two are not shown in the diagram for clarity. Tables 3 and 4 list the three principal conductivities and the three power properties, respectively. In the second item, the number six refers to direct plus converse, a concept often not used in NEI. The interactions between fluxes (third item) do not require any additional material properties since these fluxes will be included in the balance equations. For example, an electric flux increases the thermal flux due to the Joule effect without the necessity of using any additional property. The interactions between gradients (fourth item) are in general used only for measurement of first-order properties. For instance, the Seebeck coefficient (see Sect. 4.3.1) is typically measured forcing the electric flux to zero and comparing temperature and voltage gradients.
The presence of magnetic fields or mechanical stresses modifies classical NEI, requiring changes in the material properties' definitions and leading to new and more complicated interactions, such as piezoresistance, magnetoresistance, galvanomagnetism and thermomagnetism. Figure 8 shows a Heckmann diagram for first-order NEI analogous to that of Fig. 7 for EI. Again, the first-order properties are represented by triangles, see Tables 3, 4 for notation. For clarity, the interactions caused by mechanical stresses and magnetic fields are not represented although they will be included in the formulation through the Maxwell tensor plus galvano-and thermomagnetic effects (see Sects. 1 Tables 3, 4   Table 3 First-order non-equilibrium principal properties (conductivities) Summarizing, an active material simultaneously affected by mechanical, electric, magnetic and thermal fields can be subjected to 43 interactions requiring 40 first-order material properties to describe these interactions. Additional complications could arise: some of these properties can depend on the magnitude of some fields and as mentioned, if hysteretic behaviors are present, the properties are time dependent with important memory effects.

Galvanomagnetic and Thermomagnetic Interactions
Special cases studied by NEI with a magnetic field transversal to the fluxes are named galvanomagnetic and thermomagnetic interactions. They are a consequence of Lorentz forces that a magnetic field exerts on moving charges; viz. [24,33,80,101,127,141] for description and Sect. 3.3.4 for formulation of the forces. These interactions are widely present in transducers and sensors [135]; engineering and medical applications (biosensors) [68]; optimization of thermoelectric devices [101]; chemical processes and in cell-stimulant interactions [68], among others.
In spite of its practical importance, in the theoretical thermodynamics literature the magnetic field strength (or equivalently the magnetic field density or magnetic induction B) is not usually included in the multi-coupled formulation since its presence leads to the treatment of material properties as tensors, instead of scalars. But in the context of numerical methods the formulation and operation of these tensors do not present any problem.
Since the main objective of this article is to study multiphysics phenomena through numerical methods, a special formulation is developed in Sect. 4.3.2: physically the thermal and electric fluxes interact with the magnetic field due to the thermal and electric energies that moving charges carry. Therefore, in the presence of a magnetic field: • Galvanomagnetic interactions involve voltage gradients and electric fluxes, Fig. 9 left.
These interactions depend on the mechanical anisotropic level of the media; second-order interactions such as piezoresistance or magnetoresistance emerge in anisotropic materials due the dependency of properties on the intensity of B, [141]. In any case, anisotropic mechanical media are not considered in the present work for NEI, although the material can be magnetically anisotropic.
Considering under NEI only isotropic media and neglecting second-order interactions, four first-order transverse interactions due to the presence of B are defined. Assuming this induction along direction x 3 as in Fig. 10, these effects are • Hall electric flux j along x 1 (additional cause) induces voltage gradient rV along x 2 (effect). • Ettingshausen j along x 1 induces temperature gradient rT along x 2 . • Righi-Leduc rT along x 1 induces rT along x 2 .
• Nernst rT along x 1 induces rV along x 2 .
The present work studies all the interactions represented in both Heckmann diagrams, except those involving concentrations (e.g. the last two in Table 4). Concentration is not included since active materials are usually (although not always) solids or polymers. Therefore, effects such as porous media and chemical reactions are not considered. The spirit of the article is in any case to formulate the physics and numerical implementation so that any effect, known or new, can be included in the future.  Fig. 9 Physical magnitudes involved in the galvanomagnetic and thermomagnetic interactions. Principal and coupling properties (triangles) defined in Tables 3, 4 2 Literature Review The literature review has been divided into two parts: theoretical and Finite Element Method (FEM) formulations.

Theoretical Formulation Review
From a thermodynamic point of view, basically there are two theoretical formalisms (frameworks to study interactions) for the proposed multi-coupled formulation: Non-equilibrium thermodynamics used in [23,24,80] considers thermal, electric and magnetic fields only for non-polarizable/non-magnetizable materials.
The mentioned works were extended in [32,33] to incorporate the mechanical field as well as material polarization and magnetization. Dielectric/magnetic relaxations, known as Debye relaxations and interpreted as irreversibilities due to microscopic polarization/magnetization interactions, were not fully studied. The dielectric relaxation was investigated for isotropic materials in [123] and for anisotropic ones in [124]. The magnetic relaxation was investigated in [122] using the method of internal variables.
Rational thermodynamics was stated in [31,142] and transforms the thermodynamic concepts into a Continuum Mechanics framework.
Several other authors have used this formalism to obtain a multi-coupled formulation, viz. [42, 43,102] for example.
As discussed at the beginning of Sect. 1, miniaturized devices and high-frequency processes are increasingly applied in modern technologies. For these applications the local equilibrium hypothesis that will be stated in Sect. 3.1.2 is not valid, requiring newer and advanced theoretical formalisms; according to [83], the following methods can be used for the study of these cases: Extended non-equilibrium thermodynamics introduces state variables, see Sect. 3.1 for definition, and fluxes as independent variables (mixed formulation). Rational extended thermodynamics develops evolution equations for fluxes, although introducing supplementary Lagrange multipliers. Internal variable methods macroscopic representations of microscopic internal structures that are incorporated to the state variables. Hamiltonian formalisms generalization of the Poisson bracket formalism [79], expressing the evolution equations with two (total energy and dissipation) thermodynamic potentials called generators.
Using the first method, an extended approach considering thermal and electric fields without polarization was developed in [86]. The polarization effects were subsequently investigated in [35] while the mechanical field was incorporated in [87,90]. A full revision of extended nonequilibrium thermodynamics, its applications and theoretical developments can be found in [97][98][99].
Examples of multi-coupled formulations using rational extended thermodynamics and internal variable methods are presented in [44,122], respectively.

Finite Element Formulation Review
To the best of our knowledge, there has been no comprehensive attempt to propose a unified formulation for the four-field, multi-physics problem of active materials using the FEM. Most advanced publications study single interactions or phenomena; only very recently the issues of multi-physics and multi-coupling have been addressed in [137], but with micro-mechanical models. Their computational cost is very elevated and, for instance, design or optimization of devices would be very expensive if not impossible.
Separated states of art descriptions for the relevant interactions are developed in the following sections.

Piezoelectric Interaction
There are in the literature two main alternative FEM formulations for the piezoelectric problem: scalar and vector. The former uses four dof per node: three displacements and a scalar potential or voltage. A monolithic or fully coupled (coupling of the stiffness matrix) formulation was developed by [2], the first two authors of the present work

Fig. 10
Four first-order transverse galvanomagnetic (top Hall and Ettingshausen) and thermomagnetic (bottom Righi-Leduc and Nernst) non-equilibrium interactions. Causes represented by rectangles, effects by circles developed a similar one for plane strain in [110,114,128,129]; an iterative or stagger formulation was proposed by [48]. Starting from the formulation given in [2], several modifications have been introduced in the last decades: the assembled matrix was rescaled in [116] to avoid ill-conditioning while the remanent strain and polarization were included in [151]. The latter addition could be considered the first step to model the hysteretic behavior inherent to ferroelectric materials. According to [70], there are three different approaches to model hysteresis: • Thermodynamically consistent models, [72] • Micromechanical models, [92] • Models with hysteretic operators, [70] For the vector formulation the work developed by [81] uses six dof per node: three displacements and three components of the vector electric potential. The main difference between both formulations arises from the choice of the constitutive equations: the electric displacement is chosen as independent variable for the second with the disadvantage of loss of uniqueness. For this reason some gauging procedures were investigated in [131], incorporating the Coulomb gauge to the formulation developed in [81] using the penalty method. To solve the gauge problem, some authors have developed solutions based on edge vector interpolation functions, which include scalar parameters. The formulation is complicated by the fact that parameters appear on element edges, viz. [96]. The main drawback of the vector formulation is the large null-space in the discretized curl equation.

Magnetostrictive Interaction
As for the piezoelectric interaction, two FEM formulations exist, one based on the scalar (four dof) and another on vector (six dof) magnetic potentials. Furthermore, fully coupled and staggered formulations are reported in the literature.
For the scalar formulation, a fully coupled, 2-D and nonlinear FEM formulation was developed in [12]. The same authors extended the formulation incorporating the dynamic response in [11] and the eddy currents in [13]. An alternative 2-D staggered formulation was proposed in [59], assuming non-linearity only for the magnetic field. A fully coupled 3-D and dynamic FEM was formulated in [73]. The hysteretic behavior was modeled using thermodynamically consistent models in [71,76,84,85].
For the vector formulation, a fully coupled, 3-D, nonlinear and steady-state FEM was developed in [121], solving the non-linearities by the Newton-Raphson algorithm and symmetrizing the whole matrix. Similar 2-D FEM models were reported in [14,52,146]. The last reference inserted experimental magnetostrictive curve data into the tangent stiffness matrix. An alternative fully 3-D and non-linear FEM formulation, including the Maxwell stress tensor, was reported in [113]. This formulation was extended in [74,147], including the material non-linearity, i.e. the dependence of the properties with the magnetic field.

Thermoelastic Interaction
For the thermoelastic interaction description, it is necessary to develop a brief historical introduction about the theoretical formulations. In principle, the classical thermoelasticity is not compatible with the physical observations since intrinsically suffers from two shortcomings: 1. One way coupling: the mechanical field does not influence the thermal field although the thermal does with the mechanical. 2. The heat equation is simply parabolic: predicts infinite speed of propagation for the heat wave.
The first shortcoming was solved in [15] with the theory of coupled thermoelasticity. For this purpose the entropy was balanced, introducing a dissipative or entropy production term that fully coupled the governing equations. Numerically, this solution was implemented in several FEM articles. For example, [25,132] developed a 3-D isoparametric formulation and an eight-node element based on the Reissner-Mindlin plate theory, respectively. In a simpler manner and using only the one-way uncoupled thermoelasticity, staggered FEM formulations have been implemented, e.g. viz. [103]. These formulations use two steps: temperature distributions are obtained solving the uncoupled heat conduction problem and then solving the elastic problem with thermal stresses calculated from the strains caused by the temperature distribution.
The second shortcoming was theoretically addressed in [26] where the Fourier law is modified by means of an empirical parameter named relaxation time that results in a hyperbolic model. In the Continuum Mechanics community, this Cattaneo model is denoted as second sound (SS). Numerically, the influence of the SS in the thermoelasto-dynamic behavior of continuum bodies has been studied by several authors using the FEM. The difficulty of the SS modelling is the time integration scheme since numerical oscillations appear. Explicit finite differences were used in [136], enabling an improved physical interpretation. In [7], the temporal discretization was solved using discontinuous and continuous Galerkin methods (mixed method for the time integration). However, numerical oscillations appeared and in [8,9] developed a stabilization method. The Newmark-b algorithm, with optimized time steps and algorithm parameters, was used in [140,150,152].

Thermoelectric Interaction
A steady-state and non-linear 3-D FEM formulation, considering temperature-dependence properties, was reported in [49,111,115] and implemented into the research FEM code FEAP. Dynamic thermoelectric elements were implemented into the FEM commercial codes ANSYS and COMSOL in [3,39], respectively. These implementations did not consider temperature-dependent properties, but included a standard interface element to model heat convection. In addition, the commercial codes permit the study of one-way mechanical responses, feature that was used by [30,47,60]. The elasto-thermoelectric interaction, considering temperature-dependence properties, was implemented into COMSOL in [65, 66].

Other Interactions
In the specialized literature, there are works that study other interactions such as pyroelectric, thermomagnetic and galvanomagnetic.
For pyroelectric interactions, two 3-D elasto-thermoelectric FEM formulations are developed in [55,143]. The first uses standard isoparametric elements, while the second adds additional internal dof to improve the numerical results in laminated pyroelectric plates.
Several analytical solutions for thermomagnetic and galvanomagnetic interactions are available in the literature, see for example [27,36]. However, these interactions are very sensitive to geometry and to material properties, [68,105,127], and in addition they induce strong distortions of the electric and thermal fluxes. These difficulties imply that the analytical solutions are only able to study simple geometries and often are not used in sophisticated applications. For these reasons, several numerical formulations have emerged in the last two decades. The finite difference method was applied in [104,105] and the FEM was used to study the Hall effect in [21]. The mentioned FEM formulation is very simple and not fully coupled, consisting on an electric element with anisotropic electric conductivity and magneticdependency of the conductivity tensor introduced as data.

Outline of Continuum Physics
This section presents an outline of Continuum Physics, introducing the necessary concepts and equations to develop the multi-coupled governing equations composed of balance plus transport equations and of boundary conditions. Section 3.1 reviews the thermodynamic formalisms from Classical Thermodynamics (or Thermostatics) to one of the formalisms used in the present work, Non-Equilibrium Thermodynamics.

Thermodynamics
Thermodynamics is a physical science that studies transformations of energy in all its forms. The term was born in 1854 when Kelvin stated the concise definition Thermo-dynamics is the subject of the relation of heat to forces acting between contiguous parts of bodies, and the relation of heat to electrical agency.
One of the founding articles was a 300-page paper published by Gibbs, Fig. 11, in 1882: On the Equilibrium of Heterogeneous Substances. Before discussing the thermodynamic formalism and in order to clarify the formulations developed in this section, several definitions are introduced.
A thermodynamic system is a portion of matter with domain X and boundary C. The union of the system and its surrounding X 1 constitutes the thermodynamic universe, sketched in Fig. 12. There are three types of thermodynamic systems  Fig. 12 Thermodynamic universe composed of thermodynamic system: X (domain), boundary C and surroundings X 1 • Open system exchanges matter and energy with its surrounding. • Closed system exchanges energy but not matter.
• Isolated system does exchange neither energy nor matter.
The state variables S are a set of path-independent variables that describe the ''state'' in which the thermodynamic system is. Two types of state variables can be defined • Intensive I depend on the external conditions imposed on the system. • Extensive E depend on the amount of matter in the system.
As mentioned in Sect. 1.1, intensive and extensive variables represent causes and effects.
A thermodynamic process represents the energetic evolution of a thermodynamic system from an initial state to a final one. There are two types of processes • Reversible: continuum sequence of equilibrium states (idealized process) • Irreversible: any non-reversible process The discipline of thermodynamics usually encompasses different formalisms; according to [83], three can be defined Equilibrium Thermodynamics (ET) [24], studies macroscopic properties of systems in mechanical, thermal, electric, magnetic and chemical equilibrium.
This formalism should be denoted as thermostatics since it refers to time-independent and spatially homogeneous systems, but the term thermodynamics is generalized for historical reasons.
Non-Equilibrium Thermodynamics (NET) or thermodynamics of irreversible processes [33], deals with processes that depend on time and on spatial coordinates. It is based on the local equilibrium hypothesis, stated in the following subsections. Extended Non-Equilibrium Thermodynamics (ENET) [69], studies systems in which the local equilibrium hypothesis is not valid, using thermodynamic mixed formulations.
The first two formalisms will be described in the following subsections and the third is the subject of another ongoing work.
The first principle or energy conservation law states that the energy of the universe is constant; mathematically where E T is the thermal internal energy (sum of elastic, kinetic and potential) of the system and Q, W the heat and work supplied by the system, respectively. The symbol d denotes dependency on the path: the related variables are not state but path-dependent variables. The other derivative d) represents an exact differential and is used for state, path-independent, variables. Therefore, it can be established that dQ, dW represent equilibrium (reversible) and dQ, dW non-equilibrium (irreversible) parts. The first law is valid for reversible or irreversible processes and for both EI and NEI. The second principle or entropy law introduces the concept of entropy S and states that the entropy of the universe never decreases, DS ! 0.
With this concept, the first law can be rearranged as The first two terms of the right side are reversible and described by ET in the next subsection; the second pair are irreversible and studied, in addition to the first, by NET in Sect. 3.1.2.

Equilibrium Thermodynamics
ET is based on the first and second fundamental laws; the zero and third laws are not necessary in the present work.
In this formalism it is assumed that no irreversibilities are present, dQ ¼ dW ¼ 0, that is, the intensive variables must be at equilibrium. Therefore, it can be established that for reversible processes in which the entropy is an exact differential, the second law is where T e is the equilibrium (also called thermodynamic) temperature. The Gibbs equation for EI is obtained combining the first law (1) and the expression (3) where the reversible work dW is expressed as the sum of the products of extensive and intensive variables for fields other than thermal and l 0 is the vacuum permeability. This field does not contribute to dW but will do to the entropy dS in the left hand side.
Each of the partial derivatives inside the matrix represent principal and coupled material properties (Tables 1, 2) and each is calculated with the other independent variables (vector in the right hand side) maintained constant. The expression (5) can be rewritten with the symbols of Fig. 7, giving the following expression Note that intensive (I , rectangle) and extensive (E, circle) variables are mixed in the previous matrix relation; the independent variables in the right hand side are located there for a proper FEM implementation and for the use of simple compatibility expressions (see Sect. 4.2.2) to introduce the basic nodal, zero-derivative unknowns.

Non-equilibrium Thermodynamics
NET transforms the thermodynamic laws into continuum or local forms (balance equations), assuming the continuum and the local equilibrium hypotheses. In this sense, NET may be interpreted as a Continuum Physics formalism: its state variables depend on time and on spatial coordinates. Before expressing the laws in continuum forms, the reference system and the two hypotheses are defined.
Since irreversibilities are present, (3) is now rewritten to where T n denotes the non-equilibrium (also called conduction) temperature, see Sect. 4.1.4 for the disambiguation of both temperatures. The presence of the pathdependent variable dQ in the heat exchange requires the definition of T n to ensure that entropy is a state variable. Two descriptions to study the motion of the system are often used: Lagrangian and Eulerian (authors in Fig. 13), [18,41]. The material or Lagrangian X ¼ fX 1 ; X 2 ; X 3 g > description refers to the behavior of a material point (commonly used in Solid Mechanics); the spatial or Eulerian x ¼ fx 1 ; x 2 ; x 3 g > one to spatial positions, commonly used in Fluid Mechanics. In the above the symbol ðÁÞ > denotes transpose.
The consideration of the two different descriptions implies the definition of two: position X, x and displacement U, u vectors, and time derivatives. These derivatives are related by where t denotes time, P x; t ð Þ a continuum property in Eulerian description and v the Eulerian velocity. The last term on the right side of (8) is called a convective derivative, since it is closely related to the particle motion inside the system.
In the remainder of the present work a Lagrangian description is used. Therefore and according to [106], x X, u U and the convective term is avoided, an obvious advantage.
The continuum hypothesis [41], assumes that the matter in the system X is continuously distributed and fills the entire X; mathematically where q m denotes the mass and q q the electric charges densities, the latter may be positive or negative. The total mass and electric charge contained in the incremental volume DX are Dm, Dq. When DX is greater than a certain critical DX Ã , the Continuum Mechanics is an appropriated mathematical model since both densities depend on spatial coordinates and on time and not on DX. Note that the validity of this hypothesis is closely related to the size of DX Ã and, therefore, to the critical length that will be defined in (10).
In NET, the local equilibrium hypothesis [83] permits one to rewrite (4) locally for any time and for any material point, solving the restriction that establishes the time-independence of the state variables (in principle valid only in ET). A physical interpretation of this hypothesis is related to the fact that each point is in a different equilibrium state. Exchange of physical quantities among different points (or equilibrium states) is possible; furthermore, the equilibrium state of each point changes over time.
The validity of the continuum and local equilibrium hypotheses is closely related to the Knudsen and the Deborah numbers, the last defined by Reiner (see Fig. 14) where Kn is the ratio between the microscopic parameter that represents the mean free path and the macroscopic length; De is also a ratio but between the microscopic equilibration or relaxation time and the macroscopic duration of the interaction, both inside each point. Although commonly accepted, there are several important applications for which these two hypothesis are not applicable • Kn ! 1, for micro-and nano-systems, thin films, superlattices, porous media, etc. • De ! 1, for ultrasonic propagation in dilute gases, polymers, superconductors, etc.
When the condition in the previous item holds, the ENET (not an object of the present work, viz. [108,109] for applications) has to be considered.
Assuming from now on the validity of the two hypotheses, any thermodynamic variable may be expressed in continuum form as P, so that The mass balance equation states that the mass inside the system remains constant over time. According to [41] the temporal variation of the mass is The energy balance equation states that the total energy density e tt in the universe is conserved since sources and/or sinks cannot be included, see Sect. 4.1.2 for additional discussion. As will be depicted in Fig. 21, e tt is the addition of mechanical, thermal and electromagnetic parts. Therefore, the evolution of energy inside the system must be equal to the total energy flux j tt through the boundary, second equality in next equation. The first expression represents the temporal variation of the total volumetric energy density; taking into consideration (12) the second expression is obtained. Finally the divergence theorem is applied to the third expression that represents the total energy flux where n denotes the outward normal. Finally, from the previous equation and [33] the local form of energy balance is where the supra-dot indicates exact derivation with respect to time, both for partial and exact differentials due to the absence of the convective derivative of (8).
The difference between the first law and the energy balance is that in the former the energy is expressed as thermal E T and in the latter as total e tt . For this reason, the first law can include source/sink terms that represent reversible exchanges between.
The entropy balance equation is equivalent to the second thermodynamic law (3) rewritten in the continuum form again through the divergence theorem. Considering that from (11) in NET q m ds ¼ dS 6 ¼ 0, and following similar steps as before Z where j r , r are the entropy flux and entropy production, respectively. Finally, the local form of the entropy balance is given by where the entropy production must be The transport equations (defined as constitutive for ET, see Sect. 3.1.1) are obtained expressing the entropy production as where j k F are fluxes (extensive variables) and F k driving forces (intensive variables). Both types are listed in Table 6 for NEI not considering concentration phenomena, see Sect. 1.1.2 and Fig. 8. NEI variables are denoted by ðÁÞ n to differentiate them from the EI ones. The electric field can produce two fluxes: j f for conductor media and E n for polarizable media. In what follows and for notational simplicity the supraindex f (meaning flux from free sources) will be dropped and j j f , without subindex is directly the electric flux.
As mentioned in Sect. 1.1, fluxes and forces are linearly related with a good approximation by the phenomenological equations where L kl are the first-order material properties, see Fig. 5. According to [33], these properties have three restrictions: 1. Due to material symmetry and according to Curie's law (viz. [141]), driving forces cannot have more elements of symmetry than the fluxes they produce. 2. Due to the continuum second law (16), the sign of the entropy production implies 3. Due to the time reversal of the microscopic governing equations, the Onsager-Casimir (authors in Fig. 15) reciprocal equations state The last expression indicates that the material property tensors must be either symmetric or antisymmetric. Throughout the paper, L kl will represent heat and electric conductivities, Seebeck thermoelectric power and Peltier properties.
Note that restriction 1 is not valid in the non-linear region of Fig. 6, according to [83].
The conditions given in (20) are equivalent to the ones stated by the axiom of admissibility in the Rational Thermodynamics formalism, viz. [142]. For this formalism, the axiom is mathematically represented by the Clausius-Duhem inequality.
Similarly to (6), the NEI transport equations based on Fig. 8

and (19) can be established as
Instead of being included in the previous equation, the interactions due to B will be decomposed into symmetric and antisymmetric parts in Sect. 4.3.2. Also, the interactions due to the mechanical field will be considered in the EI formalism due to the absence (simplification S1 in Sect. 4.1.3.) of viscoelastic effects

Classical Continuum Mechanics
In this section the governing equations of Classical Continuum Mechanics are presented. The present work assumes a small deformation theory based on two assumptions, [106] 1. Small displacements jjujj ( jjxjj 2. Higher-order terms neglected, implying that the mechanical compatibility equation is given by where the second order displacement gradient ru has been decomposed into its symmetric r s u and skew-symmetric r sk u parts. From the assumption of small strain jjSjj ( 1 and jjWjj 2 ( jjSjj, where W is the small rotation tensor.

Linear Momentum Balance
Consider a continuum system (or body in the Continuum Mechanics framework, Fig. 16) subjected to volume f and boundary t forces. The volume forces will in general include all possible types: gravitational, centrifugal, electromagnetic etc. The body is restricted by Neumann and Dirichlet conditions on boundaries C p , C u , respectively where u, t denote the prescribed displacements and traction. The traction includes the total (Cauchy-like), [91], stress tensor. This total stress combines the EI and NEI parts defined in Tables 5, 6: The balance is obtained in a way similar to that of (13) but with Newton's law, integrating all forces in X and again applying the divergence theorem yields Z Rewriting (24) in local form Notice that this equation does not have a counterpart in thermodynamics, theory that does not deal with momentum or forces. The total stress T will be non symmetric if electromagnetic coupling is present, therefore the conservation of angular momentum is not automatically satisfied. One way to solve this difficulty when the formulation is multicoupled as in the current work, is to force the sum of T and the coupling tensor to be symmetric, as will be described in Sect. 4.1.2.

Energy Balance
Assuming for now in Continuum Mechanics adiabatic processes dQ þ dQ 0 (no exchange of thermal energy between the system and its surrounding media), from (1) dE T ¼ dW þ dW implying that the total energy is equal to the work of the external forces, both reversible and irreversible, E T ¼ W. Therefore, the mechanical internal energy (including kinetic energy) is equivalent to the mechanical power, q m _ e U ¼ À _ W. From Continuum Mechanics, W is the sum of kinetic and elastic T : S energies; therefore the mechanical energy balance in local form is given by

Classical Electrodynamics
Classical Electrodynamics also called Classical Electromagnetism is a physical theory that studies the interactions between electric charges and currents. It was developed over the course of the nineteenth century, most prominently by Maxwell, Fig. 19.

Maxwell Equations
Maxwell equations are a set of four empirical and macroscopical equations that relate electric and magnetic fields.
In vacuum, where 0 denote the vacuum permittivity; q f q , j j f , q b q , j b are the free and bound sources (electric charge density, electric flux), respectively, and 0 is now a 3 Â 1 vector with zero entries. The free sources may exist inside the material, while the bound sources are generated by the interaction of an external field on the system. The first equation in (27) is the Gauss electric law ( Fig. 17) that relates the electric field E and its sources q q ¼ q f q þ q b q , the electric charge density defined in (9). As was done for stresses, the total electric field adds the EI and NEI parts, E ¼ E e þ E n . The second of (27) is the Gauss magnetic equation with a zero right side since no magnetic monopoles exist: the magnetic field strength (or magnetic field for short) H ¼ H e þ H n is solenoidal. The third relation in (27) is the Faraday law that couples electricity with magnetism and the fourth is Ampère's law that couples the total magnetic induction B with its sources j; j b .
The bound sources characterize the response of a polarizable/magnetizable system (not vacuum, a material) by means of two fields already defined in Table 5: polarization P and magnetization M. Mathematically, viz. [67,75] these are given by The polarization and magnetization fields are a consequence of the reorientation of the electric and magnetic dipoles. When present, they influence the mechanical characteristics of the originally isotropic active materials to become transversely isotropic or even orthotropic, see Fig. 18. In the figure, it is indicated that the isotropic plane (also called basal) will be defined by x 1 -x 2 and the direction of P, M is the perpendicular x 3 .
The electromagnetic constitutive equations are then defined as where D is the electric displacement or electric induction. In vacuum and for non-polarizable/non-magnetizable materials, P ¼ M 0 and the classical expressions D ¼ 0 E and B ¼ l 0 H are recovered. The phenomenological tensors v V , v u are defined as electric and magnetic susceptibilities. Introducing (28) into (27) and taking into account (29), the Maxwell equations for magnetizable/polarizable media may be rewritten as For quasi-static conditions and according to [51,64,120], a transformation through electric and magnetic susceptibilities can change (29) to x 3 x 2 x 3 x 2 x 1 x 3 Basal plane for transversely isotropic active materials and equilibrium interactions; polarization, magnetization vectors perpendicular to the plane where I is the identity matrix and , l are the tensor permittivity and permeability of the material. Both tensors depend on many factors: mechanical and thermal states, frequency in dispersive media, field strength in non-linear media and history in hysteretic materials. These dependencies influence the development of FEM tangent matrices as will be seen in Sect. 5. The expressions (31) will be expanded in Sect. 4.2.2 to include couplings. The Maxwell equations (30) are linear; non-linearities often emerge from the constitutive equations and/or the electromagnetic-mechanical interactions.
The equation of electric charge balance is obtained combining Ampère and electric Gauss laws from (30) and knowing that the gradient of the cross product of a vector is zero, to give

Compatibility Equations
As mentioned previously and according to classical field theory, viz. [120] for example, the electric and magnetic fields may be obtained either from scalar or vectorial potentials. If in the fourth (30) , the latter valid for EI and an approximation for NEI galvanomagnetic, then where A, u are the magnetic vector and scalar potential and V, V the electric vector and scalar potential, the last one commonly called voltage. In this paper u and V are adequate for the examples considered. This also renders the compatibility conditions (to be defined in Sect. 4.2.2) very simple and therefore the numerical implementation more direct, see Sect. 5. The presence of the inverse gradient r À1 warranties that the fundamental theorem of vector calculus (Helmholtz decomposition) is verified in the presence of free charges.
For the definition of the compatibility equations, it is interesting to note that for low frequencies _ B % 0 since the derivative is proportional to the frequency, and therefore also _ A % 0, see the last of (33). In any case in the current work mechanical dynamics are more relevant than electromagnetic ones: the latter are associated with the speed of light and the former with much lower speeds related with rigidity and mass.
It is important to emphasize that from the theoretical Physics point of view the vectorial choice is more correct but from the numerical one the use of vectors is obviously more expensive and also produces several problems that, to the best of our knowledge have not been correctly solved. For instance, the electric field is completely determined with a prescription of V in a single point of the domain, while that of V requires complicated gauge techniques that not always give proper results.

Poynting Theorems
The electromagnetic energy balance is given by the Poynting theorem (Fig. 19) [64].
Four main representations of the Poynting theorem exist depending on the choice of the Poynting vector j P (viz. - Table 7). This vector points (has the orientation) in the propagation direction of an electromagnetic wave and has dimensions of power/area. According to [75], the choice of the theorem is not relevant for non-dispersive (in theoretical physics terminology) or non-dissipative (in mechanics terminology) linear materials. However, for dispersive/ dissipative materials the choice of the flux vector j P may be decisive.
The Poynting theorem may be expressed in a general form by where c P , e V u , _ r P denote a scalar factor, the density of electromagnetic internal energy and a residual term, respectively. The residual term accounts for the power density and the dispersive terms, the minus sign indicates that the electromagnetic energy is transformed into another type of energy. Table 7 shows the expressions of c P , _ r P for the four theorem representations. For conductors or nonpolarizable/non-magnetizable materials the expressions of the residual are all equal and reduce to _ r P ¼ j Á E, which is known as the Joule effect.
In the Abraham representation, see next subsection, the residual _ r P depends on temporal derivatives of the material response (temporal description); in Minkowski's representation, it depends on spatial derivatives (spatial description). The other two are mixed representations and combine temporal and spatial descriptions.

Momentum Balance
The total momentum for an electromagnetic field interacting with matter is unique. However, the duality of the electromagnetic field in matter and in electromagnetic parts may be represented by two different equations due to the corpuscle-wave duality of light. During the last century physicists and mathematicians have debated about this duality. The two main theories were proposed in [1,94], causing the Minkowski-Abraham controversy (scientists in Fig. 20). In addition, several competing theories have emerged in the last years, viz. [20] for a full review. The Minkowski and Abraham theories are closely related to the choice of the Poynting vector for the representation of the linear momentum density, G Mi , G Ab respectively, shown in Table 7. Note that both densities are equal in vacuum but not inside matter.
The Minkowski momentum balance is obtained applying the chain rule to the momentum density G Mi defined in Table 7, introducing the third and fourth of (30) and using the vectorial equality ðr Â EÞ Â D ¼ r Á D E À ðD Á EÞ I ½ À E ðr Á DÞ þ ðr DÞ Á E and the vector counterpart for the magnetic field, to obtain with According to [42], T M is called the electro-magnetic stress tensor while in other books, for instance [33], is directly the Maxwell stress tensor. Note that the proper denomination should be Minkowski stress tensor since Maxwell only developed this tensor for vacuum media. The second equation defines the Lorentz electromagnetic forces. Taking into account (29), the equality ðr g À ð r EÞ Á P and the counterpart for the magnetic field, (36) becomes The Abraham momentum balance is obtained by introducing the vectorial relation 0 l 0 ðE Â HÞ ¼ D Â B À 0 l 0 ðE Â MÞ À P Â B into (35) and taking into account the Abraham flux from Table 7 yielding The underbraced last two terms on the right side are called Abraham or ponderomotive forces in [67], taken here with a negative sign. They are highly non-linear and can be considered second-order terms of the Lorentz forces. According to [58], the ponderomotive forces can be interpreted as the difference between the canonical (undulatory phase due to the electromagnetic field) and kinetic (corpuscular phase among particles) momentum densities.
Recently, [10] has concluded that G Mi represents the canonical momentum and G Ab the kinetic one. Therefore and again according to [58], G Mi describes the wave-like and G Ab the particle-like phenomena (interaction between  20 Minkowski (186420 Minkowski ( -1909, a German mathematician who worked on geometry of numbers and mathematical physics. Abraham (1875-1922) a German physicist who was a professor of rational mechanics at the Politecnico di Milano electromagnetism and matter). Indeed G Mi is a pseudo-momentum, conserved only in homogeneous and isotropic media, [67]; also, it is antisymmetric: the moment of momentum is not conserved. Both-momenta are equal in the vacuum. The present article studies the electromagnetic interactions within the matter as in [42]. In addition, the Abraham representation guarantees the conservation of energy and linear and angular momenta, viz. [89]. Therefore, the Abraham energymomentum representation is used in the remainder.

Multi-coupled Governing Equations
The aim of this section is to develop the multi-coupled governing equations used in this article from the expressions described in Sect. 3. As discussed in Sect. 1.1 distinctions between EI and NEI must be taken into consideration, otherwise, complex interactions such as viscosity or electromagnetic relaxations must be included in the formulation. Some of these issues are considered here but others are not necessary for the objectives of the present article.
The multi-coupled governing equations consist of balance equations, transport (also called constitutive) equations and boundary conditions. For a proper understanding, this section is structured as: Sect. 4.1 develops the general balance equations considering the mechanical, electric, magnetic and thermal fields. These balance equations are obtained using the two thermodynamic formalisms to develop a comprehensible formulation. From these multi-coupled balance equations, Sects. 4.2, 4.3 develop the multi-coupled governing equations for EI and NEI respectively.
At the end of the section and for summary, Fig. 24 shows a flow chart of simplifications and particularizations used to formulate the governing equations for EI and NEI.

Multi-coupled Balance Equations
The multi-coupled balance equations are composed of energy and momentum balances. In addition, the entropy balance is considered to obtain the transport equations using the procedure described in Sect. 3.1.2.

Energy Balance
Consider a thermodynamic universe consisting of a system and its surrounding as in Fig. 12. The total energy e of this system is the sum of mechanical e U , electromagnetic e V u and thermal e T energies, the three of them due to the interaction with the surroundings, Fig. 21. According to the energy balance (14), the total energy is unique and can be expressed by the following left expression Despite not being multiplied by q m , the electromagnetic energy is a proper density since the related fields E; D; H; B are already densities. The total energy flux j tt (right expression) is the sum of thermal q and of electromagnetic j P fluxes. A possible mechanical flux is not included due to (12), viz. [33] for further details. The multi-coupled energy balance is obtained by solving for q m _ e T in (39). Introducing (26) and the Abraham representation from inserting the first line of Table 7 in (34), yields after simplifications From the underbraces, it can be appreciated that this balance equation is the NEI expression of the first law (2). For EI, the non-equilibrium terms have to be eliminated and the time variation of (4) is recovered.
The balance (40) agrees with the one given in [124] although it does not satisfy the requirements of the relativity theory; an exact relativistic formulation is developed in [33]. Here, it is assumed that the velocity of the medium with respect to the observer is small compared with the velocity of light. In Sects. 4.2, 4.3, this balance will be specialized for each interaction.

Linear Momentum Balance
According to [33], the general conservation of total momentum takes into account the mechanical and Fig. 21 Total energy contained in a thermodynamic system: sum of mechanical, electromagnetic and thermal parts. External fluxes through boundary electromagnetic momenta without the internal volume forces f , f EM , f PM , and may be expressed by adding (25) and (38) to give Note that (41) does not include body forces and therefore is a conservation equation. The reference [29] justifies the validity of dropping these forces (sources/sinks) since the formulation of the electro-magneto-mechanic system is complete. For a FEM implementation, a form similar to that of the classical linear momentum (25) may be obtained. Thus, _ G Ab is replaced by solving for r Á T M in the Abraham momentum balance (38), inserting the result into (41) and substituting f EM from (37) to obtain If other forces, for instance gravitational ones from (25), are necessary they can be directly added to (42). Except for the first term, the other terms on the right hand side are of electromagnetic nature since temperature does not contribute to linear momentum. The fourth to sixth terms on the right side are defined as Lorentz forces in Sect. 3.3.4; the last two terms can be rewritten using the definition of the Maxwell stress (37) plus the tensorial definition in the text immediately before of the definition, to give With this transformation, the application of C 1 continuous finite elements is avoided using instead much cheaper C 0 continuous functions (viz. [153] for details). The tensor T M is non-symmetric due to the tensorial product terms, and represents the electromagnetic interaction with matter; in this particular form it is deduced from part of the Lorentz forces.
Due to the non-symmetry, the angular momentum conservation is not directly fulfilled. Then, the objective is to transform the sum T þ T M into a symmetric tensor, to automatically satisfy this conservation and to facilitate the FEM implementation, see Sect. 5. There is a controversy in the literature about the non-symmetry of T þ T M . For instance, [126] concludes that the electromagnetic body forces cannot be changed into a stress tensor since this mathematical transformation has no physical meaning. But [22,91] affirmed that it is valid, and using different approaches established that T is non-symmetric, but the sum T þ T M is symmetric. In the present work, a solution similar to the one given in [91] is proposed: we split the Cauchy-like stress tensor into two parts: T ¼ T C À T M sk , where T C is the classic and symmetric Cauchy tensor. As any other tensor, T M can be decomposed into symmetric and antisymmetric parts, The terms inside the parenthesis are the vacuum contribution to the tensor, diagonal and therefore always symmetric.
With the previous split, s is now symmetric.

Entropy Balance
In (2), it is clear that both equilibrium and non-equilibrium terms are combined. Therefore, the general, multi-coupled entropy balance will require the use of the two formalisms ET, NET. Starting with the ET formalism and considering (11), the expression (4) for EI can be expressed as a function of de and ds. Differentiating (4) with respect to time, the multicoupled entropy becomes For the NET formalism, substituting (7) into (1) gives The first two terms of the right hand side would give an equation equal to (45) but exchanging T e by T n in the denominator due to the difference between (3) and (7) In this equation, the first quotient comes from d _ E T and the second from d _ W; the contribution of d _ W can be added substituting q m _ e T by its complete (from the thermodynamics point of view) expression from (40), giving Note that (40) was deduced from a complete multi-coupled continuum physics balance, that is, it includes both ET and NET contributions. The following definitions already introduced in Sects. 3.2.1 and 3.3.1 have been used • NET stress tensor: The entropy flux j r and entropy production r terms are deduced by comparing (47) and (16) to give It can be appreciated that j r in the current formulation is proportional to the heat flux. The first two terms on the right side of the bottom equation are related to thermal conduction and to Joule heating, respectively. The contribution of the last three irreversibilities to the entropy production is due to viscous interactions [33], dielectric relaxations [34] (discovered by Debye, viz. Fig. 22) and magnetic relaxations [122].
It is interesting to formulate these NEI since in some applications they are useful, however, for simplicity in the present work they are not included; therefore Simplification 1: Transport interactions due to NET stress, electric and magnetic fields are not considered.
Due to this simplification, the last three terms in the fraction of the second equation (48) are dropped. Since _ S, _ P, _ M are extensive variables present both in EI and NEI in general they will be non-zero; then, from (47), T n ¼ E n ¼ H n ¼ 0 and therefore T e T, E e E, H e H. In the remainder of the paper there will be no distinction between total and equilibrium intensive variables except for the temperature.
An important corollary of S1 is that the governing equations will be uncoupled, permitting us to develop formulations for EI and NEI separately. Thus, the multicoupled entropy balance, entropy flux and entropy production from (47) are given by The irreversibilities are here reduced to those of heat flux and Joule heating. An additional simplification is stated for both formalisms Simplification 2: No dynamic free electric charge densities are present,ρ f q = 0.
S2 implies _ D % 0 from (32) or from the Gauss law first of (30): exact for EI and a good approximation for most NEI applications in which relativistic phenomena are not present, according to [80]. In fact, a simple calculation shows (viz. [113]) that for frequencies lower than 10 5 [Hz] the conduction electric flux densities are five orders of magnitude higher than those due to _ D.

Energy Balance Modified: Equilibrium and Nonequilibrium Temperatures
Due to S1, mechanical, polarization and magnetization irreversibilities disappear from the formulation. Consequently, the energy balance (40) must be modified to eliminate the non-equilibrium terms so that only the EI variables T e , E e , H e that reversibly exchange energy remain.
The simplified multi-coupled energy balance relation is obtained from (4) and establishes that S must be an exact differential in S, E e , H e , T e since, due to the lack of irreversibilities in (1). To this end, the last in the matrix equations (6) will be transformed in Sect. 4.2.2 to give its local form in the fourth of (55). Introducing the time derivative of this equation into (46), the result into (40) and assuming that the conduction temperature is similar to the reference temperature, T n % T 0 (in Kelvin degrees), only for the thermoelastic interaction (see for instance [15,132]), the multi-coupled energy balance is where c is the heat capacity and the first-order tensors p V are the pyroelectric and p u the pyromagnetic properties, all of them to be found experimentally.
[To be exact, in (50) _ E e , _ H e should have been used but as established by S1 they are equal to the total fields.] The T 0 Biot damping terms include the thermomechanical, pyromagnetic and pyroelectric converse electromagnetic dynamic interactions as shown by the triangles numbered 7, 6, 5 in (6). The term b: _ S often is called two- way thermoelasticity and is due to the heat generated by the mechanical vibration. It is important in the design of some MEMS and NEMS devices [132]. Note that the Biot terms are not of the same nature as the transport interactions from (48) that are neglected by S1. The former ''interchanges'' energy among the four fields in a reversible manner through balance equations while the latter ''dissipates'' in an irreversible way through constitutive equations.
The first two terms on the right hand side of the first (50) represent non-equilibrium interactions (heat conduction and Joule effect, respectively) and are studied using the conductive NEI temperature T n . The related gradient rT n is a driving force, as listed in Table 6, with the heat flux as conjugate. On the other hand, the damping terms include the equilibrium temperature T e since they represent reversible interactions or intensive variables closely related to entropy.
A ''two temperature'' model is sometimes considered to study the fully-coupled thermoelastic coupling, viz. [28,40,53,148,149]. This model relates both temperatures using a constraint d [ 0 and the additional PDE The parameter d is adjusted towards experimental data, and the fictitious heat flux Q d rT n is defined for the twotemperature model. Due to the inherent difficulty of a numerical implementation of the obtained multi-coupled balance equations: linear momentum (43), entropy (49) and energy (50), these balances will be specialized in the following section for EI and NEI, taking into account one new simplification.

Multi-coupled Governing Equations for Equilibrium Interactions
This section presents the multi-coupled governing equations (balance plus boundary conditions) for EI, obtained introducing a new simplification.
S3 implies that the media to be studied by EI are polarizable and/or magnetizable, see Sect. 1.

Balance Equations for EI
The balance equations for EI are: momentum balance (43) with symmetric stress tensor, Gauss electric and Gauss magnetic laws (30), two-temperature coupling (51) and energy balance (50). Applying S3 to the first and last equations, the complete balance is In order to use a matrix formulation, the stresses T C ; T M s in some instances will be given in Voigt notation; the threedimensional mechanical (without couplings) stress-strain relation in a 6 Â 6 matrix form is expressed as T C ¼ C : S and in principal directions is given by The shear strains are defined as engineering ones. With this definition, the thermal expansion second-order tensor b introduced in (50) is calculated as where a T are the thermal expansion coefficients, equal in the directions of the basal plane but different in the perpendicular direction: a 1T ¼ a 2T 6 ¼ a 3T (see Fig. 18).

Compatibility and Constitutive Equations for EI
According to [134], eight sets of multi-coupled constitutive equations can be defined, depending on the choice of the independent state variables. In this work for both EI, NEI, from (23) and the two remaining Maxwell equations transformed in (33), the independent variables are taken as and, thus, permit a displacement-based FEM formulation. The multi-coupled constitutive equations are obtained from the electromagnetic enthalpy P, [62]. Considering • the Gibbs expression (4) • that P can be expressed as an exact differential of the intensive equilibrium variables [ Table 5, (5)] the multicoupled Clausius-Duhem inequality can be obtained; the second item follows from the fact that for EI there are no dissipations. For an FEM amenable implementation, a Legendre transformation to exchange E by P, H by M and T e by s is applied to (4): to the total thermodynamic energy e T , other energies produced by the exchanged variables must be subtracted to obtain the enthalpy result P ¼ e T À E e Á PÀ l 0 H e Á M À T e s. A similar procedure, although with different terms, is developed in [16,45,88,134,141]. In this work, we obtain where e V , e u , m are the piezoelectric, piezomagnetic and magnetoelectric couplings. Following the six-term expression of the stress tensor given in (53), the matrices e V , e u are 3 Â 6 and m is a diagonal 3 Â 3 matrix. These are shown in Appendix 4. For more information on these properties and their measurement techniques, see the piezoelectric [61] and magnetostrictive [62] standards and the Appendices. In essence, the procedure consist on assuming material linearity and a natural state (initial S ¼ E ¼ H ¼ 0; T e ¼ T 0 ) the enthalpy can be calculated as the first and second terms of Taylor series in the neighborhood of the natural state. If reversible processes are considered, all variables are exact differentials as in (5) and using the thermodynamic Maxwell relations, second derivatives of the enthalpy produce the material property tensor. For example, o 2 P=oSoE e ¼ Àe V .
Polarization P R , magnetization M R and stress T R residuals are introduced, the latter can arise for instance, from manufacturing or precompression. P R , M R often are present in insulating materials: ferroelectrics for the first due to reorientation of dipoles and ferromagnetics for the second due to magnetic memory. In this work, they are assumed to be independent on space or on variations of the corresponding field. The presence of these residuals in the formulation could permit future developments, for instance on domain switching or memory effects.
Deriving the previous equation with respect to the independent variables and including the thermal flux from (51) where the equalities from (31) have been used and where j is the thermal conductivity. Due to the two-temperature model, fifth in (55), the standard Fourier law also has to be included to calculate T e . This temperature could be inferred from the entropy, fourth equation, forcing a new nodal variable in the FEM model. Instead T n is calculated from Fourier and T e from the two-temperature model. We note that with this formulation the implementation of mixed finite element models is avoided.

Boundary Conditions for EI
For multi-coupled interactions the Dirichlet C u and Neumann C p boundaries shown in Fig. 16 are split into mechanical C U u , C U p , electrical C V u , C V p , magnetical C u u , C u p and two different thermal C T e ;T n u , C T e ;T n p parts, see Fig. 23. These boundaries must satisfy  Table 8 lists the possible boundary conditions: u, V, u are prescribed displacements, voltage, scalar magnetic potential and T e , T n temperatures for EI, NEI, respectively. Furthermore, t C denotes prescribed Cauchy tractions, q, B prescribed thermal and magnetic fluxes and D surface charge density, the boundary part of q f q . The prescribed magnetic flux has no physical sense but it is incorporated to define boundary magnetic fields in numerical computations. Note that there is no flux related to the equilibrium temperature T e since (51) is only a numerical constraint without physical meaning.

Multi-coupled Governing Equations for Nonequilibrium Interactions
Prior to specifying the multi-coupled governing equations for NEI, a new simplification S4 is introduced.
This simplification is a good approximation for conductor or semi-conductor materials studied in this paper.

Balance Equations for NEI
Considering the applicable simplifications, the balance equations for NEI are the linear momentum balance (43) without ponderomotrive forces (they depend on P; M), electric charge balance (32) with S2, the two-temperature model (51) and energy balance (50). In the last equation, the pyroelectric and pyromagnetic p V , p u effects associated to equilibrium fields are eliminated. The magnetic Gauss law from (30) is also included to obtain a scalar potential formulation.
S4 and (29) imply D ¼ 0 E, B ¼ l 0 H; then for NEI the Maxwell tensor in vacuum from (37) or (44) becomes Its magnitude is in general very small due to the small values of 0 , l 0 but again its inclusion can be relevant in some applications. We note that in NEI both stress tensors are symmetric without need of further approximations.
In this work NEI is applied to isotropic materials; then the stress-strain relations in indicial form are given by where k, l are the Lamé parameters and d ij the Kronecker delta. Also, the equal entries of the thermal expansion tensor b can be directly calculated from a T ð3k þ 2lÞ.

Transport Equations for NEI
According to the NET formalism [69], the transport equations are obtained by expressing the entropy production r from the second of (48) and subsequently introducing S1. Considering the form of (18) and (19) this yields where the four second-order tensors L kl ðBÞ include the material properties that satisfy the thermodynamic restrictions given in (20). From [80,105], the two matrices that group both types of tensors are here arranged and modified to Also, the material properties are converted into tensor entities a, q, j through Fig. 23 Essential, natural boundary conditions on boundaries C u , C p respectively for equilibrium interactions with mechanic ðÁÞ U , electric ðÁÞ V , magnetic ðÁÞ u and thermal ðÁÞ T e ;T n fields where a, q are defined as Seebeck and electric resistivity. Also, N, R, M, are the Nernst, Hall and Righi-Leduc coefficients, related to the interactions defined in Sect. 1.1.3. To operate the cross product in vectorial notation, the subindex Â indicates the reallocation of a 3 Â 1 vector into a 3 Â 3 antisymmetric matrix with zero diagonals; for instance for the magnetic induction in (60) gives Due to the dependency of the submatrix L kl on magnetic induction, for experimental measurements and to reduce the related algebra it is convenient to split (59) into symmetric and skew-symmetric parts, L kl ¼ L kl s þ L kl sk . The first part (transport properties) describe the classical phenomenological properties in the absence of B and the second the galvanomagnetic and thermomagnetic effects.
As indicated in (60) the a; q; j depend on the conduction temperature, as will be quantified in Sect. 7; in what follows the non-equilibrium argument of T n will be omitted to simplify the relations. The antisymmetry of the second parts in (60) comes from the form of B Â . The explicit matrix expressions for the three relations are The above definitions satisfy the reciprocal relation (21) or equivalently j sk ðBÞ ¼ j sk ðÀBÞ > etc., assuring that the entropy always increases, as expressed in (20). To be used in the transport equations, the tensor of electric conductivities are defined as c À1 ¼ q. Due to the inversion process c cannot be divided into symmetric and skew-symmetric parts, although the particular case of the absence of magnetic field directly gives c ¼ q À1 s ¼ I=q.
Each term of (59) or equivalently from (60) clearly represents a physical effect. In the absence of B and from (58) after dropping the tensorial notation • At the diagonals, rV=q, j rT n , the heat and electric laws stated by Ohm and Fourier. • At the top off-diagonal, a T n j, the Peltier and Thomson effects described in Sect. 1. • At the bottom off-diagonal (along with the second product in the low diagonal), a rT n =q, the Seebeck effect.
and in presence of B as defined in Sect. 1.1.3 • From the skew term on the right side of a in (60) the Ettingshausen and Nernst effects, coupled and both controlled by N. • From the skew term on the right side of q in (60), the Hall effect. • From the skew term on the right side of j in (60), the Righi-Leduc effect.
To sum up, considering the transport equations (58), the mechanical and magnetic constitutive (55) and the twotemperature model (51), the transport (third, sixth) and constitutive (the rest) equations for NEI are where the permeability for conductor or semiconductor materials is assumed to be the vacuum l 0 , since magnetization interactions are not present. The electric and thermal transport equations, third and sixth in (63), are now interpreted from a physical point of view • For the electric, the first term on the right side represents electric conduction, the second Seebeck effect • For the thermal, the first term represents heat conduction; the second Peltier effect

Boundary Conditions for NEI
The Dirichlet and Neumann boundary conditions are shown in Table 9. The only magnitude different from those of Table 8 is j, the prescribed electric flux.

Finite Element Formulation
A variational formulation within the FEM is developed in the present section geared to implement in a computer program the multi-coupled governing equations for: (1) EI, using (52), (55) plus Table 8, (2) NEI, using (56), (63) plus Table 9 (see Fig. 24). Both finite elements (FE) are implemented in the research program FEAP [138] from the University of California at Berkeley; this program provides several dummy routines (user elements) to implement new modular elements using programming languages Fortran or C.
The standard compatibility equations are listed in (54). The use of a magnetic vector potential A would modify the magnetic compatibility to B ¼ r Â A, see the first (33). In addition, the fourth Maxwell equation (30) would be Ampère's law r Â H ¼ j þ _ D. But the magnetic scalar potential u is a good approximation for domains without electric fluxes. This situation applies to EI for which S3, and to NEI for which S4 hold, as explained in Sects. 4.2, 4.3. However, for a complete and fully multicoupled formulation the use of the vector potential A should be considered.

Finite Element Discretization
A standard isoparametric formulation is used in which the element shape functions N are expressed in terms of parametric coordinates n ðn; g; fÞ. At any point of the parametric coordinates in the element domain the cartesian coordinates and the dependent variables along with their time derivatives are discretized (approximated) with standard shape functions as (viz. [153]) where a repeated subscript a implies summation over nodes numbered 1; . . .; nel of a particular finite element and the supraindex ðÁÞ h denotes FEM approximation. For the mechanical dof, the three displacements at any point in the element domain expressed in parametric coordinates are uðnÞ ¼ fu 1 ðnÞ; u 2 ðnÞ; u 3 ðnÞg > .
The vector a U a ¼ fa U 1a ; a U 2a ; a U 3a g > groups the displacements at a local node a with coordinates n a ; N a N a ðnÞ is the value of the anode related isoparametric shape function at the generic location n.
The gradients are also discretized with the same shape functions where the FEM gradient matrices for mechanical, vectorial (electric, thermal or magnetic) and galvano ''compatibility'' are expressed as Table 9 Summary of Dirichlet (essential), Neumann (natural) boundary conditions for nonequilibrium interactions Dirichlet Neumann

Eq. (55)
Eq. (63) In the previous equations N a;j denotes the partial derivative of N a with respect to the coordinate x j . For the definition of these matrices, the small displacement strain tensor definition after (23) and the shear strain order given in (53) have been used. Similarly to N a , B s a B s a ðnÞ, B a B a ðnÞ are the matrix values of the a-node related derivatives of the shape function at the generic location n.
Finally, the variations and their gradients are also discretized as

Finite Element Formulation for Equilibrium Interactions
Consider the system shown in Fig. 16; in the related governing equations described above many practical material or constitutive non-linearities are contemplated. The material linearity can be due to that of most constitutive equations and the geometric one to the small displacement assumption. The non-linearities of the model emerge from the Maxwell stress tensor and from some material properties that depend on the field intensity. From the numerical point of view non-linearities are addressed by the Newton-Raphson algorithm and time integration by standard Newmark schemes [50,100].

Weak Forms for EI
As described in [153], see Fig. 25, while the original partial differential equations (also named strong forms) include second derivatives, the weak forms (or integral equations) only involve first spatial derivatives of the dof's.
According to standard variational methods reported in the previous reference, the weak forms are obtained multiplying each balance equation (52) by variations of the appropriate dof's du, dV, du, dT e , dT n and integrating over the domain X. The final weak forms (68) are obtained by applying the divergence theorem to the gradients of all equations and enforcing the Neumann boundary conditions from Table 8. This yields After this equation and following the traditional FEM, the tensor notation will be dropped in this section and substituted by vector notation including the transpose symbol.

Discretized Constitutive Equations for EI
At any generic point n of a particular finite element, the stresses, fluxes and fields from (44) and (55)   where the repeated subindices ðÁÞ b play the same role as that of ðÁÞ a in (64).
In what follows, both stresses T C , T M s will be expressed either in a 6 Â 1 vectorial form or in a 3 Â 3 tensor form, as required for the context. For clarity, the argument n will be dropped.

Discretized Residuals for EI
The Galerkin approach permits us to convert the continuous formulation given in (68) into an amenable (for numerical analysis) discrete formulation using the approximations in (64) to (67).
Following the standard FEM method [153], the integral domains are discretized into elements: R X dX P e R X e dX, the residuals are computed element by element and assembled. The resulting residuals are expressed as The sign has been changed in the second to fourth equations to obtain a better conditioned multicoupled matrix in Sects. 5.2.4, 5.2.5. The previous boundary integrals are nonzero only if the element has one or more sides with prescribed magnitudes in nodes that lie on the boundaries C p . In the residuals the following forces and fluxes must also be discretized. From (29), (38), (50) where the antisymmetric P Â , M Â have the same profile and signs as B Â in (61). The time derivatives _ B Â , _ M Â can be directly calculated.

Tangent Stiffness Matrices for EI
The tangent ''stiffness'' matrices (using a Continuum Mechanics terminology) are obtained by linearizing the residuals (70) with respect to the dof. The non-zero contributions from these residuals are listed in this subsection, and in general, these submatrices are non-symmetric. For each pair of local nodes a; b of the element and for the dof's I; J U, V, u, T e , T n .
The mechanical residual gives the direct and coupling submatrices There is no coupling with the conductive (non equilibrium) temperature since it is only present in the heat fluxes and not coupled with the stress, electric or magnetic fields. The electric residual couples with itself and also with the other fields: this residual only contains the variable induction D that in turn is coupled with mechanical, electric, magnetic and equilibrium thermal dof For the same reasons, the magnetic residual gives The diagonal submatrices K VV ab , K uu ab will be positive since the corresponding derivatives of D, B are negative.
Finally, the two thermal residuals produce the non-zero direct and coupling submatrices Note that the first stiffness shows the unusual characteristic of not including a material property; this is due to the nature of (51), as mentioned it is a constraint rather than a constitutive equation.

Tangent Capacity Matrices for EI
The tangent ''capacity'' matrices (with the terminology of Heat Transfer) are similar to the tangent stiffness but the derivatives are with respect to the first time derivative of the dof Capacity terms arise from the ponderomotive forces of the mechanical residual and are expressed by The general non-zero terms of the non equilibrium temperature residual from the Biot and standard heat conduction terms give

Tangent Mass Matrices for EI
In the tangent ''mass'' matrices (again with the Continuum Mechanics terminology) the derivatives are with respect to the second time variation of the dof.
The only second derivative with respect to time is related to the mechanical dof, that is For simplicity and from the several choices for mass matrices, we use a diagonal-lumped version. The three types of matrices: stiffness, capacity and mass, include partial derivatives that are developed and listed in Appendix 1.

Assembled System for EI
The previous matrices and residuals are assembled element by element to give a linearized algebraic system, with the nodal 7Â7 matrix form of (79). In this matrix, the parameters c 1 , c 2 and c 3 arise from the time integration formula used. For the Newmark method 2 where c N and b N are user-chosen parameters, viz. [153]. The multi-coupled assembled matrix for EI is nonsymmetric, consequence of the lack of reciprocity of some of the couplings. From a numerical point of view, this asymmetry is a drawback for the efficiency of the system solving; special asymmetric equation solvers must be used to guarantee quadratic rates of convergence. This difficulty is avoided using the FEAP command UTAN, viz. [138].
At each non-linear step the nodal unknowns are updated from the previous ones calculated in (79) with Finally, all integrals from Sects. 5.2.3 to 5.2.6 are calculated numerically at Gauss points.

Finite Element Formulation for Nonequilibrium Interactions
Consider again the system shown in Fig. 16; the NEI of this system are fully defined by the balance (56) and constitutive (63) equations and corresponding boundary conditions from Table 9.
To be used in the FEM residuals, (63) are discretized using (64), (65). The discretized compatibility expressions (69) will be used again. All shape functions, material properties etc. in the right hand sides of the following equations will be evaluated at n but the argument is omitted for simplicity. The nodal values are directly related with the corresponding node b.

Discretized Residuals for NEI
As done for EI in the previous section, the weak forms are obtained multiplying the balances (56) by the variations from (67) and integrating over the FE domain. Applying the divergence theorem to the first term of all equations, the Neumann boundary conditions from Table 9 and the FEM approximations from (64), (65) the residuals are Again the signs of the second, third and fourth residuals have been changed. The mechanical residual R U a is a 3 Â 1 vector entity since includes three mechanical dof's or displacements. Consequently, five residuals and seven dof are present. Both stress tensors must have a 6 Â 1 dimension.

Tangent Stiffness Matrices for NEI
As explained in Sect. 5.2.4, the tangent ''stiffness'' matrices are obtained deriving the residuals with respect to the dof's. There will be contributions from the five residuals (80), although several of the derivatives will be zero. The mechanical residual produces the non-zero direct and coupling submatrices Besides the principal K UU ab , the mechanical field interacts with the electric and magnetic fields and the thermal (both equilibrium and non-equilibrium) through the Lorentz forces (terms with j Â B). In addition, mechanical and thermal fields are also coupled since the stress depends on temperature, first of (63). The derivative of j Â with respect to a U b (first equation) is zero since piezomagnetic effects are not included.
The electric and magnetic residuals produce the following non-zero submatrices In the previous expression, K V u ab represents the galvanomagnetic interaction due to the magnetic dependency of the electric flux developed in (63).
For the magnetic dof, the only non-zero term is the principal interaction, third in (82). The nullity K u T n ab ¼ 0 is due, for conductor materials, to the absence of pyromagnetic interactions. In addition K u V ab ¼ 0: the electric flux alters the magnetic field according to the fourth Maxwell law from (30), however, the choice of a scalar magnetic potential implies that B only depends on u. A similar approximation was also assumed in [104,105]. The thermoelectric K VT n ab interactions are observed again.
The derivation of the fourth residual R T e a directly gives the first two submatrices from (75); those of the fifth residual are The thermal and electric fields are coupled by several thermoelectric interactions in K T n V ab that are described in detail in Sect. 4. Also, K T n u ab takes into account the thermomagnetic interactions given by a, j in (60).

Tangent Capacity Matrices for NEI
The ''capacity'' matrices represent the dissipative interactions, and are obtained deriving the residuals with respect to the time first-derivative of the dof. The following submatrices are non-zero The first expression represents the thermoelastic damping that was introduced in [15] to study the thermomechanical interactions, and arises from the thermal fluxes that have been originated from the volumetric strain variations, Biot's damping term. The heat transient is represented by the two-temperature matrix C T n T e ab .

Tangent Mass Matrices for NEI
The only non-zero ''mass'' matrices are related with second derivatives and represent mechanical inertias. With the reallocation of the shape functions described for EI, we obtain the lumped matrix

Assembled Matrix for NEI
As mentioned in Sect. 5.2.7, the stiffness, capacity and mass sub matrices are grouped for each node in the 7 Â 7 matrix (86) are implemented in the research code FEAP.

Numerical Examples for Equilibrium Interactions
The basic validations for EI are divided into linear and nonlinear interactions (in the sense of cause-effect responses, see Fig. 5). The linear interactions are closely related with the constitutive equation (55) and the non-linear ones with the Maxwell stress tensor from the same equations.
Although the complete study of EI often involves the analyses of complicated geometries (see for instance [112]) in this work and for demonstration purposes a simple parallelepiped will be used, with a material that couples elastic electro, magneto, and thermal fields. Table 10 lists the validations selected that were already introduced in Tables 1 and 2. Cases named electric and magnetic susceptibilities, piezoelectric, piezomagnetic, magneto electric, pyroelectric, pyromagnetic are linear interactions; the other electrostriction and magnetostriction are non-linear. Other interactions such as elasticity, heat capacity and thermal expansion are not completely present since they are basic and already have been validated in ongoing work from [117,118] for NEI.
For all cases, the 2-D geometry shown in Fig. 26 is considered, with l 3 ¼ 1:14 Â 10 À3 [m] and l 2 ¼ 3l 3 . For simplicity and to avoid rigid body motions, all Dirichlet boundary conditions are set to zero at the vertical and horizontal bottom sides. On the top side, the prescribed Dirichlet electric, magnetic and thermal boundary conditions for each case are given in Table 10. Neumann conditions are automatically zero in all the boundaries.
Since the voltage, magnetic and temperature scalar unknowns are potentials, an analytical solution for the potential distribution can be calculated, viz. [130] in which the Laplace transformation technique is applied. Denoting Ç to a generic scalar potential that in this context can be T e , V, or u where @ ¼ pð2k À 1Þ=l 2 , and Ç is the particular scalar field prescribed at the top face. According to the compatibility Table 10 Cases and prescribed boundary conditions at the top face from Sect. 4.2.2 and the third of the constitutive (55), the three components of the field associated to each Ç are N ¼ ÀrÇ , therefore Figure 27 shows the potential distribution given by (87), assuming the prescribed potential to Ç ¼ 10. This solution satisfies well the natural boundary conditions with k ¼ 90 terms, except close to the prescribed surfaces where small oscillations appear.

Case I
The electric susceptibility interaction relates electric field (cause) and polarization (effect) through the medium permittivity, see the second row of Table 1.
From a numerical point of view and for this validation, only the voltage dof is required. Therefore, the constitutive equation of the electric displacement (fourth in (55)) reduces to D ¼ Á E and the relationship between electric field and polarization (29) top (in absence of residual polarization, as in the following cases) becomes P ¼ ð À 0 IÞ Á E. The polarization along the x 3 direction in Fig. 26 is therefore given by Figure 28 shows the component P 3 along the direction x 2 at x 3 ¼ l 3 =2 as shown in Fig. 26; this line is used for the rest of EI cases. The polarization is negative due to the sign of the electric field, from larger to smaller voltages and follows a non-linear distribution according to (88). The values of 0 , , as will be done for others in the following cases, are taken from Appendix 4. As in the remaining EI cases, analytical and numerical results are very similar, except close to the sides-a finer mesh in these areas would improve the solution.
The electric field E 3 (and the proportional polarization P 3 ) increases very quickly along x 2 close to the vertical sides due to the closeness of the equipotential lines of Fig. 27 in this zone. In the center, these lines are almost parallel, therefore P 3 is constant and maximum.

Case II
The magnetic susceptibility interaction relates magnetic induction (cause) and magnetization (effect) through the permeability of the medium. For this validation, only the scalar magnetic potential is required. The constitutive equation of the magnetic induction (fifth in (55)) is reduced to B ¼ l Á H and the relationship between magnetic induction and magnetization (in absence of residual magnetization, as for the rest of EI cases) is from (29) Therefore, the vertical magnetization (Fig. 26) is given by where Ç u, N 3 H 3 for this case. Figure 29 shows the component M 3 , and as in the previous case, M 3 is nonlinear and negative.

Case III
The magnetoelectric interaction (direct or converse) relates magnetic and electric fields through the magnetoelectric coefficients m. In the present validation, the direct interaction is considered, prescribing a magnetic field (cause) and calculating the electric displacement (effect). Again from (55) the reduced constitutive relation is D ¼ m Á H, therefore D 3 ¼ m 33 N 3 . In this equality it is implicit that Ç u, N 3 H 3 . The component D 3 is again non-linear and negative since D has the same sign as H before and this magnetic field is negative. For this validation, two dof's-voltage and magnetic scalar potential-are considered. Figure 30 shows the numerical distributions of scalar magnetic potential (top) and coupled induced voltage (bottom) for the magnetoelectric interaction.
As expected, the magnetic potential has the same distribution as that of Fig. 27. The maximum gradient (maximum magnetic field) clearly concentrates at the upper corners, and due to the coupling, a strong voltage gradient (Fig. 30 bottom). This voltage distribution is very similar to the one generated by two point electric charges applied at two internal points close to the corners.

Case IV
In this case, the pyroelectric interaction relates the thermal and electric fields through the pyroelectric coefficients p V ; again, two interactions (direct and converse) are possible. Assuming that the initial temperature is T 0 ¼ 0, the constitutive equation for the electric displacement from (55) For this validation, two dof's (equilibrium temperature and voltage) are required and the analytical solution is (87) for coordinates x 2 ; x 3 corresponding to the top face. Figure 31 shows the electric induction D 3 (identical to , distribution similar to that of Fig. 27 for T n along a horizontal line, but scaled by the pyroelectric coefficient p V 3 and with sign changed. The positive sign is due to the that of T e and the related coefficients p.

Case V
The pyromagnetic interaction relates thermal and magnetic fields through pyromagnetic coefficients, with the two possible interactions: direct and converse. The constitutive equation for the magnetic induction last of (55) reduces to B ¼ p u T e again with T 0 ¼ 0.
For this validation, two dof (temperature and magnetic scalar potential) are required and the analytical solution is where Ç T n in the top face as before. Figure 32 shows the component B 3 , with the same distribution as in Fig. 31 but exchanging p V 3 by p u 3 and consequently D 3 by B 3 .

Case VI
The piezoelectric interaction relates electric and elastic fields; depending on the field considered cause, the piezoelectric interactions can be direct or converse. This   (Fig. 26). Line analytical, circles finite element results case validates the direct piezoelectric interaction for which the electric field is the cause and mechanical stress the effect. Four dof's are required in this validation: three mechanic displacements and voltage. Therefore, the stress constitutive, first of (55), reduces to T C ¼ Àe V Á E. The Cauchy stress field is then given by The applied electric field generates non-linear stress distributions along directions out-of-plane x 1 and vertical x 3 , Fig. 33 shows T C 11 (solid line) and T C

33
(dashed line). The different signs and magnitudes of the two stresses are due to those of e V 13 and e V 33 , see Appendix 4. These distributions are the common ones for piezoelectric materials polarized in the x 3 direction.

Case VII
Electrostriction is a non-linear, electro-mechanical interaction modeled by the superposition of two terms s : a linear one given by the constitutive (without residual stresses) (55) and a non-linear represented by the symmetric part of the Maxwell stress tensor (44), giving The simulation of this interaction requires four dof's: three displacements and voltage. The objective is to investigate the contribution of the Maxwell stress tensor to the total one. Therefore, from both (89), neglecting the piezoelectric interaction and the Cauchy stress where Ç V, and N 3 E 3 . Figure 34 shows the component T M s33 . This validation is completely non-linear due to the quadratic term of the Maxwell tensor in (90). Three iterations of the Newton-Raphson algorithm are necessary to converge to the solution. From a physics point of view, these iterations are corrections of the Cauchy tensor shown in Fig. 33.
Although in this case the correction is five orders of magnitude smaller than the main stress, the Maxwell component depends quadratically on the applied electric field E 3 and on the permittivity 33 . In some applications, T M s33 can be of the same order or even larger than the Cauchy stress. Note in Fig. 34 that the slope of the curve tends to zero at the lateral sides, a typical characteristic of the estrictive behavior.

Case VIII
The piezomagnetic interaction relates magnetic and elastic fields, under direct or converse interactions. This validation is similar to Case VI replacing electric by magnetic field, or alternatively electric voltage by magnetic potential. The solution for the direct interaction is therefore given by where Ç u, and N 3 H 3 . Figure 35 shows the non-linear, out-of-plane Cauchy stress T C 11 (solid line) and vertical T C 33 (dashed line). Due to the equal sign and similar magnitude of the piezomagnetic constants, both stress components are similar. As in Fig. 33, the values at the vertical sides are small but not zero, due to the small gradient of H 3 at these sides. It can be appreciated that the magnitude of the stresses are not negligible, considering the low magnetic potential prescribed and that active materials often are brittle and/or of low strength.

Case IX
Magnetostriction is a non-linear, magneto-mechanical interaction that as described in Case VII superposes T ¼ T C þ T M s , and is studied reducing the constitutive equation to From a numerical point of view, the simulation of this interaction requires four dof's: three displacements and magnetic scalar potential.
From both (91), neglecting the piezomagnetic interaction and the Cauchy stress, the contribution of the Maxwell stress tensor to the total stress is given by where Ç u, and N 3 H 3 . Figure 36 shows the stress component. Again, this interaction is non-linear and the Newton-Raphson algorithm is used for resolution. The Maxwell stress correction is more significant than in the electrostrictive Case VI, although still three orders of magnitude smaller than the Cauchy main contribution. This fact can be the reason for which the magnetostrictive interaction is almost always considered in the study of practical devices. Other interesting result from this validation is that the non-linear stress is always negative, independently of the sign of the applied magnetic field.

Numerical Examples for Non-equilibrium Interactions
In this section, a set of cases to validate most of the NEI formulation described in Sect. 1 is presented. For this purpose, a sample of indium antimonide (InSb), a metalmetalloid alloy showing a thermoelectric behavior, is considered. The InSb thermoelectric properties depend on temperature and are taken from [105] but to facilitate the obtention of validating direct solutions, constant values at an approximated average temperature ð T n h þ T n c Þ=2 are used, see Table 12.
The Hall, Righi-Leduc and Nernst coefficients are also dependent on temperature, viz. [105]. Again for simplicity the non-linearity of these magnetic properties in the FEM formulation is not considered, In Table 11 the cases and corresponding boundary conditions are listed. Symbols with an overbar represent the prescribed quantities that define the case and the subindices refer to the hot and cold faces and to Cartesian directions. A p-type semiconductor is considered with the corresponding positiveness of the Seebeck and Nernst coefficients.
The geometrical dimensions of the InSb sample are depicted in Fig. 37: l 1 ¼ l 3 ¼ 1:14 Â 10 À3 [m], l 2 ¼ 3l 1 , the last two as in Sect. 6. These dimensions are typical of commercial devices [93], except for l 2 being multiplied by a factor of 3 to partially avoid boundary phenomena produced by the magnetic fields.
In total, there are seven cases related to NEI. The properties of InSb show in Table 12 are taken from [63,105]. To isolate each effect, only the relevant coupling properties are set different to zero.

Cases I and II
Case I is a first-order Fourier effect, a non-equilibrium interaction that relates temperature gradient (cause) and heat flux (effect) through the thermal conductivity j firstorder property, see Fig. 8 and Table 3. Mathematically, the classical heat conduction model is obtained considering a ¼ 0 in the last (63), remaining only the NEI temperature dof. Figure 38 left shows this distribution, linear due to the absence of heat sources. The corresponding constant heat flux can be calculated analytically with q 3 ¼ Àj DT n =l 3 ¼ À2:37 Â 10 5 [J/sÁm 2 ], value very similar to the constant one obtained with the FEM everywhere in the sample (not shown). The negative sign indicates that the flux travels in conductive form from the hot to the cold side. Ohm Peltier-Seebeck  The transverse Righi-Leduc of Case II is due to the presence of a magnetic field B 1 perpendicular to the heat flux, as in Fig. 10. In this situation, the heat conductivity becomes a tensor entity as in (62) bottom. Two dof (temperature and magnetic potentials) are present. The magnetic field slightly alters the linear distribution produced by Fourier conduction. This alteration is only visible for very high values of B. As in other effects with presence of B 1 , the alteration is antisymmetric due to the cross product inherent to the magnetic field and more evident at the lateral sides x 2 constant. Figure 39 shows the heat fluxes along x 3 (direct, top) and x 2 (coupled, bottom). In the first, the Righi-Leduc is slightly relevant near the corners and in the rest of the domain the Fourier dominates with the value of q 3 given in the previous paragraph. In an infinite domain this alteration would be not noticeable. In the second, the induced heat flux along x 2 is maximum near the lateral faces. Not considering edge effects, its value is similar to the approximated q 2 ¼ ÀM B 1 DT n =l 3 ¼ À5:36 Â 10 3 [J/sÁm 2 ]. The maximum Righi-Leduc flux is about forty times smaller than that of Fourier. However, it can increase when the magnetic field is higher, and detrimental temperature concentrations will appear due to the finiteness of the geometry.
As in the rest of cases where B 1 is present, the magnetic field has been prescribed with an equivalent potential applied on the frontal x 1 -constant face with the parallel one prescribed to zero. To directly prescribe B 1 (for instance with an interface element, viz. [112]) would produce the same fluxes at the center but somewhat different at the corners, due to the curvature of the magnetic field.

Cases III and IV
The Ohm and Hall effects are qualitatively (but not quantitatively) equivalent to those of Cases I and II if temperature and heat flux are substituted by voltage and electric flux.
Ohm's law relates voltage gradient and electric flux through the electric conductivity, as shown in Fig. 8 and Table 3. From the third (63) and again not considering thermoelectric coupling, a ¼ 0, the law is obtained. The analytical flux is j 3 ¼ Àc DV=l 3 ¼ 2:15 Â 10 6 [A/sÁm 2 ], very similar to the constant one obtained with FEM everywhere in the sample (not shown). The Hall effect is a transverse interaction due to the presence of a magnetic field perpendicular to the electric flux, see Sect. 1.1.3. In this situation, the electric resistivity becomes a tensor equivalent to the second (62). The boundary conditions given Table 11 are used. Voltage distributions for these cases are shown in Fig. 40 with linear (left) or quasilinear (right) distributions. Note that the effect of the Hall effect is to disturb the isolines close to the x 2 -constant faces.
In the Hall effect Case IV, the variability of V generates vertical and horizontal electric fluxes with j 3 ¼ 2:36 Â 10 3 [A/m 2 ], positive since the voltage difference from 0 to 0.1 [V] is against the positive axis; this flux is almost two orders of magnitude smaller than that of Ohm's. The symbol jqj is the determinant of the tensorial  62) center. If in the equations the coupling R or the flux B 1 are zero, Ohm's law is recovered, therefore as an approximation the value of the scalar electric resistivity is directly taken as q ¼ 1=c from Table 11. The Hall voltage distortion produces a significant transverse Hall electric flux j 2 ¼ À7:1 Â 10 4 [A/m 2 ], larger than the -a priori-main vertical flux. The numerical distributions of these fluxes are shown in Fig. 41. Although good approximations between analytical and numerical fluxes in the center of the sample are observed (j 3 ¼ 2:4 Â 10 3 ; j 2 ¼ À7:1 Â 10 4 at the central node), strong concentrations of one order of magnitude appear at the edges. Again these concentrations are in diagonally opposite edges due to the rotational nature of the magnetic field.

Cases V and VI
As mentioned in Sect. 1, the Peltier, Seebeck and Thomson effects describe NEI between thermal and electric fields. Mathematically, the Seebeck effect is represented by the coupling term on the third equation (63); the Peltier by the coupling term in the last equation. Thomson will not appear in this case due to the use of constant properties.
Whereas the Seebeck effect can be uncoupled by specific boundary conditions, Peltier's cannot due to an additional Joule effect, last term in the last (56), created by the induced electric currents. As discussed in Sect. 1.1 the Joule effect is an irreversible interaction between fluxes: not described by material constants but implicitly present in the energy balance equation. Due to the impossibility of decoupling Peltier effect, Cases V and VI combine two effects to study Peltier-Seebeck and their Ettingshausen-Nernst transversal effects, respectively.
For Case V, hot and cold face temperatures are prescribed equal to eliminate direct conduction vertical thermal flux, see Table 11; an electric flux is directly created by the vertical variation of V. The distributions of temperature and voltage are quadratic, the first symmetric and the second almost due to a small voltage increase. It is not trivial to analytically calculate the direct or induced fluxes with the transport equations (63) since a Joule effect is also present in the balance (56) and rT n is locally non zero. Extensive thermoelectric verifications were done in [112,115] and will not be repeated here.
The relevant results of this case are shown in Fig. 42. In the left, the non-linear temperature distribution is due to Joule with an increment of 2.27 C concentrated in the middle. Had the complete temperature dependency of (92) being considered, the nonlinearities would have been more pronounced. in the right figure, the linear vertical flux is shown. The positive values are due to the plus sign of a, implying that heat is taken from the cold side (bottom) towards the hot side (top) in a refrigeration process, overcoming Joule in the middle. Ettingshausen-Nernst is a secondary effect again created by a transverse magnetic field B 1 . In this situation, the quadratic temperature distribution of Case V is significantly distorted and increased. The distortion is asymmetric in the two lateral sides, see Fig. 43 left, extraordinarily increasing Joule's effect in the left and decreasing the initial temperature in the right. The voltage distribution is nearly undisturbed and remains basically linear; the same can be said for the vertical heat flux, with a variation of 2 % in its maximum and minimum (both not shown). The temperature increment is important since also increments stresses and specially concentrates them, viz. the ongoing works [117,118] for detailed studies on mechanical couplings. A symmetric induced horizontal heat flux, Fig. 43 right, appears; is one third the vertical flux, and the fact that these fluxes are perpendicular to each other can significantly increase or decrease the performance of a thermoelectric device, viz. [145].
A symmetric and significant induced horizontal electric flux j 2 is created by this effect, see Fig. 44. It changes sign from the hot to the cold face and again the rotational nature of the magnetic field increases the effect in one lateral side.

Case VII
From a mechanical point of view, the Lorentz forces defined in (38) are body forces due the circulation of an electric current in the presence of a magnetic field perpendicular to it. As mentioned for Joule, this is an interaction present in the balance equation (56) and not in the material properties. In this equation these forces are represented by the last term in the first equation. They are analogous to the Maxwell stress tensor but present in conductive media. Also and as discussed previously, the Lorentz effect causes the galvanomagnetic and   Figure 45 shows the horizontal T C 2 Cauchy stress caused by an electric flow prescribed in the x 3 direction through potential gradient and by a magnetic field in the x 1 direction (see Table 11) with the plane x 3 ¼ 0 clamped. The deformation is increased six orders of magnitude, to visualize the distortion. As can be appreciated, the stress is small but can increase if the magnetic field is higher and the dimensions much smaller. In these situations, it would be important to consider Lorentz forces in the design of micro-devices.
The deformation is similar to that of a short cantilever block under a volumetric horizontal force against x 2 , that is, a combination of shear (prevalent) and bending. The stresses are negative (compression) close to the edge x 2 ¼ 0 and positive (traction) close to x 2 ¼ l 2 according to the sign of the force. As predicted by the theory of elasticity, the stress maximum is located where the slope of the deformation changes sign. At the built-in face these stresses are equal at the edges (maximum restraint) but of opposite signs for mechanical equilibrium in the x 2 direction.

Closure
In this article a formulation of multicoupled nonlinear finite elements has been presented. Four fields have been coupled: mechanical, electrical, magnetic and thermal, the last one both from equilibrium and non-equilibrium thermodynamics. Therefore, a multiphysics theory has been used, to develop multi-balance equations: balance of linear and angular mechanical momentum, balance of linear electromagnetic momentum, balance of electric charge, balance of energy and entropy.
A revision of the state-of-the-art of numerical studies (under small displacements and strains) of active materials is also included. In general, these studies do not couple more than two fields and are devoted to rather simple, from the theoretical point of view, practical applications.
A number of simplifications are introduced in the formulation, allowing the distinction between equilibrium and non-equilibrium interactions. Roughly, the former is applied for the study of non-conducting materials and the later for conducting materials.
Equilibrium interactions are conservative: there is no increment of entropy and the related constitutive equations can be derived from a potential, the electromagnetic enthalpy. The temperature, that by definition introduces irreversibilities, is introduced by a two-temperature model that assumes small increments. Nonequilibrium interactions are non-conservatives, and the related constitutive equations (also called transport equations) are deduced from the entropy production term of the system.
In both types of interactions non-linearities are present. In equilibrium, due to the Maxwell tensor and in nonequilibrium due to the Joule effect and the temperature dependency of the material properties. These interactions have been implemented with two special finite elements into the research code FEAP. Standard fully integrated eight-node, 7 degrees-of-freedom, isoparametric elements, Newton-Raphson and Newmark-b algorithms have been employed.
After formulating the tangent matrices, residual vectors and their related partial derivatives, a number of simple validations are presented. The objective of the article is not the study of real devices, but these cases permit the validation and assure the correct running of the elements.
The formulations and finite elements described are currently applied in ongoing works by the authors to more complicated applications, including pulse dynamics, further couplings and optimizations.
The EI derivatives in (72) to (74) have to be explicitly calculated and are listed in the first part of this appendix.
First, the derivatives of the Cauchy stress tensor from (72) are The first equation is a 6 Â 3 matrix; the second to fourth are directly 6 Â 1 vectors. The derivatives of the symmetric part of the Maxwell stress tensor from the same stiffnesses are The first equation cannot directly follow the matrix multiplication convection, as oT M s =oa U b in (72) must be interpreted as a 6 Â 3 matrix composed of three 6 Â 1 vectors. ! Therefore, in the right hand side the derivatives of D, B are also with respect to each mechanical dof with a resultant dimension of 3 Â 3, giving one of the components of the previous expression to be converted into the 6 Â 1 Voigt notation. For instance, for the second mechanical dof, the result is The rest of the Maxwell tensor derivatives are directly 3 Â 3 matrices. Again for (72), using the third of (71) the derivatives of the ponderomotive forces are where _ B Â , _ M Â can directly be calculated from (29). Similar to the first equation of this Appendix, the first derivative of f PM is implemented as three 3 Â 1 column-vectors. The other three derivatives are directly 3 Â 1 vectors. The subindex Â implies the allocation of a 3 Â 1 vector into a 3 Â 3 antisymmetric matrix, see Sect. 5 and (66).
The derivatives of D in the Maxwell stress expressions and in (73) are from its definitions 3 Â 1 vectors. Due to the form of the first two equations (71), these derivatives are equal to those of P present in f PM (third of (71)), except for the one with respect to a V b . The same must be done for the derivatives of _ D with respect to the time derivatives of the dof.
and those of D > or D Â (as B > ; B Â in the next equations) are directly the transpose or the cross form of the results. Some of these expressions will be used for the derivatives of the capacity matrices. The remaining derivatives are The derivatives of B and its time derivative for the stresses, ponderomotive forces and (74) are Since it is useful in a hysteresis simulation of ferroelectric materials, it has been assumed that the properties , l may vary with the corresponding field although not with time. In order to preserve the physical sense of laboratory tests, and since the matrices of , l are diagonal (see Appendix 4), the derivatives are performed term by term and may be stored in a 3 Â 1 vector; for instance the material permittivity derivative results in The remaining derivatives are The partial derivatives from capacities (76) represent heat dissipation due to electromagnetic dynamics; at the microscopic level, they can be related to friction among dipoles.
Finally, the derivatives of the two-way couplings from (77) using (50) are

Appendix 2
Related to NEI, the derivatives of the tangent stiffness tensors are obtained using the chain rule. For (81), the Cauchy tensor derivatives are equal to the first and fourth from EI in Appendix 1. The derivatives of the Maxwell vacuum tensor are from (63) oT Again Voigt notation may be applied to these derivatives to obtain a 6 Â 1 vector. The non-zero derivatives of the electric flux are Following the same procedure, the non-zero derivatives of the thermal flux are For both fluxes and from (60) the derivatives of the Peltier and conductivities' matrices are, for the symmetric parts In the previous expressions the equalities have been used, where O is any second-order tensor, symmetric or not. Similar expressions to the second equation can be deduced for q, j. The total derivatives da=dT n , dq=dT n , dj=dT n are directly calculated from their expressions in (92). The derivatives of B with respect to a U b , a T n b are zero since neither piezomagnetic nor pyromagnetic interactions are considered; in addition, due to the use of magnetic scalar potential oB=oa V b % 0. The only non-zero term is As mentioned in Sect. 4.2.2 for a T , the material properties reflect the transverse isotropy of Fig. 18. The exception in the listed data is in p V , p u , for which properties are assumed isotropic due to the lack of data.