A novel mixed and energy-momentum consistent framework for coupled nonlinear thermo-electro-elastodynamics

A novel mixed framework and energy-momentum consistent integration scheme in the field of coupled nonlinear thermo-electro-elastodynamics is proposed. The mixed environment is primarily based on a framework for elastody-namics in the case of polyconvex strain energy functions. For this elastodynamic framework, the properties of the so-called tensor cross product are exploited to derive a mixed formulation via a Hu-Washizu type extension of the strain energy function. Afterwards, a general path to incorporate nonpotential problems for mixed formulations is demonstrated. To this end, the strong form of the mixed framework is derived and supplemented with the energy balance as well as Maxwell’s equations neglecting magnetic and time dependent effects. By additionally choosing an appropriate energy function, this procedure leads to a fully coupled thermo-electro-elastodynamic formulation which benefits from the properties of the underlying mixed framework. In addition, the proposed mixed framework facilitates the design of a new energy-momentum consistent time integration scheme by

yields bending type deformations.In the latter, however, Coulomb forces are responsible for the electrically induced deformations of the EAPs, which can result in more complex shape changes compared to their ionic counterparts.
Among electronic EAPs, dielectric elastomers (DEs) have demonstrated remarkable electrically induced actuation properties facilitated by their lightness, fast response, biomimetism, and low stiffness properties. 2,3DEs are indeed capable of exhibiting massive electrically induced deformations, as demonstrated by experimental studies, 4 where electrically induced area expansions of 1962% in DEs have been reported.As a result, DEs have been identified as ideal candidates for their use in the field of soft robotics, where cutting-edge technological developments start evidencing the onset of a paradigm switch: Traditional hard robotic systems will be replaced by soft robotic solutions, 5 especially in applications requiring safe interactions with humans or, when human avatars are demanded in extremely hazardous environmental conditions. 6However, the applications of DEs are not limited to the field of electrically induced actuation, as DEs have been successfully applied as Braille displays, deformable lenses, haptic devices, and energy generators, to name but a few. 7dvances in the field of soft robotics can be accelerated through high fidelity finite element simulations for which a reliable material characterization of these materials, and of DEs in particular, is of major importance.Recently, Mehnert et al. 8 devoted great effort to the mechanical and coupled electro-mechanical characterization of so-called very high bond (VHB) polymers.This characterization was extended to account for the thermo-viscoelastic behavior in these materials.The experimental works in References 8-11 enabled a comprehensive characterization of the thermo-electro-viscoelastic properties of particle filled silicone EAPs.Crucially, the characterization conducted on these experimental studies has led to subsequent works with deep roots in continuum mechanics, thermodynamics, and even micromechanics.In Reference 12, the analytical strain energy function is formulated based on thermo-electro-mechanical invariants and internal variables (accounting for viscoelastic effects).A calibration of the analytical constitutive model was carried out by fitting their response to the data generated in the experimental work using optimization techniques.Crucially, as proven by the aforementioned works, the performance of DEs is highly sensitive to temperature variations and its consideration is of paramount importance.This is for instance the case, when understanding the dynamics of DE-based soft robotic underwater applications, 13 where the temperature of the surrounding fluid might vary.Accordingly, these findings clearly justify the necessity of embedding temperature effects into the constitutive models of DEs.
Analytical constitutive models should not only approximate the available data of the material response (whether available from experiments or analytical computations or computational homogenization), but also comply with physical requirements, that is, objectivity, material symmetry and so forth. 14,15Further mathematical assumptions have to be embedded by realistic constitutive models.It is well assumed that these must satisfy the ellipticity or rank-one convexity condition. 16This condition prevents the formation of extremely localized deformation modes, guarantees the propagation of real travelling plane waves in the material, 17 and ensures the well-posedness of the underlying boundary value problem. 18A sufficient condition complying with the ellipticity condition is the polyconvexity 16,19 of the strain energy function (analytical constitutive model).][22] In addition to reliable and physically-mathematically sound constitutive models, long-term time dependent finite element simulations of DEs require the use of stable and robust spatial and temporal discretizations.4][25][26][27] The underlying reason for this lies in their thermodynamic roots, as they are endowed by construction with the discrete analogue of the conservation properties of the continuum, namely the conservation of total energy, total linear momentum, and total angular momentum.Consistency of these methods, namely their ability to conserve (or dissipate for nonreversible constitutive models) the total energy of a system in agreement with the laws of thermodynamics, 23 is attained by replacing the (exact) partial derivatives of the strain energy function (or other types of energy like, for example, the internal energy function) with respect to its arguments with their carefully designed algorithmic counterparts.
Recently, Betsch et al. 28 proposed a novel EM scheme in the context of nonlinear elasticity, by taking advantage of the concept of polyconvexity and the use of a novel tensor cross product pioneered by de Boer 29 and re-discovered in the context of nonlinear continuum mechanics by Bonet et al. 30,31 In essence, the authors in Reference 28 proposed the consideration of three discrete derivatives which were used to form an algorithmic version of the second Piola-Kirchhoff stress tensor.These three discrete derivatives represent the algorithmic counterparts of the work conjugates of the right Cauchy-Green deformation tensor, its cofactor, and its determinant.This strategy leads to a simplified expression of the algorithmic second Piola-Kirchhoff stress tensor, compared to those obtained by the classical approach. 23,24Later, this work has been extended to multiphysics scenarios such as thermoelasticity, 32,33 nonlinear electro-mechanics, 34,35 and more recently to thermo-electro-mechanics. 36 In Reference 36, a simple thermo-electro-mechanical constitutive model is 10970207, 2023, 10, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7209by Swansea University Information, Wiley Online Library on [01/06/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License considered in which the thermal effects are coupled exclusively with the volumetric deformations, excluding any direct interaction between thermal and electrical effects.The present work aims at incorporating more generic and sophisticated constitutive models.In addition, it advocates for a multifield formulation where not only displacements, electric potential, and temperature are regarded as unknown fields.Rather, strains, electric displacement, and suitable stress-type Lagrange multipliers, are incorporated as unknown fields, relying on a stable spatial discretization of all fields involved.In that regard, a noteworthy contribution lies in the introduction of a general path to incorporate nonpotential problems for mixed formulations, where this path is not limited to the considered area of thermo-electro-elastodynamics.
The manuscript is organized as follows: Section 2 describes the multiphysics initial boundary value problem (IBVP) associated with the underlying governing equations for DEs in nonisothermal scenarios.Section 3 introduces the necessary elements of algebra associated with the tensor cross product operation.In this context, the tensor cross product is used to conveniently rewrite the cofactor of the deformation gradient.Furthermore, this section presents the analytical constitutive model considered.In addition, a Hu-Washizu type variational formulation for the purely mechanical elastostatic case is presented.From its stationary conditions, the associated multi-field strong form is derived.This strong form serves as a starting point for the subsequent extension to the nonconservative multiphysics case of interest, where electrical and thermal effects are incorporated.Section 4 presents the weak form of the multi-field multiphysics IBVP presented at the end of Section 3. Additionally, the conservation properties of the proposed continuous formulation are examined.Section 5 illustrates the algorithmic treatment of the weak form in Section 4, where derivatives of the analytical constitutive model are replaced with their algorithmic counterparts to embed the desired conservation properties into the resulting EM time integration scheme.Section 6 presents the spatial discretization of the algorithmic weak form shown in Section 5. Eventually, numerical tests are conducted throughout Section 7 with the aim of testing the spatial and temporal convergence properties of the proposed multi-field EM time integrator, as well as its long-term stability.Finally, concluding remarks are provided in Section 8.

NONLINEAR CONTINUUM THERMO-ELECTRO-MECHANICS
We consider a deformable continuum body  with boundary  subject to time t ∈  = [0, T], where T ∈ R + .We distinguish between the reference configuration  0 ⊂ R 3 and the current configuration  t = ( 0 , t), where  ∶  0 ×  → R 3 denotes the bijective mapping that maps material coordinates X ∈  0 to its spatial counterparts x = (X, t) ∈  t .

Initial boundary value problem
The IBVP for thermo-electro-mechanics is obviously composed of contributions from the mechanical, the thermal, and the electrical fields and is introduced in the following.

Elastodynamics
The basic kinematic quantity in elastodynamics is the deformation gradient which maps infinitesimal line elements dX ∈  0 to their spatial counterparts dx ∈  t .The IBVP for elastodynamics is given by The above set of equations is comprised of the balance of linear momentum (2) 1 , where  0 ∶  0 → R is the material mass density, the dot is the first Piola-Kirchhoff stress tensor, and B ∶  0 ×  → R 3 denotes the prescribed volume force per unit undeformed volume  0 .Furthermore, the Neumann boundary condition is given by (2) 2 , where N ∶  T  0 ×  → R 3 denotes the material outward normal vector to the Neumann boundary  T  0 ∈  0 , and T ∶  T  0 ×  → R 3 is the prescribed Piola stress vector.Equation (2) 3 denotes the Dirichlet boundary condition with the prescribed displacement field  on    0 ∈  0 .The whole boundary for the mechanical field is comprised of both introduced boundaries, which may not overlap, such that Finally, (2) 4 -(2) 5 are the initial conditions for configuration and velocity with  0 and v 0 being the initial configuration and velocity, respectively.

Thermodynamics
The IBVP for thermodynamics is given by where (4) 1 denotes the local energy balance with  ∶  0 ×  → R + being the absolute temperature field and  ∶  0 ×  → R is the entropy density.Furthermore, R ∶  0 ×  → R denotes a prescribed heat source per unit undeformed volume  0 and Q ∶  0 ×  → R 3 denotes the Piola heat flux vector per unit undeformed volume  0 .Equations (4) 2 and (4) 3 are the Neumann and Dirichlet boundary conditions, respectively.Therein, Q ∶  Q  0 ×  → R denotes a prescribed rate of heat transfer across a unit of undeformed area applied on  Q  0 ⊂  0 .Essential temperature boundary conditions  are applied on    0 ⊂  0 .The boundaries for the thermal field in (4) need to satisfy Eventually, the initial thermal field is prescribed in (4) 4 with the initial temperature  0 .

Electrostatics
Restricting our focus to capacitor-like DEs allows us to neglect magnetic and time-dependent electrical effects.Thus, it is possible to reduce Maxwell's equations to Therein, D 0 ∶  0 ×  → R 3 represents the Lagrangian electric displacement vector,  e 0 ∶  0 ×  → R is the prescribed electric volume charge per unit of undeformed volume  0 , and  e 0 ∶    0 ×  → R denotes the electric surface charge per unit of undeformed area    0 ⊂  0 .Furthermore, E 0 ∶  0 ×  → R 3 represents the Lagrangian electric field vector and Φ ∶  0 ×  → R is the electric potential field which is prescribed with Φ on  Φ  0 ⊂  0 .The boundaries are composed

Constitutive equations
The behavior of body  0 is characterized by the chosen energy density function.For many applications it is important to incorporate coupling effects between the fields, that is, to consider interactions between the mechanical, thermal, and electric fields.In doing so, we herein assume an energy function  =  (F, D 0 , ) (cf.Reference 36).Considering the second law of thermodynamics, the constitutive equations must hold for the first Piola-Kirchhoff stress tensor P, the entropy density , and the Lagrangian electric field vector E 0 .Furthermore, the constitutive equations are supplemented by the Piola heat flux which is governed by Duhamel's law of heat conduction.Therein, K is the material thermal conductivity tensor and  ∶  0 ×  → R 3 denotes the material gradient of the absolute temperature field Assuming a suitable energy  , and K to be a positive semi-definite second-order tensor, the constitutive equations ( 8) and ( 9) satisfy the Clausius-Duhem inequality and can thus be regarded as thermodynamically consistent.

ADVANCED MODELING ASPECTS WITHIN THE CONTINUOUS DESCRIPTION
In this section, we outline advanced modeling aspects within the continuous description.We want to emphasize that many of these aspects are not necessarily limited to the thermo-electro-mechanical system considered in Section 2, but can be applied to other problems as well.More precisely, we present a polyconvexity-inspired framework based on the tensor cross product operation, thermo-electro-mechanical constitutive models, a Hu-Washizu type mixed formulation for nonpotential (coupled) systems and an alternative energy balance formulation to facilitate the design of the desired energy-momentum consistent time integration scheme.Eventually, we propose a multiphysics, mixed IBVP for thermo-electro-elastodynamics.

Polyconvexity-inspired framework
In order to find a polyconvexity-based energy density function describing the constitutive behavior of a DE, the tensor cross product operation between second-order tensors is introduced.Afterwards, kinematics and constitutive equations are reconsidered by the beneficial properties of the tensor cross product.

Tensor cross product
The tensor cross product was introduced in Reference 29, exploited in the context of nonlinear continuum mechanics by Bonet et al. 31 and subsequently extended to multiphysics problems (cf.References 32-36).It enables to simplify algebraic derivations, especially when polyconvex constitutive models are considered.
The tensor cross product between two second-order tensors A, B ∈ R 3×3 can be defined as Therein, Einstein's summation convention is employed to pairs of repeated indices, where i, j, , , a, b ∈ {1, 2, 3}.In addition,  ijk denotes the third-order permutation tensor.A noteworthy feature of the tensor cross product is that it allows to redefine the cofactor of a second-order tensor A ∈ R 3×3 as Furthermore, the determinant of a second-order tensor A ∈ R 3×3 can be written in an elegant manner as Some further helpful properties of the tensor cross product can be found, for example, Reference 36.

Kinematics
In addition to the deformation gradient F in (1), alternative kinematic measures can be defined.More precisely, its cofactor H ∶  0 ×  → R 3×3 can now be defined with the help of the tensor cross product as in (12), namely This allows to map infinitesimal material area elements dA ∈  0 to their spatial counterparts da ∈  t .Furthermore, with the help of equation ( 13), the Jacobian determinant of F can be written in terms of the tensor cross product as The Jacobian determinant maps infinitesimal volume elements of the reference configuration dV ∈  0 to their spatial counterparts dv ∈  t .For an objective (material frame indifferent) representation of the energy density function, we introduce the symmetric right Cauchy-Green strain tensor In addition to that, the cofactor of C can be defined in a similar fashion as in (14), that is, where G ∶  0 ×  → R 3×3 .Eventually, we introduce the determinant of C in analogy to (15), that is, where C ∶  0 ×  → R + .

Energy density function
We assume that the behavior of DEs can be characterized by means of the energy density function where Example 1. Inspired by the works, 37,38 we exemplarily consider the following additive decomposition of the energy density function  (C, G, C, D 0 , )* † as where In the above, the mechanical parameters a, b, c ∈ R + and d = 2 (a + 2 b) have been employed.Moreover, the permittivity of vacuum is denoted by  0 = 8.8541 × 10 −12 AsV −1 m −1 and  r ∈ R + is the relative permittivity of the considered medium.Furthermore, following.36 Therein, , e ∈ R + are thermo-mechanical parameters and  R ∈ R + is the reference temperature.In addition,  t () ∶ R + → R represents the purely thermal contribution, 36 defined as where  ∈ R + denotes the specific heat capacity.Finally, the factor following the work of Reference 38.Notice that multiplication of f  and  em aims at a fully coupled model (cf.References 38 and 37), which contrasts with our previous publication, 36 where f  = 1.Accordingly, the constitutive equations are given by where S ∶  0 ×  → R 3×3 denotes the second Piola-Kirchhoff stress tensor.In addition, assuming thermal isotropy, Duhamel's law of heat conduction (9) reduces to Fourier's law of heat conduction.Accordingly, with K = k 0 I we obtain where k 0 ∈ R + denotes the coefficient of thermal conductivity with respect to the reference configuration.
Remark 1.The proposed constitutive model ( 21) within Example 1 can be considered as polyconvex only up to a certain temperature  cr .More precisely,  must satisfy the conditions for this to be true.However, as can be seen in Reference 36, the critical temperature is typically relatively high, so we will not address this topic further.
Remark 2. Mehnert et al., 37 considered a similar additive decomposition to that in equation ( 21) within Example 1, where the main difference resided in the higher nonlinearity associated with the temperature dependent factor f  (), defined therein as where a 1 , a 2 ∈ R are material parameters.Such nonlinearity was introduced with a two-fold purpose.First, to introduce a dependence of the specific heat capacity with respect to deformations.This can be seen from the definition of the specific heat capacity, that is, Second, a definition of f  () as in ( 29) ensures a nonlinear dependence of both stresses and electric fields with respect to temperature, which, according to equation (8), would depend on a linear fashion with respect to f  (), comprising the linear term ∕ R and the nonlinear term g().

Mixed formulation for nonpotential coupled systems
In order to develop a multi-field, mixed formulation for nonpotential systems, we carry out two steps.First and foremost, we start with a variational Hu-Washizu type formulation for elastostatics (cf.Reference 28) and convert it to its associated (mixed) strong form.In a second step, we augment the (mixed) strong form to a mixed thermo-electro-mechanics strong formulation.
A noteworthy feature of the above procedure is that it is not restricted to thermo-electro-mechanical coupled problems but does also work for various coupled systems including nonpotential contributions.

Variational Hu-Washizu type formulation for elastostatics
For the elastostatic case under consideration, the potential energy can be defined as with its internal and external contributions respectively, where Ŵ ∶ R 3×3 → R denotes the strain energy.In the above, the solution function  ∈ V  is employed with space where H 1 denotes the Sobolev space.The principle of stationary potential energy states that the first variation of the potential energy disappears if equilibrium prevails (cf.Reference 40).The principle of stationary potential energy is equivalent to the procedure of deriving the weak form directly from the strong formulation of a problem.By introducing the objective kinematic set within a cascade manner § as the potential energy can be re-expressed as an augmented 7-field Hu-Washizu type version where the above kinematic set is enforced (cf.Reference 28).Therein, the kinematic set  = (C, G, C) and the set of Lagrange multipliers  = ( C ,  G , Λ C ) have been introduced to simplify notation.Furthermore, the solution functions C, G,  C ,  G ∈ V A , with generalized space and C, Λ C ∈ V A , with generalized space are used, where S denotes the space of symmetric second-order tensors and L 2 is the space of square integrable functions.Imposing stationarity conditions to (35), that is, variation with respect to the independent variables and setting the variation to zero yields the mixed variational formulation for elastostatics where the equations above have to hold for arbitrary  ∈ V 0  , with C, G,  C ,  G ∈ V 0 A , and C, Λ C ∈ V 0 A .The variational formulation (38) can be converted to its strong form using basic algebraic operations.This strong form consists of the equations Div(2 and the corresponding boundary conditions

Mixed thermo-electro-mechanical strong formulation
To extend the BVP in ( 40)-( 41) to the case of thermo-electro-elastodynamics, we proceed as follows.First, we take into account the inertia term within the balance of linear momentum, that is, in equation ( 40) 1 .Second, we consider the IBVP of thermodynamics ( 4) and the BVP of electrostatics ( 6) to incorporate thermal and electrical effects.Finally, we employ the energy density function  =  (C, G, C, D 0 , ) from ( 21) to achieve coupling between the different physical processes.The above steps yield the final multiphysics, mixed IBVP ¶ consisting of the equations the boundary conditions

Energy formulation to facilitate discrete energy consistency
In order to facilitate the design of an energy-momentum consistent time integration scheme, the local energy balance on a continuum level in (42) 3 is reconsidered as which was already proposed in Reference 34 and applied in the context of thermo-electro-mechanics in Reference 36.

WEAK FORMULATION
Based on the strong formulation given in (42), but exchanging the local energy balance in (42) 3 with ( 45), the associated symmetric weak formulation of ( 42) is obtained after multiplying each equation with a suitable test function, integrating over the domain  0 , and applying basic algebraic operations.This yields where the external thermal and electrical potentials are given by  and C, Λ C ∈ V A must hold for the solution functions, with the remaining function spaces defined by

Balance laws
In this section, the balance laws of the coupled problem given by ( 46) are examined.In particular, total angular momentum and total energy are considered.To this end, a homogeneous Neumann problem is assumed by applying

Conservation of angular momentum
The total angular momentum of the continuum body with respect to the origin of the inertial frame is given by To verify conservation of angular momentum, we take the following admissible values for the test functions (46) using arbitrary but constant values for  ∈ R 3 .With the above, (46) 1 yields Accordingly, the conservation of total angular momentum can be written as Furthermore, we introduce the skew-symmetric matrix ζ ∈ R 3×3 by  × a = ζa with property ζT + ζ = 0, such that (46) 2 with (52) can be written as In the above the total external mechanical torque exerted by body and surface loads is introduced as Hence, for vanishing external mechanical loads, the total angular momentum is conserved.The total energy of the thermo-electro-mechanical body consists of the total kinetic energy and (cf.Reference 36) Accordingly, conservation of total energy is achieved if Therein, the derivatives of the energies with respect to time can be written as In order to verify the conservation of total energy, we replace the test functions (w v , w  , w  , w Φ , w D 0 ) in (46) with and obtain Furthermore, we apply (46), which leads to Eventually, we use (w (46), which yields Finally, we add equations (60) 2-5 to (61) 1-3 and subtract (60) 1 and (62) 1-3 , respectively, such that we obtain the desired result Accordingly, for vanishing external loads, the total energy is a constant of the motion.

TIME DISCRETIZATION: ENERGY-MOMENTUM SCHEME
In this section, we conduct the time discretization of the weak equations ( 46) derived in the previous section.For that, we divide the interval  into subintervals [t n , t n+1 ] ⊂ .In doing so, we assume that the values of the fields , v, , Φ, D 0 , , and  are known at time t n and unknown at time t n+1 .For completeness, we denote the fields at time t n by (•) n and at time t n+1 by (•) n+1 .Furthermore, we introduce the average value and the time step size The objective is to develop a numerical one-step method that permits to compute the values of the unknown fields at time t n+1 .For this purpose, we choose the time discretization for the weak equations (46) 1−5 .In addition, we propose for the weak equations (46) 6−8 and for the weak equations (46) 9−11 .
In (66), the time-discrete versions of the external mechanical, thermal, and electrical contributions are given by where Therein, ⟨•, •⟩ denotes the inner product and Remark 3. If  i is a scalar quantity, the definition of the discrete derivatives is simplified.It is then reminiscent of Greenspan's formula (cf.Reference 45) and becomes As shown in detail in Reference 36, the discrete derivatives possess two important properties: 1. Directionality: This property is crucial to show the conservation of energy and states that 2. Consistency: One can show that the discrete derivatives are well defined in the limit ||Δ i || → 0, where Remark 4. Note that D  i  is not defined for ||Δ i || = 0.This case may occur, for example, during the first iteration of Newton's method.To overcome this potential numerical difficulty, it is reasonable to compute the discrete derivatives for ||Δ i || → 0 as for the midpoint rule, that is, The validity of (75) follows from (74).

Semidiscrete balance laws
In this section, we show that the proposed time integration scheme in equations ( 66)-(68) satisfies the conservation laws for any time step size.As customary in doing so (cf.Reference 32), we assume the homogeneous Neumann case.

Conservation of angular momentum
The total angular momentum at time t n+1 and at time t n is defined as The total angular momentum is conserved if its variation from t n to t n+1 vanishes, namely To show that the time-discrete version of the weak form satisfies this requirement, we assume that , where  ∈ R 3 is arbitrary but constant.Inserting this value into (66) 1 results in where ζ is again a skew-symmetric matrix and Considering the properties of ζ, we finally arrive at Thus, as long as M m,ext n+ 1 2 = 0, for example, when there are no external mechanical torques acting on the body, we can conclude that the proposed formulation conserves the total angular momentum.

Conservation of energy
The total energy at time steps t n+1 and t n can be written as with the corresponding kinetic energies and The total energy is conserved from t n to t n+1 if its variation vanishes, that is, To show that the proposed time-discrete formulation conserves the total energy, we start by assuming that in equation (66

and thus obtain
In addition, we replace the test functions (w  , w Φ , w D 0 ) in (66 Inserting the constitutive relations and subtracting (66) 4 from (66) 5 then yields The enforcement of the constraints on position level in the equations of (68) can be transferred to the velocity level, as shown in Reference 28.More precisely, if the equations in (68) hold at both time t n and time t n+1 , one can show that also holds.Next, we consider the directionality property of the discrete derivatives (73) as well as the relations (87) and (88).By making the additional assumption that the test functions (w  C , w  G , w Λ C ) in (89) take the admissible values and Inserting the spatial approximations into the time-discrete weak form in (66)-(68) yields where A n e e=1 denotes the assembly operator and where the nodal residuals are defined as

R a,e
In the above, B a ∈ R 6×3 is the standard nodal operator matrix and (•) V refers to Voigt's vector notation of symmetric stress-related quantities.Note that we have omitted the superscript (•) h in the representation of the residuals (97), for clarity.

NUMERICAL EXAMPLES
In this section, we investigate the properties of the novel mixed thermo-electro-mechanic formulation numerically in a variety of examples.Therein, three different finite element families are employed for space discretization.These are: The nodal points of the different finite element mesh types are given in Figure 1 for convenience.We employ the enhanced Mooney-Rivlin material model proposed in Example 1 for the numerical investigations.In that regard, the material parameters provided by Table 1 are used for all subsequent examples.

Patch test
The purpose of the widely used patch test is to verify that a finite element formulation is capable of correctly reproducing homogeneous states of stress -even with distorted meshes.This is a fundamental condition that any admissible finite element formulation has to fulfill.To show that the novel mixed formulation passes this patch test requirement, we consider a cube-shaped body  0 = [0, 1] 3 m 3 , which has an initial temperature of  0 = 293.15K and is investigated with two different finite element meshes (see Figure 2).The initial coordinates of the inner nodes of the distorted mesh can be

F I G U R E 2
The regular (left) and the distorted (right) mesh employed in the patch test.
found in Reference 46.We employ H2 c H1 d elements for both the regular and the distorted mesh.For the sake of simplicity, equivalent results for the H1 c H0 d and P2 c P1 d elements are not shown in Figures 2-5.
The simulation is performed for a total time T = 1 s with a time step size Δt = 0.1 s and a tolerance of Newton's method of eps= 1 × 10 −6 .During the simulation, the body is held on three faces to prevent translational and rotational movements and is compressed displacement-driven to half of its original height.These specifications can be expressed with respect to the initial configuration by the mechanical Dirichlet boundary conditions and The compression causes the body to expand in the X 1 and X 2 directions and to heat up, which is due to the chosen material model.Note that the patch test presented here is a quasi-static example, since we neglect the transient effects for the mechanical and the electrical part, but consider them for the thermal part.Thus, the proposed integrator is active and the discrete gradients are in use.Furthermore, we choose  0 = 0. We subsequently investigate, whether the formulation is able to reproduce a homogeneous state of stress as well as a homogeneous temperature distribution and a homogeneous

F I G U R E 5 Electric potential distribution Φ [V]
resulting from the patch test for the regular and the distorted mesh.
distribution of the electric potential, as it can be expected for the patch test.We use von Mises stress  vM as stress measure, where the required Cauchy stress tensor can be obtained by As Figures 3-5 show, the novel formulation satisfies the patch test requirement.More precisely, a simulation with the given parameters leads to a homogeneous stress, temperature, and electric potential distribution for both the regular and the distorted mesh.

Analytical convergence analysis
The idea of this example is to numerically verify the spatial order of all fields for the finite element families H1 c H0 d , P2 c P1 d , and H2 c H1 d (cf.References 35,36,47).For this ad-hoc manufactured example all time effects are neglected, such that we have a thermo-electro-elastostatic problem with stationary heat conduction.Accordingly, the proposed integrator is not employed for this example.We again consider a cube-shaped body  0 = [0, 1] 3 m 3 .For the analytical convergence analysis we assume solutions of the form for the primary fields, where Γ k = 0.01 k, k ∈ {1, 2, 3}, θ = 10 K, and Φ = 100 V.These assumed solutions allow to compute the first Piola-Kirchhoff stress tensor, the Lagrangian electric displacement vector, and the Piola heat flux vector # analytically.Eventually, we compute the body force, charge density, and heat source analytically via the strong form given by (42) as where all time effects are neglected.The above analytically computed values are subsequently used for the numerical simulations.In that regard, Dirichlet boundaries are imposed at the faces of the body for all fields in agreement with the assumed analytical solutions (101).In particular, considering Figure 6, outer faces (patterned surfaces) of the body are imposed with Dirichlet boundary conditions by  a , Φ a , and  a , given by (101).The simulation is performed for a total time T = 1 s with a time step size Δt = 0.1 s and a tolerance of Newton's method of eps= 1 × 10 −6 .Snapshots of the final mesh with von Mises stress, electric potential distribution, and absolute temperature field are depicted in Figure 7 for P2 c P1 d and H2 c H1 d elements.Eventually, the analytical solutions in (101) for all fields (•) a are compared to the numerically computed solutions (•).For this purpose, the h-convergence rate is computed by employing the L 2 norm of the error e (•) , given by where To this end, the h-convergence results are depicted in Figure 8.As expected, at least p + 1 convergence is observed for all fields for the considered elements.

Rotating cross-shaped body
The objective of the first transient example is to illustrate the conservation properties as well as the order of accuracy of the novel mixed formulation.For this purpose, we consider a cross-shaped body, whose geometry and finite element mesh is given in Figure 9.The finite element mesh consists of a total of 104 H2 c H1 d elements.We further specify electrical Dirichlet boundary conditions of the form  where Apart from that, no other Dirichlet boundary conditions are imposed.However, due to a prescribed initial velocity the body rotates around the X 3 -axis during simulation.The system is simulated for a total time T = 10 s with a time step size Δt = 0.05 s.Furthermore, we set the tolerance of Newton's method to eps= 1 × 10 −3 and the initial temperature of the body to  0 = 293.15K. F I G U R E 9 Geometry and finite element mesh of the cross-shaped body.
A simulation with the given parameters yields the total angular momentum evolution depicted in Figure 10.Looking at the differences of the total angular momentum from one time step to the next, it is evident that the novel mixed framework is consistent with respect to the total angular momentum for both the midpoint rule (MP) and the energy-momentum scheme (EM).As Figure 11 shows, this does not apply to the total energy.While the energy is conserved for t > 0.5 s (dashed line) when using the energy-momentum scheme, using the midpoint rule results in an energy blow up.
Next, we investigate the order of accuracy of the energy-momentum scheme.For this purpose, we compute the L 2 norm of the error, which accumulates in the time interval 0.5 s ≤ t ≤ 0.6 s, for different time step sizes and for all primary variables (, , Φ).The L 2 norm of the error is computed according to (103).However, for this example, we replace the analytical solution with a reference solution that is computed with a very small time step size (Δt = 2.5 × 10 −5 s) due to the lack of an analytical solution.As Figure 12 shows, the energy-momentum scheme inherits, as expected, the second-order accuracy of the midpoint rule (slopes of the regression lines are (, Φ, ) = (2.0449,2.0127, 2.0242)).

Microfluidic pumping device
The purpose of the final example is to evaluate the proposed energy-momentum scheme within a realistic long-term application.
A sophisticated technical realization of temperature sensitive DEs are microfluidic pumping devices.Such microfluidic pumping devices are employed, for instance, as medical implants, for example, for micro injection of drugs (see

F I G U R E 10
Left: Evolution of the total angular momentum resulting from the rotating cross-shaped body example over the course of the simulation computed using both the energy-momentum scheme (EM) and the midpoint rule (MP).Right: Incremental differences of the total angular momentum from one time step to the next.

F I G U R E 11
Left: Evolution of the total energy resulting from the rotating cross-shaped body example over the course of the simulation computed using both the energy-momentum scheme (EM) and the midpoint rule (MP).Right: Incremental differences of the total energy from one time step to the next.Reference 37).The microfluidic pumping device to be considered is closely related to the one presented in Reference 48 and 37 and likewise does not consider fluid-structure interaction.The full shape of the cylindrical pumping device with H2 c H1 d mesh is shown in Figure 13.
Radii and dimensions of main body and outlet nozzle and boundary conditions of the pumping device are provided by Figure 14.According to the symmetry only one eighth of the device with corresponding symmetry boundary conditions is considered during the simulation (see also Figure 14).
For the pumping purpose two elastomer layers are sandwiched between two compliant electrodes at top and bottom walls, respectively (see Figure 14, right).In particular, the pumping motion is controlled by a Dirichlet boundary condition of the electric potential field via where Φ 0 = 0.75 GV and .., n cyc (see also Figure 15).Furthermore, the heating of the pumped fluid is idealized by a thermal Neumann boundary condition imposed at all internal walls and controlled by The simulation is performed with a total simulation time of T = 100 s, a time step size of Δt = 0.1 s, and a tolerance of Newton's method of eps= 1 × 10 −6 .Accordingly, more than three pumping are performed.
A mesh for one eigth of the device is comprised of 320 H2 c H1 d elements with a total of 9290 displacement, electric potential, and absolute temperature unknowns (size of the system of linear equations to be solved after static condensation) is employed.
Snapshots of the motion with electric potential, absolute temperature, and von Mises stress results are shown in Figure 16 for half of the pumping device.Due to the chosen setting, large deformations, temperature, and electric potential evolutions can be observed.After the heating phase and during constant phases of the applied potential fields (see Figure 15) the total energy is a constant of motion.To this end the total energy and the energy difference are shown in Figure 17.As can be observed the energy is perfectly reproduced by the proposed EM integration and is numerically stable during the long-term simulation with periodic loading.Furthermore, the energy difference is bounded by the chosen tolerance of Newton's method.

CONCLUDING REMARKS
In this article, we have a novel framework for coupled thermo-electro-elastodynamics.Building upon previous works in the field of elastodynamics, 28 electro-elastodynamics, 34,35 and recently upon thermo-electro-elastodynamics, 36 we have proposed a new tailor-made multi-field, mixed finite element formulation for the multiphysics context considered.It is worth mentioning that the approach we have chosen to address the problem is not limited to thermo-electro-elastodynamics, but can be applied to design mixed formulations for nonpotential (multi-field) systems in general.
The starting point is a variational mixed formulation for elastostatics that is subsequently converted to its associated mixed strong form.In a second step, the mixed strong form is then extended to a mixed coupled strong form by supplementing the equations with the desired physics in form of additional initial boundary value problems and a suitable coupling within, for example, the energy density function.These two steps can be applied for a variety of different systems, including those with nonpotential contributions.
Our present work is based on a mixed elastostatic formulation, which is specifically designed for the case of polyconvex strain energy functions (cf.Reference 28).This formulation is particularly elegant due to the cascading introduction of independent strain fields and the use of a rather unknown algebraic operation, namely the tensor cross product.We were able to extend the benefits of this formulation to the case of thermo-electro-elastodynamics by following the aforementioned path.In this context, we have incorporated a state of the art polyconvexity-based thermo-electro-mechanical constitutive model that characterizes the behavior of dielectric elastomers.This model embeds a full interplay between the three physics involved, namely, thermo-electro-mechanics.This is in contrast with a more simplified constitutive model considered in our previous publication, 36 where the thermal field was coupled exclusively with volumetric deformations.
Furthermore, we have shown how the advantages of energy and momentum consistent time integration schemes can be transferred to multiphysics problems, in particular thermo-electro-elastodynamics.
Finally, we have included a series of numerical examples to investigate the proposed framework.More precisely, we have evaluated the time and space convergence properties as well as the long term stability of the formulation.In doing so, we were able to emphasize the advantageous characteristics of our formulation.
The source code used for the finite element computations is implemented in Matlab under MIT license and is available at https://github.com/kit-ifm/moofeKIT.Version 1.00 of the code, used in this article, is archived at Reference 49.
ENDNOTES * The model given in (20) is typically known as a modified entropic electro-elasticity model.
† It is important to remark that the proposed framework is not restricted to the suggested energy function given in (20).
‡ The polyconvexity properties of (23) are discussed in more detail in Appendix A. 39 § G is bounded to independent variable C instead of C() = (∇ X ) T ∇ X .The same applies for the independent variable C which is bounded to the independent variables G and C, respectively instead of G(C()) and C().¶ By additionally enforcing constraint (6) directly, we also yield a mixed formulation with regard to the electric part, which is known as 'hybrid finite-element model' (see Reference 41) and commonly applied in electrostatics due to some advantages (cf.References 35,36,42).# A detailed survey of the analytical computations is provided in Appendix B. || Note that summation convention applies to pairs of repeated indices in (A1).
which replaces (97) 1 and (97) 2 , such that we eventually obtain the vector of nodal residuals For the sake of simplicity, the specific contributions of K e  , K e  and so forth are not provided here.In particular, K e This static condensation procedure eliminates the number of unknowns for the system of linear equations to the minimal set Δq e  with the reduced residual vector Re and the associated tangent matrix Ke .Accordingly, the computational effort to obtain the solution is significantly reduced.After computation of Δq e  , we are able to compute Δq e  via relation (A13).

APPENDIX B. ANALYTICAL CONVERGENCE ANALYSIS EXAMPLE: ANALYTICAL COMPUTATIONS
The assumed analytical solutions in (101) build the starting point of the analytical computations, that we need to provide for the numerical example in Section 7.2.In the following, we provide the related computations in order to compute the source terms which are neccessary to plug in (102): • With (1), ( 16), (17), and (18) we are able to compute the kinematic relations 2 ) e 2 ⊗ e 2 + (1 + 2 Γ 3 X 2 3 ) e 3 ⊗ e 3 , (B1) ) ) • With ( 6), (10) we compute the gradients of the electric potential and the thermal field E a 0 = −∇ X Φ a (X),  a = ∇ X  a (X). (B5) • With ( 21) and E 0 =  D 0  we are able to compute the Lagrangian electric displacement vector ) .

2 w
Φ dA.The proposed discretization in time conforms to an EM scheme.An essential part of this scheme is the utilization of the so-called discrete derivatives D C  , D G  , D C  , D D 0  , and D   , which were introduced in Reference 43 and can be interpreted as algorithmic or time-discrete counterparts of  C  ,  G  ,  C  ,  D 0  , and    .According to Reference 44, the concept of discrete derivatives is a further development of the original idea from Reference 23, limited to St. Venant-Kirchhoff material.The five discrete derivatives can be defined, using the abbreviated notation 

F I G U R E 1
Nodal points of the different finite element families.The nodes of the continuous fields are represented by bullets •, whereas the nodes of the discontinuous mixed fields are indicated by crosses ×.

TA B L E 1
Material parameters employed for the numerical examples.

F I G U R E 3 F I G U R E 4
Von Mises stress distribution  vM [Pa] resulting from the patch test for the regular and the distorted mesh.Temperature distribution  [K] resulting from the patch test for the regular and the distorted mesh.

F I G U R E 6
Dirichlet boundary conditions (left), initial P2 c P1 d (mid), and H2 c H1 d (right) meshes employed for the analytical convergence analysis example.F I G U R E 7 Von Mises stress  vM [Pa] (left), electric potential Φ [V] (mid), and absolute temperature  [K] (right) for P2 c P1 d (top), and H2 c H1 d (bottom) meshes, respectively.

10970207, 2023, 10 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/nme.7209by Swansea University Information, Wiley Online Library on [01/06/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E 8 Spatial convergence analysis for H1 c H0 d (left), P2 c P1 d (mid), and H2 c H1 d (right) elements for the analytical convergence analysis example.

F I G U R E 12 L
2 norm of the error resulting from the rotating cross-shaped body example, plotted for different time step sizes and all primary variables.FI G U R E 13Mesh of the microfluidic pumping device.

F I G U R E 14
Top view (left) and side view (right) including mechanical (symmetry lines and hatching), electrical (Φ 1−4 ), and thermal (Q) boundary conditions employed for microfluidic pumping device example.F I G U R E 15 Heating of the pump with Q (left) and electric potential fields Φ 2 and Φ 3 for the pumping process (right) employed for microfluidic pumping device example.F I G U R E 16 Snapshots of the microfluidic pumping device with electric potential field Φ [V] (left), the absolute temperature field  [K] (mid), and von Mises stress  vM [Pa] (right) at t = 0 s (top), t = 6 s (mid), and t = 22 s (bottom).

F I G U R E 17
Figure 15 (right).