A computational framework for magnetically hard and soft viscoelastic magnetorheological elastomers

This work deals with a comprehensive theoretical and numerical framework that allows the modeling of finite strain agnetorheological elastomers (MREs) comprising mechanically soft nonlinear elastic–viscoelastic polymer phases and agnetically hard (i.e. dissipative) or soft (i.e. purely energetic) magnetic phases. The framework is presented in a general anner and is implemented using the finite element method. Two software implementations are developed, one using FEniCS and the other in Abaqus. A detailed analysis of the numerical schemes used to model the surrounding air is made and their pros and cons are discussed. The proposed framework is used to simulate two model geometries that are directly relevant to recent applications of MREs. The first two-dimensional example simulates a mechanically soft beam consisting of a single wavy-chain of hard or soft magnetic particles. The beam is subjected to transverse magnetic actuation loads that induce important vertical deflections. Despite the overall small local strains in the beam, a significant viscoelastic effect is observed when high-rate magnetic fields are applied. A torque model for the particles is also used to analyze the beam geometry and is found to be in relatively good agreement with the rest of the approaches for small actuation fields. The second example discusses the rotation of a three-dimensional ellipsoid embedded in a cubic elastomer domain, while the ensemble lies inside a larger cubic air domain. Non-monotonic uniaxial and rotating magnetic fields are applied leading to complex, non-monotonic rotations of the ellipsoidal particle. The hard and soft magnetic cases exhibit significant differences, whereas viscoelasticity is found to induce strong coupling with the magnetization rotation but not with the dissipative magnetization amplitude. Extensive supplementary material provides all details of our implementations as well as animated visualization of results. © 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Magnetorheological elastomers (MREs) or magnetoactive polymers (MAEs) are most often designed as composites comprising a soft elastomer matrix carrying micron-sized magnetic inclusions. While both of their ingredients, elastomers and magnetic particles, have been available for many decades, it took surprisingly long until MREs were originally fabricated and studied. In fact, landmark studies on establishing a rigorous theoretical framework can be already found in the middle of the twentieth century [1][2][3]. Further cornerstone contributions accounting also for electrodynamics and small strain magnetoelasticity are due to Pao and Hutter [4], Maugin [5], Eringen and Maugin [6], James and Kinderlehrer [7], DeSimone and James [8] among others. These early theoretical contributions on general coupled magneto(electro)mechanics predate the first well-known experimental works on MREs [9]. In a related work [10], the same authors provided dipole-interaction-based models for the explanation of some of their experimental findings. Despite these early achievements, it took more years for MREs to find their place in the mechanics and materials communities. The works of Dorfmann and Ogden [11], Dorfmann and Brigadnov [12], Dorfmann and Ogden [13], Kankanala and Triantafyllidis [14] mark the rise of continuum magnetoelasticity, specially designed at larger scales relevant to MREs, in the context of which one also has to mention Kovetz [15] as well as Ericksen [16] and Bustamante et al. [17]. These re-initiating contributions sparked a vast amount of works fostering the understanding of MREs both at the microstructural (scale of the particles but not below at the level of magnetic domains) and the macrostructural level. The former is obviously important since magnetic particle-particle interactions and their resulting effect on the overall material response depend on the microstructural distribution of particles in the composite. One approach to efficiently capture the effect of these interactions on the macroscopic behavior is homogenization, where we mention the works [18,19] on analytical as well as [20][21][22][23][24][25] on numerical homogenization, respectively. Nonetheless, only recently such homogenization efforts [19] have led to practical macroscopic material models for magneto-elastic composites [26,27]. In the same context we also refer to Kalina et al. [28] and the more recent works based on dipole interaction by Garcia-Gonzalez and Hossain [29].
In addition to the above-mentioned microscopic effects, purely macrostructural magnetomechanical coupling plays an equally important role for the overall response [30,31] and in particular underlying instabilities [32,33]. Indeed, these omnipresent structural effects render a major challenge in the experimental determination of material properties [34][35][36]. This is one more reason for the importance of multiscale studies to understand the multitude of interacting effects in MREs.
What has largely been neglected by most research studies up to now is the effect of viscosity of the elastomer matrix on the overall magnetomechanical response on the one hand and the effect of ferromagnetic hysteresis at the particle and macroscopic scale on the other hand. As notable exceptions, we mention [37][38][39][40] 1 for the former and [41,42] for the latter effect. Keip and Sridhar [43] addressed the point of dissipative ferromagnetism at the particle level with the framework of micromagnetics [44]. However, while this approach provides great insight into the fine-scale details of domain wall motion, its computational cost renders it prohibitive for the simulation of macroscopic MREs.
We also mention the recent contributions considering already magnetized elastic [45][46][47] and viscoelastic [48] bodies. However, neither of them accounts for the evolution of magnetization which sets physical bounds on the actual magnetization of bodies and thus is also of great practical importance. Hence, there are still a number of open tasks in the modeling of both viscoelasticity and magnetic hysteresis and, in particular, their interplay. Motivated by this and also by recent application-oriented developments [49,50] on the so-called hard 2 MREs (h-MREs), the present work aims to contribute to the understanding of the effect of filler particles exhibiting ferromagnetic hysteresis on an otherwise nonlinear elastic-viscoelastic elastomer. For that purpose we present a thermodynamically consistent computational framework based on the formalism of generalized standard materials [51,52] accounting for viscoelasticity and ferromagnetic hysteresis inspired by Miehe et al. [53] and Kumar and Lopez-Pamies [54].
The implementation of the material models is carried out using the newly developed open source material modeling toolbox Materiaux [55] paired with FEniCS [56,57] as well as a user-element (UEL) code based on Abaqus. The two implementations are provided via zenodo.
The numerical examples in the present work discuss individual particles in a mechanically soft, but otherwise non-magnetic matrix and a slender structure made of stiffer elastomer comprising a chain of spherical magnetic particles. In either case, once remanently magnetized, an h-MRE sample will behave much differently from a s-MRE (i.e. comprising magnetically soft particles, such as iron, exhibiting negligible dissipation). To our best knowledge, there does not exist any macroscopic continuum or homogenized material model for nonlinear elastic-viscoelastic h-MREs. The examples in the present contribution are thus chosen to explicitly account for a simplified particlematrix structure of MREs. Nevertheless, the proposed framework and accompanying software provide the tools to obtain such models, e.g., when combined with numerical homogenization techniques along the lines of Lefèvre et al. [26] and Mukherjee et al. [27] on s-MREs and Mukherjee et al. [42] on h-MREs.
This work is structured as follows. The fundamentals of the magnetomechanical theory are presented in Section 2. The discretization of the magnetomechanical initial boundary value problem as well as numerical and computational aspects are covered in Section 3. In Sections 4 and 5, we provide a set of representative numerical simulations of magnetic particles carried by an elastomer and subject to complex magnetic loads. These examples demonstrate the effect of viscosity and ferromagnetic hysteresis compared with purely elastic and magnetically non-hysteretic responses. For that purpose, we apply appropriate material models with realistic parameters and also account for the free space (air) around the bodies under consideration. A special Section 4.1 is devoted to the modeling of the surrounding air. Section 6 concludes this contribution and provides an outlook to future work.

Theory
In this section, we introduce the building blocks of the computational framework. We begin with the fundamentals of magnetomechanics at finite strains before presenting the magnetomechanical initial boundary value problem (IBVP) considered in this work. In the very last subsection we present the material models employed, which are representative for viscoelasticity and ferromagnetism. Thermodynamic consistency of the models is achieved by rigorous application of the framework of generalized standard materials (GSM; [51]).
Given the complexity that will arise simply as a result of the combination of models for various phenomena, we strive to keep the discussion of each physical contribution rather compact and provide instructive references along the way.

Fundamentals of magnetomechanics
In this subsection, we briefly introduce the fundamental fields and equations of magnetomechanics based on the pioneering works of Tiersten [2,3] and Brown [1] whose theory has been specialized for magnetoactive (magnetorheological) elastomers by Dorfmann and Ogden [11] (see also equivalent formulations by Kankanala and Triantafyllidis [14] and Ericksen [16]).
We consider the deformation y of a body V embedded in R 3 . The current position x ∈ V t of each point with initial coordinates X ∈ V is given by x = y(X, t) and the displacement field is given by u(X, t) = y(X, t) − X. We furthermore assume to have a smooth extension of y from V into R 3 such that the notions of "current" and "initial" (configurations, positions) on V extend to all space. For convenience in this work, we take the initial configuration as reference and accordingly refer to quantities defined on this initial configuration as referential or Lagrangian to distinguish them from their current or Eulerian counterparts. We follow the widely adopted convention to use uppercase symbols and operator names to indicate their referential nature. Based on y we have the tangent map known as deformation gradient, which is the fundamental ingredient for strain measures.
On the side of magnetostatics in the case of vanishing currents, we may employ the auxiliary magnetic scalar potential ϕ(X) as a primary variable of which the negative gradient is the referential magnetic h-field H, i.e.
This definition readily implies that Curl H = 0, i.e., Ampere's law for magnetostatics in absence of free currents is automatically fulfilled. The choice for ϕ(X, t) as primary magnetic variable in absence of free currents is motivated by the computational efficiency of the resulting finite element formulation and for this reason is widely used in the computational magneto-mechanics community. Nonetheless, alternative formulations based on the vector potential [17], which are particularly efficient in the context of instabilities, have been proposed in the last years such as those by Kalina et al. [24] and Danas [25] in two dimensions as well as by Miehe et al. [23] and Dorn et al. [58] in three dimensions. In this context we also highlight the work of Semenov et al. [59] on the numerics of vector potential formulations for three-dimensional electromechanics, that are perfectly parallel to those for magnetomechanics. The corresponding magnetomechanical balance equations read Div S + f 0 = 0 and Div B = 0.
The first equation corresponds to the balance of linear momentum in terms of the total first Piola-Kirchhoff stress S and the mechanical body force density per unit referential f 0 . The second one is the magnetic Gauß law in terms of the referential magnetic field B.
In the present framework, the internal Helmholtz free energy density 3 W is assumed to depend on F and H but also on a set of internal variables denoted for the moment as I. As a result, the local dissipation density D is given as the difference of the external power P minus the rate of change of the internal energyẆ , i.e., Owing to the arbitrariness ofḞ andḢ, one may apply the standard Coleman-Noll-Gurtin argument to arrive at the well-established constitutive relations We note that the energy density W (F, H, I) is of energy type with respect to F(u) and I but of co-energy type with respect to H(ϕ). Thus, solutions of corresponding boundary value problems are of saddle-point nature. While this is a drawback with respect to purely energetic formulations, the greater numerical efficiency renders the present approach favorable in most scenarios. For a complete energy type formulation for dissipative magnetomechanical systems the reader is referred to Kalina et al. [41] and Mukherjee [61].

Remark 1.
It is recalled here that the Eulerian fields, σ (Cauchy stress), b and h may be written in terms of their Lagrangian counterparts by a push-forward operation leading to the well-known equations In this configuration, one may also define the Eulerian magnetization vector m, which is non-zero only in a magnetizable body, as Here, µ 0 is the magnetic permeability in vacuum to be specified below. The current m is used to describe the magnetic response of solids and will be used extensively in our work too. Nonetheless, we note that in the present work, it does not constitute an independent variable but is merely a post-processed field in the deformed configuration and is obtained by Eq. (7).

The abstract magnetomechanical initial boundary value problem
In this section, we explicit the equations introduced in the previous sections in the context of an initial boundary value problem (IBVP). The (referential) domain of the problem is denoted as V ∞ . In the most general case, V ∞ = R 3 but for finite element computations, one typically resorts to truncated domains As different V i may have different (in nature) material properties and internal states, we use i as a superscript to distinguish between them. The fields u(X, t) and ϕ(X, t) are defined for all (X, t) ∈ V ∞ × T , where T denotes the total time interval. From u(X, t) and ϕ(X, t), one directly obtains F = 1 + Grad u and H = − Grad ϕ according to Eqs. (1) and (2), respectively. 4 The equations to be satisfied by u(X, t) and ϕ(X, t) are given by (3). In the present work, for simplicity, we ignore mechanical body forces such as gravity. Hence, we have Div S(F(u), H(ϕ), I) = 0 and Div B(F(u), H(ϕ), I) = 0, 3 What we call Helmholtz free energy density is sometimes referred to as amended in the sense that it also contains magnetostatic contributions from non-magnetic media, e.g., air, or just vacuum [17,60]. It is thus not bound to the presence of matter. 4 Henceforth, we will write F(u) and H(ϕ) when we insist to emphasize from which states {u, ϕ} the derived quantities F and H are computed or when we rather regard them as operators, e.g., in variational principles. as governing equations for u and ϕ, where S(F, H, I) and B(F, H, I) are obtained from (5). In addition to (8), we impose Dirichlet boundary conditions on u and ϕ π (u(X, t)) = π (u(X, t)) for X ∈ S u (9) ϕ(X, t) = ϕ(X, t) = −X · H ∞ (t) for X ∈ ∂V ∞ .
Here, π denotes a projection operator such that one can apply conditions only on certain components of u on a surface S u and We note that, in the present work, for simplicity, we apply boundary conditions for ϕ only at the outer boundary ∂V ∞ ; nonetheless more general conditions may readily be considered as in [58,62]. In addition, guided from practical problems at hand, we also consider trivial Neumann-type boundary conditions on S and B (see, e.g. [22,63]). The framework presented thus far is identical to that of quasistatic non-dissipative magnetomechanics, which may be written in terms of a variational principle that corresponds to (8) and (9) [17,30].
In turn, in the context of dissipative magneto-mechanics additional equations are required. Specifically, given the fields u and ϕ that comply with the field equations and boundary conditions, the rate of internal statesİ is governed by the typical evolution equation deriving from the generalized standard materials (GSM) formalism. These equations take the form for smooth and non-smooth dissipation potentials as respectively, where ∂İ is the subdifferential and the dissipation potential D is convex inİ [64]. It is worth noting at this point that some domains in V ∞ have no internal state, e.g. the empty space (or air) surrounding a magnetomechanical solid. Nevertheless, there is still magnetic energy present even in free space and needs to be taken into account in actual IBVPs. The previous set of field and constitutive equations and boundary conditions can be combined in a variational principle for the unknowns of the problem {u(X, t), ϕ(X, t)} (short: {u t , ϕ t }) alongside with the evolution equations for the internal variables I i (X, t) (short: respectively. In the first expression, the admissible sets K(t) and G(t) read The initial values of the internal states I i t=0 have to be specified explicitly, whereas u t=0 and ϕ t=0 result from (12) at t = 0. The evolution of internal states affects the constitutive responses such that not only the I i but also the fields u and ϕ may undergo changes even under constant boundary conditions. Thus, the evolution of internal states is indeed coupled with the external primary fields. The individual energy densities W i as well as the dissipation potentials D i will be specified further in Section 2.3.

Remark 2.
In subdomains V i that correspond to empty space or air, the deformation map requires special treatment [22,33,63]. For details we refer to Section 3.3.
Remark 3. The Euler-Lagrange equations corresponding to (12) and the evolution equation (13) can also be obtained from a variational principle along the lines of Hackl [65], who describes the corresponding procedure for finite-strain elastoplasticity.
The magnetic load is applied in terms of an external field b ∞ t . The outer domain boundary ∂V ∞ is sufficiently remote from ∂V such that the magnetic field perturbations caused by the presence of V p to practically vanish at ∂V ∞ . Note further that the Eulerian field b ∞ t is applied and not the Lagrangian one, which may be obtained if required by use of (6).

Specialization to a particle-matrix system
In order to proceed further with more concrete examples, we specialize our study to the important case of mechanically rigid, hard-magnetic (i.e. dissipative) particles (occupying the domain V p ) embedded in a soft nonmagnetic, nonlinear elastic-viscoelastic elastomer matrix (occupying the domain V m ). The ensemble (occupying the domain V = V m ∪ V p ) is then immersed in an air domain V a = V ∞ \ V, all forming the overall domain V ∞ . This problem corresponds to "mesoscopic" magnetorheological elastomeric structures and materials. While a macroscopic treatment without the resolution of individual particles would be computationally more efficient, we do not take this path yet, as there exist currently neither enough experimental data nor sufficiently general continuum or homogenized models for hard magnetic nonlinear elastic-viscoelastic MREs. 5 Moreover, the goal of the present study is to provide a more general numerical framework including the air domain. This framework then may be easily specialized to conduct homogenization calculations without air and periodic boundary conditions similar to the studies of Javili et al. [66], Miehe et al. [23] Kalina et al. [67], Mukherjee et al. [27] and Mukherjee et al. [42].
In this regard, the fairly generic setting under consideration that specializes the abstract domain V ∞ = ⋃ i V i from the previous section is shown in Fig. 1.
Standard isotropic invariants. Without loss of generality with respect to the general formulation, we consider that the materials under study are isotropic both mechanically and magnetically. For this reason and to keep our analysis tractable but at the same time general, we focus on a minimum set of isotropic invariants which are then used to propose energy densities for each phase. The magnetomechanical invariants given below ensure material frame indifference and isotropy of the energy density functions. Specifically, we retain the isotropic mechanical invariants with denoting the right and left Cauchy-Green tensors, respectively. Furthermore, we consider the magnetomechanical invariants 6 Additional invariants of similar form will be introduced later to deal with the corresponding dissipative contributions below.
In what follows, we specify the material models employed for the air, the elastomer matrix and the magnetic particles.

Air domain
The "air" domain V a = V ∞ \ V only has a magnetostatic vacuum energy density (per referential volume), which reads This energy is a function of I 5 and J only, which implies that the solution for h is not affected by the local deformation of any non-magnetic medium. The corresponding Lagrangian magnetic field and total first Piola-Kirchhoff stress as well as their Eulerian counterparts are obtained as and respectively, where 1 is the second-order identity tensor. In order to avoid severe mesh distortion in the air, especially in the present context of large transformations of MREs, a proper definition of the deformation map y(X, t) and thus F in V a is required. This constitutes a highly non-trivial issue that is discussed in detail in Section 3.3.

Elastomer matrix
The elastomer matrix occupies the domain V m (see Fig. 1) and is assumed to follow a (quasi-)incompressible viscous constitutive law. We employ a simplified version of the generic GSM-based viscoelasticity approach of Kumar and Lopez-Pamies [54]. In that work, the authors rigorously demonstrated the conditions that are necessary for material frame indifference and material symmetry and those details are not repeated in this study. In addition, using that framework, Ghosh et al. [69] have recently proposed an explicit, analytical homogenization model for nonlinear elastic-viscoelastic particle-reinforced elastomers and as a result the proposed numerical framework in this study can be readily adapted in future implementations to include such homogenization models. Notwithstanding, any other viscoelastic model, e.g., [70][71][72][73][74], may be readily used if required by relevant experimental data.
Specifically, we use a viscous elastomer material model of a Zener-Maxwell-type as depicted in Fig. 2. This description assumes a multiplicative decomposition of the total deformation gradient F in terms of an energetic, elastic part F e and a dissipative, viscous part F v , such that For later use, we define next some important quantities derived directly from the previous multiplicative decomposition. In particular, we introduce the corresponding right and left Cauchy-Green elastic tensors C e and B e , as well as the right and left Cauchy-Green viscous tensors respectively, as This further implies that invariants of C e can be expressed as invariants of CC v-1 . In direct analogy to (16), this leads us to the invariants It is noted from the second equation above that J = J e J v > 0 is a constraint on the various volumetric total, elastic and viscous parts and its implications will be discussed later in this section. The matrix energy density is then assumed to derive from an additive decomposition of an equilibrium part Ψ eq , a non-equilibrium part Ψ neq and the vacuum or "background" magnetic contribution present in non-magnetic media Ψ vac from (19), such that W m (I 1 , J, I e 1 , I 5 ) = ρ 0 Ψ eq (I 1 , J ) + ρ 0 Ψ neq (I e 1 , J e ) + Ψ vac (I 5 , J ), where I 1 and J are given by (16), I e 1 and J e by (24). Subsequently, we choose to work with hyperelastic-type energy functions of a form given by [75] ρ 0 Ψ eq (I 1 , and In this expression, both Ψ eq and Ψ neq are considered quasi-incompressible with G ′ ≫ ∑ r G r and g ′ ≫ ∑ r g r . A more detailed discussion of this point is carried out below.
Following Le Tallec et al. [76] and [54], we write the dissipation potential as where d v is the viscous strain-rate (i.e. the gradient of the velocity field) and d v d = d v − 1/3 tr d v 1 denotes its deviatoric part. The two coefficients η K and η J denote the shear and bulk viscosities, which in the present work are assumed constant but more complex descriptions may be readily considered as discussed in [54].
Remark 4. The dissipation potential (28) may be shown to be quadratic inḞ v [54] and thus be strictly convex for constant η K > 0 and η J > 0. As a result the second law of thermodynamics is trivially satisfied, i.e., ( It is noted that as a direct consequence of the definition of d v in (28), the symmetric C v is chosen in the present work as the internal variable and thus an evolution equation forĊ v needs to be determined from the GSM relation Obviously, other choices are possible. However, one should always make sure to satisfy the second law of thermodynamics mentioned previously. After some cumbersome algebra that is provided in Appendix B, we obtain in the compressible case, the following evolution equatioṅ This equation results in a compressible C v and thus J v ̸ = 1. Again, we recall the constraint J = J e J v . Nevertheless, in problems involving elastomers, one may simplify this last equation further by assuming incompressibility of C v . When more general quasi-incompressible solids are considered, such an assumption has only practical interest since it implies taking the limit η J → ∞ and use η K to fit any experimentally relevant data. In this incompressible limit, the evolution Eq. (30) becomeṡ This last relation leads readily to an incompressible (up to numerical approximations to be discussed in Section 3) C v and thus J v = 1. As a result, J = J e irrespective of the hydrostatic parts considered in the energy W m . This further implies that if J ̸ = 1 then one needs to choose a potential W m such that it allows J e ̸ = 1. Otherwise, the kinematic constraints resulting from the multiplicative decomposition (22) cannot be satisfied. For this reason, in the present study, we will consider a large but finite g ′ in the definition of W m in (27) thus allowing J e = J ̸ = 1 in general (but fairly close to unity). The value of G ′ in (26) and g ′ in (27) of course do not need to be the same.
Non-dissipative limit. This corresponds to η → 0, i.e. the dissipation potential vanishes, but then C v → C such that B e → 1, which causes the non-equilibrium contribution to vanish as well. Representative response. The response of a representative model obtained from the proposed viscoelastic framework is shown in Fig. 3. In the graph we observe three different regimes. For low stretch ratesλ ≤ 1.0 s −1 (blue and orange curves) the response is practically non-dissipative and governed by the equilibrium energy density alone. For intermediate stretch rates 10 1 s −1 ≤λ ≤ 10 2 s −1 (green and red), one can observe significantly different responses during loading and unloading indicating a considerable viscous dissipation of energy. For even faster stretch loadinġ λ ≥ 10 3 s, the area between the loading and unloading path decreases again because the loading is too fast for any significant viscous relaxation to happen. Hence, the response of the model at hand is governed by the sum of the equilibrium and the non-equilibrium energy density, whereby the latter acts like just another elastic contribution because of negligible viscous relaxation and thus negligible dissipation. We refer to Section 3.1.1 for details on the algorithms used to actually compute the response.

Ferromagnetic particles
The ferromagnetic particles (V p in Fig. 1) considered in this work are assumed to be mechanically rigid (even though that is not necessary in the context of the present general framework) by comparison with the elastomer matrix. Thus, they may undergo rotations but no significant deformation [41]. In view of this, we employ the model of Mukherjee et al. [42] 7 in the limiting case of rigid particles, i.e. c → 1 in their terminology and some main features of this work are briefly recalled here for completeness. Specifically, the model is based on a magnetic internal vector-valued state variable H r that is defined at an intermediate stretch-free configuration. The corresponding Eulerian and Lagrangian counterparts are then given by with R and U being the rotational and the stretch components of F, i.e. The fields h r and H r contribute to the "full" fields h and H, respectively, as which implicitly defines the energetic h e and H e . In the present setting of (quasi-) rigid particles, i.e. U ≈ 1, this intermediate configuration coincides with the initial configuration. In order to ensure material frame indifference and isotropy, the model detailed below is given in terms of the invariants I 1 , J and I 5 already defined in (16) and (18), respectively, and in terms of the additional invariants Since U = C −1/2 , the former invariant can also be expressed as I er 5 = h · RH r . The abstract form of the GSM relation for the dissipative ferromagnetic model reads ∂ W p ∂H r (I 1 , J, I r 5 , where ∂Ḣr denotes the subdifferential [64,77] with respect toḢ r and B r is introduced as the energetic dual to H r .
The energy density W p is assumed to be of the form where Ψ mech is chosen to follow a simple Neo-Hookean law since the particle is fairly rigid, such as For the particles to be practically rigid in comparison to the elastomer matrix, one simply needs to choose sufficiently large values for the moduli G p ≫ ∑ r G r and G ′ p ≫ G ′ . Specific values are reported later in the results sections. The magnetic energy density Ψ mag reads ρ 0 Ψ mag (I r 5 , I 5 , I er 5 ) = ρ 0 Ψ e,mag (I 5 ) + ρ 0 Ψ r,mag (I r 5 , I er 5 ) where χ e is the magnetic susceptibility in absence of dissipative ferromagnetic evolution whereas χ r corresponds to the susceptibility added 8 due to dissipative ferromagnetic evolution. Therefore, we refer to them as the "energetic" and "remanent" susceptibilities of the particles. The saturation parameter m s coincides with the saturation magnetization if χ e = 0. The function f , controlling the magnetic saturation behavior, is chosen of the form so that it is sufficiently smooth and twice differentiable for all x > 0. Notably, f ′ (x) leads to an inverse saturationtype response which tends to ∞ as x → 1. This property, in turn, leads to a saturation type ferromagnetic response, which will be shown later in this section. We note at this point that we omitted the coupling energy density that leads to intrinsic magnetostriction and is present in the more general model of Mukherjee et al. [42]. The reason is simply that the magnetostrictive strains in the present mechanically-stiff particles are negligible. The vacuum energy density Ψ vac is that given by (19).
Next, we define the (non-smooth) dissipation potential D p, * in terms of its Legendre-Fenchel dual where the magnetic switching surface Φ 9 is given as with b c being the coercivity of the magnetically hard particle. The actual dissipation potential D p is then obtained from the dual potential as whereΛ is the Lagrange multiplier 10 enforcing Φ(B r ) ≤ 0. The optimality conditions of (43) arė of which the latter three are known as the Karush-Kuhn-Tucker (KKT) conditions. Inserting W p from (37) and D p from (43) in the generic evolution principle (13) yields the variational principle with stationary conditions governing the evolution of the internal state and the auxiliary variables B r andΛ. For later use, we introduce the objective function of the optimization problem (43) as such that 8 Their relation to the total susceptibility is actually more complicated than simply adding them up. In the limit of very small χ e , the total susceptibility is at first order equal to χ e + χ r . For a more detailed description we refer to Mukherjee et al. [42]. 9 It is perhaps helpful to note that the switching surface in magnetism is treated in a similar fashion to the yield surface in elasto-plasticity. The reader is referred to Mukherjee and Danas [78] for more details. 10 Note that we writeΛ following a convention borrowed from plasticity whereΛ is referred to as (plastic) consistency parameter and the "dot" indicates thatΛ indeed has the physical unit of a time rate. However, when identified as a Lagrange multiplier, there is no physical significance of a time integral ofΛ. In this sense we point out that in many works on GSM and also at some instances in plasticity theory, e.g. Lubliner [79],Λ is simply written as Λ.  Non-dissipative limit. The magnetic response in the non-dissipative limit of the above constitutive behavior corresponds to the limit b c → 0 and may be written as withf Specifically, considering b c = 0 we obtain the condition B r = 0. Hence,Ḣ r vanishes (see (46) 2 ) leading to D p (Ḣ r ) = 0, i.e., no dissipation. Finally, substituting B r = 0 into (46) 1 leads to an algebraic relation between H r and H. After some simplifications in accordance with the assumption of quasi-rigid particles we arrive at (49), i.e. an expression for Ψ mag only in terms of I 5 .
Representative response. The uniaxial response of a realistic ferromagnetic material is depicted in Fig. 4. For details on the algorithms that have been employed to compute the response, we refer to Section 3.1.2.
For extensions of the present model that allow to capture more complex magnetic responses in the case of negligible deformation we refer to Mukherjee and Danas [78]. In this regard we also mention the possibility to superpose the responses of a set of instances of the model above in order to fit experimental data [41,[80][81][82].

Discretization and numerics
In this section, we present the discretization of the magnetomechanical IBVP as well as certain algorithmic details. We first specify the updating algorithms for the viscoelastic and magnetic internal states. Then, we briefly summarize the global time-discrete problem. Subsequently, we turn to the spatial discretization aspects using the finite element method.

Time discretization
In what follows, subscripts "n" and "n + 1" indicate quantities evaluated at the time steps t = t n and t = t n+1 respectively. In addition, superscripts "k" and "k +1" correspond to iterations for the solution of a nonlinear problem for {u, ϕ} and internal states. For example, u k+1 n+1 refers to the value of the deformation map u at the (k+1)-th iteration for the computation of solution of the global IBVP. Note that we use u instead of u k+1 n+1 and so on when we refer to variables rather than concrete solutions.
For simplicity we restrict ourselves to lowest order approximations of rates. In particular, we introducė In what follows, any "dotted" quantity is to be understood as time-discrete.

Time discretization of the viscoelastic (rate-dependent) evolution in the elastomer domain
For the time discretization of the viscous evolution Eq. (31), we employ a Crank-Nicolson scheme. Using i.e. the right-hand-side of (31), the time-discrete evolution equation reads as We note that this is a nonlinear equation for C v that brings along a certain numerical effort. It can be solved by a standard Newton-Raphson algorithm in a straight-forward manner. Depending on the size of the boundary-value problem, the additional cost required for the solution of the local (spatially point-wise) nonlinear-equation is more or less negligible when compared to that of the overall solution time of a global nonlinear boundary value problem, e.g., a FEM simulation. As a simpler alternative, we have also successfully employed the explicit fifth-order Runge-Kutta scheme discussed in Kumar and Lopez-Pamies [54] and Ghosh et al. [69] and the differences were found to be negligible for the examples considered here.
Remark 5. We note that for the given incompressible dissipation potential D m , which in the present model explicitly depends on C v , it is not possible to derive the second-order accurate Crank-Nicolson scheme from an incremental variational principle, but only a semi-implicit first-order scheme. Nonetheless, we have found that the implicit scheme is inaccurate (in particular with respect to viscous incompressibility) and leads to substantial convergence problems. In any case, in the context of incompressible problems, the fifth-order Runge-Kutta explicit scheme is found to be robust and easily applicable with regard to maintaining unity for J v .

Time discretization of the ferromagnetic (rate-independent) evolution in the magnetic particles
In the case of the ferromagnetic rate-independent evolution, we face an optimization problem with an inequality constraint. A fully implicit discrete version of (44) (or (46)) using ∆Λ = ∆tΛ is represented by the system Since the switching surface Φ(B r ) is a function of B r only, the evolution Eqs. (54) can be obtained as the stationary conditions of the time-discrete version of (45) and (48), i.e.
respectively. In contrast to the viscoelastic evolution, the scheme selected for magnetic evolution is of first order only. This choice is justified by the non-smoothness of the evolution problem at hand. Both the "direct" (54) and the "optimization" (55) ferromagnetic evolution problem can be solved by active-set algorithms [83] of which the well-known return-mapping schemes [84,85] first devised for plasticity are a special case. In the present case, we opt for a rather generic approach outlined in Algorithm 1. Note that in this algorithm the condition on the trial state ensures that ∆Λ ≥ 0 thus allowing to drop the inequality constraint and employ a standard Newton-Raphson scheme to solve the nonlinear evolution equations. In a computer implementation, one might check whether ∆Λ ≥ 0 is indeed fulfilled afterwards. Moreover, B r and ∆Λ are only auxiliary results of which the history is not of concern. Storing them is thus not necessary but might save some iterations through improved initial "guesses" depending on the actual nonlinear solver scheme and implementation. Algorithm 1: Abstract scheme for the local ferromagnetic evolution at an outer iteration k + 1 For the validity of the optimization approach (55), it is required that the function L p does not explicitly depend on H r but only through the rateḢ r . While this holds in the present setting, it might not hold automatically for more complicated ferromagnetic evolution models. Any explicit dependence on H r then has to be substituted by H r n . We remark also that the evolution equation for the internal magnetic state in terms of F and H (and history states) is practically independent of F in the present setting of quasi-rigid particles where C ≈ 1.

Global time-discrete problem
With the material models for the individual domains or phases defined in Section 2.3 and the discrete evolution laws given in Sections 3.1.1 and 3.1.2, we provide next the time-discrete IBVP for the particle-matrix-system under consideration. We consider the generic IBVP given by (12) and (13) as a starting point. In the present setting, those lead to the (incremental) system where F = F(u), H = H(ϕ) and the internal states C v and H r correspond to solutions of the discrete evolution equations in Section 3.1.1 and Section 3.1.2, respectively. Accordingly, in (56), C v = C v (F(u)) and H r = H r (H(ϕ)) are functions of the external fields u and ϕ. History states at time t = t n , while not explicitly mentioned, are of course relevant to the solution procedures for C v and H r . Also, the auxiliary internal states B r and ∆Λ are not mentioned in the global problem above as they do not appear explicitly in the energy densities W i but only in the local evolution equations.
For a complete overview of the solution of the time-discrete IBVP, we provide Algorithm 2 that summarizes the global scheme.
Thereby, we remark that the iterative solution of the time-discrete system (56) is conducted in individual steps for the solution of the external and internal states (lines 9 to 11). In particular, we emphasize that the solution increments, which eventually yield the updated external state {u, ϕ} k+1 n+1 in line 9, are obtained from the linearization of (56). The internal states are updated in separate steps afterwards by a solution of the respective nonlinear evolution problems. For a fast convergence of the loop started in line 8, the aforementioned linearization has to be consistent in the sense that changes in the internal states due to changes in the external states must be accounted for. This is known as algorithmically consistent linearization and is outlined in Appendix A.
Algorithm 2: Generic overall solution procedure for time-discrete problem Data: ∆t, t start , t end , u start , ϕ start , C v start , p v start , H r start , B r start , ∆Λ start 1 n ← 0; 2 t n ← t start ; 3 {u n , ϕ n , C v n , p v n , H r n , B r n , ∆Λ n } ← {u start , ϕ start , C v start , p v start , H r start , B r start , ∆Λ start }; 4 while t n < t end do 5 t n+1 ← t n + ∆t;

Spatial discretization
The initial boundary value problems considered within the present work can be discretized by classical finite element approaches for magnetomechanics. In the numerical examples we restrict ourselves to simplicial meshes and employ second-order Lagrange elements for the displacements u and also for the scalar-valued magnetic potential ϕ [22,30]. In general, we employ a quadrature of degree two but under-integrate (degree zero) the volumetric (dilatation) terms in the stored energy W m . The discrete values representing the internal states will be associated with finite element quadrature points of the degree-two quadrature rule. The actual implementation is done in the Python front-end of the finite element toolbox FEniCS [56] named dolfin together with materiauXdolfin from the Materiaux material modeling kit for the convenient handling of evolution equations and internal states. In addition, we have also developed Abaqus user elements (UELs). Both the FEniCS and the Abaqus implementation are provided as supplementary material. Furthermore, code examples showing the complete chain from the creation of a material model up to the basic definition of a discrete nonlinear FE problem are provided as supplementary documents hyperelastic model.pdf, viscoelastic model.pdf and ferromagnetic model.pdf.

Deformation in free space
In the Lagrangian finite element formulation of magnetomechanics one requires an appropriate extension of the deformation map y(X) from the actual bodies V into the surrounding empty space V a , whereby "appropriate" essentially means det F > 0. Since empty space has practically vanishing elastic properties, one can assign any reasonable soft elastic material law to the empty space to guide the adaption of the free space to the deformation and motion of the embedded magnetoelastic, deformable bodies. A robust and clean approach to this problem is the one employed by Psarra et al. [32], Psarra et al. [33] and Mukherjee et al. [42]. This method essentially uses (non-local) multipoint constraints that control the displacements in empty space based on the displacements at the boundary of the body ∂V.
While this approach is supported by the FE software Abaqus via the *EQUATION command if the constraints are linear or by means of either a non-local penalty method or *MPC command if the constraints are nonlinear, it is not directly supported by the current stable release (2019.1.0) of the FE toolbox FEniCS. Instead, with FEniCS, one may follow the approach of Pelteret et al. [63], which removes all mechanical stiffness from the air domain V a and fixes the discrete deformation within V a . In turn, a suitable extension of the deformation in V to V a is computed in a completely separate step. Our approach of such a staggered scheme is detailed in Section 3.3.2. In comparison with [63], our version is designed for larger displacement increments and thus allows for larger load steps, which significantly accelerates the overall solution procedure in most scenarios. Both the non-local multipoint constraints and the staggered scheme are by construction free of any artificial stiffness contribution in empty space and thus regarded as accurate. There are, however, important differences between the two approaches. On the one hand, the staggered scheme fails in certain problems, such as large transformations induced in beam bending problems with strong coupling between magnetic field and magnetoelastic structures (e.g. in the beam example presented in Section 4), a situation that in general poses no problem for the non-local multipoint constraint approach. In turn, the latter may be not trivial to implement for complex geometries (such as the one considered in [58]), whereas the staggered approach does not suffer from such a limitation.
A third possibility -used extremely often in the literature -is to directly equip empty space with a fictitious auxiliary "material" law, i.e. treat empty space as an extremely soft non-magnetic elastic medium [22]. This, however, introduces an additional stiffness to the original problem, which can be significant depending on the relative stiffness of the MRE structure. In fact, the artificial stiffness is limited from below due to ill-conditioning of the discrete system of equations and, depending on absolute stiffness values, also due to spurious deformations in regions with strongly non-uniform fields [86,Section 9.3]. Thus, we rather regard this approach as a fall-back solution, albeit very useful in complex geometrical problems to get a qualitative estimate of the process.
As a further possibility, we mention the work of Liu et al. [87] who solve in an alternate manner the decoupled mechanical and magnetostatic subproblem until convergence in a staggered manner is achieved. By that, they apparently avoid the spurious magnetomechanical coupling in free space and thus are able to use a much softer elastic material in the empty domain. In turn, they have to pay the cost that comes along with important coupling effects. It eventually depends on the problem, whether such an approach is advantageous or not.
Below we present essential details of three important ways to extend the displacement field into the free space surrounding the deformable magnetomechanical bodies, namely (i) non-local multipoint constraints, (ii) the staggered scheme with a separate deformation phase for the free space, and (iii) the treatment of empty space as a very soft non-magnetic solid.
In addition, in 4.1 we assess the accuracy of approach (iii) through a comparison with approach (i). In this case of large beam deflections and thus large air deformations the staggered scheme becomes highly unstable and is not shown. Instead, in Section 5, the staggered scheme is used successfully to simulate a rotating 3D MRE ellipse subjected to non-monotonic and rotating magnetic fields. 11

Non-local mesh-coupling
In order to obtain a numerically feasible "deformation" in the free space, we apply a linear constraint on every node in the air domain V a following the approach described first in [33]. This implies that all air nodes move relative to the deformable body and any spurious deformation modes induced by the magnetic non-uniformities become irrelevant since they are canceled out by the imposed deformation. To achieve this, we first define a distance coefficient for each node in V a . This requires to find the closest node on the boundary of the deformable body ∂V to the node in V a . This is easily achieved by a simple nearest neighbor search algorithm. The distance coefficient is subsequently defined to be where L Air is the side length of the free-space box 12 and n is the number of active constraints and by definition equals the number of mechanical degrees of freedom in V a , i.e., the subscript i = 1, 2 in two-dimensions and i = 1, 2, 3 in three-dimensions. 11 We note here that the staggered and the non-local constraint scheme have shown the same accuracy in the two-dimensional beam problem but considering a stiffer beam material to reduce the overall deformations. 12 More general, L Air can (or has to) be replaced by any value large enough to ensure that d Subsequently, we employ a constraint between every pair of nodes, the one at the boundary of the deformable solid and that in the air, which reads This constraint is a linear constraint between the displacement degrees of freedom of the body boundary ∂V and those of the free space V a . Obviously a more complex non-linear constraint may be imposed if necessary. The linear constraint may be imposed using various different approaches. In the earlier works of Psarra et al. [32] and [33], a penalty method has been implemented. In this work, given that the constraints are linear, we employ a direct elimination technique by use of the pre-built *EQUATION command in Abaqus. This way the constraint is exactly satisfied without the use of any penalty parameter. Various checks that have been done show that the elimination technique and the penalty method proposed in the earlier studies delivers practically the same results. More importantly, the proposed method allows for large transformations of the air without compromising the convergence rates and accuracy of the magnetomechanical solutions.

Staggered approach
In this approach, the discrete deformation in empty space V a is fixed during the numerical solution of the coupled boundary value problem and only adapted to the deformation on ∂V in a completely separate step [63]. As a consequence, the solution domain of u in Line 9 of Algorithm 2 is not V ∞ but actually V whereas ϕ is still solved for in V ∞ . The adaptation of y or, respectively, u in V a then requires a correction of the magnetic potential ϕ in V a , possibly followed by re-computation of the magnetomechanical equilibrium. The simplest implementation of this approach first solves the nonlinear magnetomechanical problem and then, if required, computes a new freespace deformation. This comes, however, with the caveat that since the boundaries of physical bodies move but the neighboring ones in empty space do not, the boundary deformation is limited. This results in certain constraints on the size of load steps or time increments before the deformation of free space is adapted.
This problem even increases with increased spatial resolution of the problem geometry, i.e. finer meshes. As a remedy, the free-space adaption step can be performed before the assembly of the residual vector and the system matrix of the coupled problem such that free-space deformation happens within the nonlinear solver loop. The abstract free-space deformation problem to be solved before 13 Line 10 in Algorithm 2 is given as where u| V a and u| V indicate a restriction to the disjoint domains V a and V, respectively. Thus, the solution space Note that it is not required to solve this problem to a high accuracy since its sole purpose is to prevent material interpenetration, i.e. to maintain det F > 0. The next step after the free-space motion is to correct or restore the magnetic equilibrium, which can be approximately done by a projection that minimizes the L 2 -distance between the current h-field before and after free space motion. In the case of strong structural coupling, that is when the shape of the magnetoelastic body or structure is highly sensitive to magnetic fields, even tiny changes in the current configuration of the free space affect the deformation of the body. Then it is paramount to repeatedly equilibrate the coupled problem and (re-)adapt the deformation in free space before proceeding to the next load step. Otherwise, both the accuracy and the robustness of the staggered scheme can suffer dramatically.
In principle, the choice of the actual material model is quite arbitrary. In the present work we employ an energy density of the form Ψ mm (I 1 , I 3 , X) = w(X) with G 1 = G 2 = 1, G ′ = 0.5 and a weight function where d el (X) is the "diameter" (in FEniCS/UFL parlance) of the undeformed finite element containing X. A variant of this approach consists in applying w only to G i but not G ′ .

Remark 7.
We note that the additional computational effort to solve the free-space deformation problem and the subsequent correction of ϕ is usually negligible since the linear systems involved are much smaller and much easier to solve than those for the coupled problem. Moreover, free-space deformation is not always necessary but may be triggered by some heuristic criterion. In the present work, we initiate free-space motion as soon as the maximum norm of the discrete residual of the free-space motion problem (59) reaches a certain threshold or in case of material interpenetration.

Auxiliary energy density of free space
The ad-hoc auxiliary energy density given below is proposed in order to minimize the resulting perturbation of the original problem, i.e. minimize the additional stiffness. For this purpose we employ a neo-Hookean-type energy of the form (61) where we let the weighting function w(X) depend on the volume V el 0 (X) of the undeformed finite element containing X, i.e.
such that the stiffness decreases with element size, whereby V ref is a certain reference element size and w min and w max refer to the lower and the upper bound of w(X). This heuristic proposition is appropriate for meshes that are finer in the vicinity of magnetic bodies and gradually become coarse at larger distances from them. We note that finite element discretizations may suffer from spurious magnetomechanical deformations in very soft non-magnetic media as described in [86, Part II, Section 9.3.1]. Thus, an auxiliary energy density with varying coefficients as the one above is only an attempt to balance the required stiffness to cope with spurious magnetomechanical effects with a minimal perturbation of the original problem. It thus has certain limits, in particular, for boundary value problems of mechanically very soft magnetic media and compliant magnetic structures. As we will see, this approach may lead to quantitative errors as demonstrated by a comparison with the non-local multipoint constraint approach and thus must be used with great care. In turn, it is easier to employ, especially in domains where the free-space geometry is extremely complex (see for instance the modeling of an MRE device by Dorn et al. [58]).

2D numerical examples: a two-dimensional magnetoelastic beam
In this section and the following one, we present a set of examples which demonstrate the capabilities of the proposed numerical framework and its implementation. At the same time, the examples demonstrate the combined effect of viscoelasticity and dissipative ferromagnetism as a motivation for future material modeling efforts. We note that all material models employed are implemented for the general three-dimensional case being totally agnostic of the spatial dimension of the initial boundary value problems below.
In all the simulations shown in the remainder of this section, we consider a plane strain setting and employ the same set of material parameters, whereby we consider four material combinations as shown in Table 1. The  comparison of these cases provides important insight into the practical effect of the dissipative behavior of both the magnetic particles and the elastomers. The precise forms of the individual material laws can be found in Section 2.3. The non-dissipative (hyperelastic) variant for the elastomer is simply the hyperelastic part Ψ eq of the behavior, whereas the non-dissipative magnetic response is given by (49). The actual parameters employed are provided for each example below. The first example considers a chain of circular magnetic particles embedded in a slender beam. This geometry is inspired by Zhao et al. [45] and Mukherjee et al. [42]. However, in contrast to these earlier studies, in the present contribution, we numerically resolve a small, finite number of particles instead of considering a (macroscopically) homogeneous beam mainly due to lack of a complete ferromagneto-viscoelastic continuum material model. The structure under consideration is depicted in Fig. 5. The center coordinates of the thirteen particles are chosen arbitrarily but are fairly uniformly spaced (see Table 2).
As indicated in Fig. 5, those instances of the beam that carry magnetically hard particles are either premagnetized in x 1 -or x 2 -direction with a field of b ∞ = b ∞,pre . During "actuation" of the (premagnetized if magnetically hard) beam, the field b ∞ = b ∞,op is applied in the perpendicular direction such that the angle between premagnetization and actuation field is 90 • in either case. The premagnetization field reaches a peak value of b ∞ = 2 T before linear decrease to zero. During this phase, we block all displacements at the boundary of the beam ∂V, at the Table 3 Material parameters for the particle-chain example.

Domain
Function Parameters (19) magn: µ 0 = 4π × 10 −1 µm T A −1 outer boundary of the domain ∂V ∞ and within the "empty" space V ∞ \ ∂V. For magnetically soft beams, the premagnetization phase is irrelevant and thus is omitted.
In the subsequent actuation phase, the field magnitude b ∞ = b ∞,op is linearly increased starting from zero at a rate ofḃ ∞ = 5 Ts −1 . In this phase, the beam is only clamped at its left end as indicated in Fig. 5a. The maximum operation load is chosen such that the beam shows a normalized deflection of |u 2 /l| ≈ 0.5 whereby u 2 is evaluated at point "A", the right-top corner of the beam (Fig. 5a). Limiting the deflection prevents crashes of the simulation due to excessive mesh distortion in the empty space. Because of substantially different responses for magnetically hard and soft cases, the above described approach results in actuation fields of b ∞ ≤ 0.06 T and b ∞ ≤ 0.8 T, respectively. For both the magnetically hard and soft cases, the loading stage is followed by a relaxation phase with a duration of 1 s before unloading. For the cases under consideration, this was sufficient for a practically complete relaxation. The unloading rate is of the same magnitude as the loading rate, i.e. |ḃ ∞ | = 5 T s −1 . In the present example, we employ an elastomer phase whose stiffness is somewhat stiffer than that considered by Zhao et al. [45] and Mukherjee et al. [42] but at the same order. Table 3 reports the material parameters for the individual phases.

Effect of the modeling of free space
In order to assess the error introduced by the auxiliary air stiffness, we compare the deflection responses in the hard case with premagnetization in x 1 -and operation field in x 2 -direction. For this, we use as a reference accurate solution the non-local (NL) constraint scheme (see Section 3.3). Note that for the aforementioned choice of material parameters, the staggered scheme is rather unstable and is not shown. For the non-uniform (NU) air stiffness (61), we employ G 1 = G 2 = 40 Pa in free space, which is less than one percent of the elastomer shear moduli in the beam. For the weight function (63), we employ parameters w min = 1 × 10 −5 , w max = 1 and V ref = 1 × 10 −4 . We emphasize that these choices are heuristic and mesh dependent. In addition, a uniform (U) weight function scheme is also shown in an attempt to reveal the relative stiffness induced by the presence of a soft air model. Fig. 6 shows that the additional stiffness (see (61)) of the free space (air) domain does not cause any significant error provided that it is employed with an appropriate weight function. Instead, the remaining curves that use the trivial weight w(X) = 1, i.e. uniform air stiffness, lead to significant deviations from the NL zero air-stiffness model. We note again that due to involved numerical issues, it is not possible to simply further reduce the uniform air stiffness. Fig. 7 exhibits the employed non-uniform weight function. It clearly displays a unit weight only inside a boundary layer around the beam. Outside this layer, where the non-uniform magnetic field contributions are expected to be much smaller than in the vicinity of the beam, the stiffness is reduced drastically. The beam itself is not part of the air domain and thus not shown.

Hard versus soft magnetic particles
In the following, we investigate the effect of soft versus hard magnetic particles upon the beam deflection using the NU stiffness implementation with the aforementioned parameters, which turned out to be considerably faster than our Abaqus implementation. At this point, it is useful to note again that the staggered approach to free-space deformation (see Section 3.3.2) did not work in the example under consideration. The approach suffered from the extreme sensitivity of the beam's deformation with respect to the magnetic solution, which is slightly perturbed by each adaption of deformation in free space. Our staggered implementation has not been designed with such issues in mind and further research has to be done along these lines.
We begin the discussion of the results with the magnetically hard cases hard and visc/hard with premagnetization in x 1 -direction and actuation field in x 2 -direction, as depicted in Fig. 8. We observe a significant effect of the elastomer viscoelasticity, which is not only due to the chosen value of the viscosity η = 40 kPa s but also due the high sensitivity of the beam with respect to the applied field which leads to a rather high effective loading rate. Besides that, the selected contours at instances "1" to "4" show that the magnetization (yellow arrows) follows the center line of the (deflected) beam. This indicates that the deflection of the beam is a result of its bending stiffness/compliance, i.e. a structural property, rather than due to local shearing. Furthermore, we note that the bending process is driven by the magnetic torque that is applied by b ∞ on the magnetized particles and which is of greatest magnitude for an angle |∡(b ∞ , m)| = 90 • . For additional validation and in order to put these results into context with recent works on h-MREs [45,88,89], we compare our fully coupled PDE approach to simpler torque models (presented briefly in Appendix C) in Fig. 9. Since the interaction of particles is rather week for the example under consideration, the difference is mainly explained by the mismatch of the saturation magnetization and the actual magnetization in the full PDE model. The latter, albeit practically uniform within each particle, has a magnitude of only 0.6 MA m −1 (saturation value is 0.67 MA m −1 ) despite the fact that a premagnetization field magnitude of b ∞ = 2 T has been applied. The reason for this is that in actual magnets upon unloading from the fully magnetized state, a part of the magnetization decreases following the magnetic laws presented in the context of Fig. 4b. Therefore, a key weakness of the torque model approach arises from the difficulty that the remanent magnetization is in general not known a priori -even for the rather simple structures considered here. In a sense one needs to solve the premagnetization problem as accurately as possible in order to feed the torque model. As a remedy, one might come up with an intermediate approach that exploits the approximate uniformity of magnetization for the ferromagnetic evolution. The situation of course becomes non-trivial when non-uniform magnetization is to be expected upfront or when particles come closer to each other as a result of large structural deformations as might be the case for magnetoelastic grippers. In such cases, the solution of the coupled PDE system with ferromagnetic evolution and -importantly -also the free space is paramount to obtain meaningful results.
Next, we investigate the converse premagnetization/loading case, i.e. b ∞,pre in x 1 -and b ∞,op in x 2 -direction with the corresponding results shown in Fig. 10, where we only consider the case visc/hard. In this example, the Fig. 10. Deflection response of a viscous elastomer beam comprising magnetically hard particles with vertical premagnetization under horizontal applied field. The contours at instances "1" to "4" refer to the current magnetic field magnitude, the yellow arrows indicate the magnetization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) deflection is downwards which is a direct result of the orientation-switched premagnetization and loading directions in comparison with Fig. 8. Apart from changed directions, the deformation response is essentially the same as before, which confirms the torque-driven nature of the bending problem for the premagnetized beams under consideration. Nevertheless, the magnetic field contours at instances "1" to "4" in Fig. 10 show obvious differences to Fig. 8 due to different (pre)magnetization and actuation field directions.
In addition to the magnetically hard cases above, we show for comparison the response of a magnetically soft beam (case visc/soft) in Fig. 11 with the applied field in x 2 direction. We chose instance "1" to be not at zero applied field but near b ∞ = 0.4 T to better visualize the magnetization in the still practically undeformed beam. As one can see, the magnetization is perfectly aligned with the x 2 direction, which is the direction of the applied field. As the loading increases, we reach instance "2" where the beam is still magnetized almost perfectly in x 2 direction. Thus, the deflection of the beam at hand cannot be driven by magnetic torques exerted upon the particles. Instead, even beyond "1", there is only a rather small deflection observed which results from the imperfect wavy particle alignment. At later stages, approaching "2", sufficiently deflected portions of the beam experience a compass-type torque, which drives the bending further. We observe a kink at point "2" where the magnetic field has reached its maximum value of b ∞ = 0.8 T such that the beam enters the viscous relaxation regime. At instance "3" the relaxation is complete, which can be nicely observed in a plot of the deflection over time instead of applied field, which is deferred to Fig. 13. After the beam has come to rest at "3", the effect of viscosity is again clearly visible during unloading from "3" to "4" and also in the following viscous relaxation at vanishing applied field. We emphasize that the maximum applied (actuation) field in the magnetically soft case has a magnitude of more than ten times of that in the magnetically hard cases shown previously.
The responses for all three loading cases are compared in Fig. 12. Therein, we point out the (almost) "mirrored" responses for the visc/hard(hard) cases depending on the premagnetization and actuation field directions as well as the instability that governs the response of the magnetically soft beam during loading. In addition, we show for completeness a green line that corresponds to cases that show no deflection and only minor axial deformation (not shown explicitly).
We close this section by investigating the deflection response in terms of |u 2 |/l as a function of time in Fig. 13. The loading and unloading responses of the beam with magnetically hard particles are shown in Fig. 13a, where one observes the immediate bending reaction (solid) to changes in the external field (dashed). The inset shows that the kink in the deflection graph happens at the same time as the loading stops to increase, i.e. when the b ∞ graph continues horizontally. From this point on, the change in the deflection is entirely due to viscous relaxation.
In Fig. 13b, we show the response of the beam with magnetically soft particles during loading and unloading. Therein, one observes a delayed reaction during loading, while such a delay cannot be observed during unloading. This is because the deflection initially results from the imperfect particle alignment. Later, once a certain deflection has been reached, compass-type torques exerted on the deflected portions of the beam further drive the deformation. Fig. 11. Deflection response of a viscous elastomer beam comprising magnetically soft particles under vertical applied field. The contours at instances "1" to "4" refer to the current magnetic field magnitude, the yellow arrows indicate the magnetization. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 12. Comparison of deformation responses of viscous elastomer beams comprising magnetically hard or soft particles under different magnetic loading directions. The "red", "marine blue" and "violet" responses are detailed in Figs. 8 to 11. The labels "S1" and "S2" correspond to premagnetization and operation steps, respectively. The respective arrows indicate the direction of the applied field during these steps. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

3D Numerical examples: An ellipsoidal particle in a viscoelastic matrix
As a second geometry, we consider a single ellipsoidal magnetic particle embedded in an elastomer cube of side length l = 1, which itself is surrounded by an empty space box of side length L = 6 as depicted in Fig. 14. The angle between the major axis of the ellipsoid and the 1-3-plane is π/4. The major and minor axes of the ellipsoid have lengths 0.6 and 0.3, respectively. Furthermore, θ denotes the angle of the applied field b ∞ defined at the x 1 -x 3 -plane. Note that the x 3 -component b ∞ 3 = 0 throughout this section. The material parameters are summarized in Table 4, where we highlight that the matrix material is chosen to be rather soft with a "total" shear modulus less than 100 kPa. The arrows next to the labels "S1" (premagnetization) and "S2" (operation) indicate the direction of applied magnetic field in the respective phase. The time axis is broken to exclude excessive relaxation periods. Table 4 Material parameters for the embedded-ellipsoid example.

Domain
Function Parameters We consider two representative cases, (i) non-monotonic uniaxial magnetic loading in Section 5.1 and (ii) a nonmonotonic, rotating magnetic loading in Section 5.2. In both cases, the adaptation of the free space deformation to that of the boundary of the body is achieved by the staggered scheme, as described in Section 3.3.2.
The evaluation of useful quantities associated with the ellipsoid, such as its rotation and magnetization is performed near the origin "O" (0, 0, 0). Note that the ellipsoid does not deform but only rotates such that all Eulerian magnetic quantities may be considered approximately uniform in the interior of the ellipsoid under a uniform applied magnetic field b ∞ (but see discussion about this point in [19]).

Non-monotonic uniaxial applied magnetic field
The magnetic loading in this example is along the x 2 -direction only. The magnitude of the applied field b ∞ is linearly increased over time up to a maximum value of 2 T. Then, it is held for a relaxation time of 5 s before it is decreased to zero again as detailed in Fig. 15. In order to visualize the effect of viscosity, i.e. the rate-dependence of the responses, we consider two different loading rates,ḃ ∞ = 0.5 Ts −1 ("slow"; see Fig. 15a) andḃ ∞ = 5 Ts −1 ("fast"; see Fig. 15b).
We first discuss the response of magnetically hard ellipsoids. Fig. 16a shows contours at six different instances for a non-viscous matrix material (case hard) and "fast" loading. Interestingly, the magnetization is almost perfectly aligned with the applied field (vertical direction) throughout the loading phase while the ellipsoid rotates slightly. Only from instance "5" to "6", we observe a slight rotation of the magnetization vector. This rotation occurs during unloading, where the magnetized ellipsoid rotates back to its initial position. Figs. 16b,c provide details of the rotation evolution of the ellipsoid and the magnetization vector. In both subfigures, one may observe the effect of matrix viscosity and loading rate, whereby the "slow" loading rate is already quite close to the non-dissipative limit. Obviously, rotation of the ellipsoid is rather small and the alignment of the ellipsoid does not have a strong effect on the direction of the magnetization evolution. While this result might be surprising, it is perfectly in line with experimental observations [90]. As can be seen from Fig. 16d, viscosity does not affect the typical hysteretic magnetization response in terms of the magnitude m, which is a rather important result for the understanding of coupling between magnetic and mechanical properties in such materials.
Next, we discuss the results for magnetically soft ellipsoids under the same loading conditions. Fig. 17a shows contours at six different instances for viscous matrix material (case visc/soft) and "fast" loading. The magnetization is almost perfectly aligned with the applied field direction. This time, the ellipsoids show a much more important rotation such that the misalignment between the major axis and the magnetization is much smaller than in the magnetically hard case (c.f. Fig. 16). As for the magnetically hard case discussed previously, the effect of viscosity on the ellipsoids rotation is rather pronounced. In turn, in this example, we observe no significant effect on the magnetization magnitude as well as on the resulting magnetization direction. In this example of a soft magnetic particle, one thus may safely consider a negligible coupling between magnetic response and viscoelasticity of the matrix.

Rotating magnetic field
Let us now consider the case where the ellipsoid is first pre-magnetized by a procedure similar to the one discussed in Section 5.1, whereby the field magnitude is not decreased to zero but to b ∞ = 0.5 T. Once this value is reached, the field b ∞ is held for five additional seconds for relaxation. Then, the field angle θ is changed from 90 • to 270 • , while the magnitude b ∞ remains unchanged as depicted in Fig. 18. The loading rates for the magnetic field magnitude and the field angle areḃ ∞ = 5 Ts −1 andθ = 90 • s −1 , respectively. A first "premagnetization" phase ("S1") with a maximum field of b ∞ = ∥b ∞ ∥ = 2 T is followed by an "operation" phase ("S2") with a field of 0.5 T with changing angle θ (see Fig. 14b) for t ≥ 1.7 s. The relaxation periods where b ∞ and θ are held constant have a duration of 0.5 s.
We focus here only on viscoelastic matrix materials such that the two material combinations under consideration are visc/hard and visc/soft (see Table 1).
The main features of the responses of both cases are depicted in Fig. 19, where part (a) shows deformed configurations for the case visc/hard and part (f) shows deformed configurations for the case visc/soft. In both subfigures the contours correspond to the current magnetic field b, yellow arrows indicate the magnetization m, blue arrows the applied field b ∞ and thin white arrows the orientation of the major axis of the ellipsoid. In Fig. 19a, we see that the ellipsoid undergoes a significant rotation from around 40 • to almost 90 • with respect to the vertical direction as the applied field rotates. During this process, the angle between the major axis and the magnetization of the ellipsoid changes only weakly. Instead, the angle between the magnetization direction and the applied magnetic field increases rapidly once the applied field rotation starts. In a sense the hard magnetization of the ellipsoid resists to the applied external field b ∞ . The peak rotation of the magnetically hard ellipsoid is not observed at the maximum rotation of the applied field attained at 180 • but some time before -between instances "4" and "5" -where the angle between the applied field and the magnetization is ∡b ∞ − ∡m = 90 • . Near this instant, the magnetic torque reaches its maximum. Beyond that angle, the magnetic torque decreases and the mechanical response of the matrix pushes back the ellipsoid by some small angle. Note that the magnitude of the magnetization also changes significantly during the field rotation phase "S2". In particular, in Fig. 19e, one can see a kink in the yellow graphs. This indicates the onset of a magnetic evolution: beginning from a field angle θ ≈ 115 • the magnitude of the magnetic internal variable decreases and its angle varies non-monotonically. The respective details can also be observed in Figs. 19b to 19e.
The behavior of the magnetically soft ellipsoid is rather different as can be seen from the deformed configurations and contours in Fig. 19f but also from the graphs in Figs. 19b to 19e. In this case, the magnetization rotates much more since no energy is dissipated during the rotation of the magnetization with respect to the ellipsoid. As a result, the ellipsoid initially follows closely the rotating applied field and magnetization but eventually loses track of the magnetization near instance "3", albeit only weakly. At this point, the magnetic torque on the ellipsoid changes direction. As the applied field reaches its final orientation at instance "6", the magnetization is reversed with respect to instance "1" but the mechanical state of the system at "6" is practically identical to "1". This behavior of the magnetically soft case is explained by the fact that the magnetic torque exerted on the ellipsoid is due to the shape anisotropy of the ellipsoid. 14 Hence, the system is invariant with respect to a simultaneous reversal of applied field and magnetization. What essentially distinguishes "1" from "6" is only a small effect of viscosity near these two states.
Finally, we point out that in the magnetically hard case, the evolution of the magnetization is a manifestation of changes in the internal magnetic state H r as shown in Fig. 20. Therein, we observe changes in both the magnitude H r and the angle of H r with respect to the vertical direction.

Conclusion
The capabilities and the generality of the computational framework and its implementation in two different software packages, FEniCS and Abaqus are demonstrated by carefully selected examples inspired by prospective locomotion and actuation applications. Specifically, the effect of schemes used to model the surrounding air are also discussed in some detail and their advantages and disadvantages are clearly identified. In particular, we find that the non-local constraint scheme is able to provide the more accurate responses, yet it may require a more case-to-case analysis to implement in complex confined geometries (such as air spaces between a magnet and an MRE or an internal material void space). The staggered scheme, where the air and solid are solved separately, is able to simulate accurately more complex geometries but tends to become unstable and non-convergent in certain cases of strong structural magnetomechanical coupling (such as in the beam deflection problem). An ad-hoc scheme using a non-uniform air stiffness distribution is able to deal with all cases, but it is geometry dependent and mesh dependent and requires an independent error estimation (for instance by using one of the two previously discussed schemes). Finally, the classical approach of a weakly stiff air, which is the easiest of all to employ, is shown to lead to quantitative inaccuracies. Yet, it is a powerful tool that allows for a proper qualitative analysis of the system. Finally, a simple torque model for hard MREs is also used, where it is shown to have a good qualitative agreement but in the case of realistic hard magnetic responses it can lead to quantitative errors since the actual magnetization value is unknown in the magnetic phases. This last approach also requires an independent calculation of the premagnetization response by some other numerical or analytical means and is only of use in very small magnetic actuation fields, where the response remains magnetically non-dissipative.
In the first set of examples, we discuss the magnetomechanical response of a slender viscous elastomer beam carrying circular particles that are arranged in an imperfect wavy chain pattern. For magnetically hard particles, we consider different directions of premagnetization combined with a magnetic loading transverse to the premagnetization direction. For these cases, we observe that the magnetic torque exerted by the applied field upon the magnetized particles governs the deflection of the beam. In contrast, such local magnetic torques did not play any role for the beam with magnetically soft particles. In both cases, we find significant viscous effects despite the rather low viscosity parameter and the small underlying local strains in the beams. This is directly related to the fast loading rates applied, which are of potentially great practical relevance in rapid actuation scenarios. In a sense, the obtained results imply the necessity of viscoelasticity in the design of such beam actuators.
In the second set of examples, we present simulations of an elastomer with an embedded magnetic ellipsoid exposed to a magnetic field. We study the response and the magnetization process under uniaxial magnetic loading where we observe significant qualitative differences between elastomers with magnetically hard and magnetically soft particles. Moreover, we study extensively the effect of viscosity of the elastomer for two different loading rates upon the magnetization response. This effect is found to be more significant for the magnetically hard case than for the magnetically soft one, especially for the orientation of the magnetization vector and less for its amplitude. This is a somewhat surprising result that demonstrates the complex interactions of mechanically and magnetically-induced dissipative phenomena in h-MREs. In the second loading case, we consider a rotating magnetic field whereby the magnetically hard ellipsoid is premagnetized beforehand. Then, both the magnetically hard and soft ellipsoids initially follow the rotation of the applied field. Nevertheless, the magnetically soft ellipsoid snaps back to its initial position after some critical value is reached due to the symmetries of the geometrical ellipsoidal shape. By contrast, the (premagnetized) magnetically hard ellipsoid continues to follow the applied field much further until the peak torque is reached.
From our numerical results, one deduces that both the magnetization history in the case of magnetically hard particles and the viscosity in the case of sufficiently fast magnetic loading play a significant role for the overall response of MREs. This observation does not only hold for device-type geometries (e.g., the beam problem) but also for more material-type systems (e.g., the ellipsoid embedded in a matrix). Moreover, both dissipative phenomena interact in a non-trivial way such that carefully designed numerical simulations are paramount for a deeper understanding and accurate modeling of both soft and hard MREs.
By virtue of the results presented here, the proposed computational framework can be easily used to carry out numerical homogenization studies in the context of microstructurally-guided modeling of magnetically soft and hard magnetorheological elastomers similar to the works of Mukherjee et al. [27] and Mukherjee et al. [42], respectively. Such a work is left for a future study.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. as a function of the external state, i.e. I n+1 = I n+1 (u n+1 , ϕ n+1 ) which has implications on the linearization of the problem with respect to {u n+1 , ϕ n+1 }. For what follows, we introduce the notation From the latter we obtain which can be substituted for ∂I/∂S| k n+1 · ∆S k n+1 in (A.3). The resulting linearization is symmetric if and only if which depends on the chosen parameterization of both W and D, which governs E. The above symmetry condition is necessary and sufficient for the existence of an incremental (reduced) potential Π that governs the evolution of the internal state I as presented in Sections 3.1.1 and 3.1.2.

Appendix B. Viscous evolution equation
In this section, we carry out the algebra leading to the viscous evolution equation for a compressible (30) and incompressible (31) viscous dissipation potential. We start by considering the non-equilibrium energy function (27) and the viscous dissipation potential in (28). Before proceeding further, we first write down some intermediate results necessary for the subsequent algebra. Those read where use of relations (24) has been made.
To proceed further, we rewrite the dissipation potential as where J i jkl = δ i j δ kl /3 is the hydrostatic fourth-order projection tensor. In the following, we evaluate the two terms involved in the GSM equation (29). Using index notation and C v as the state variable, we obtain the intermediate result To get this expression, we have used the intermediate result ∂d v r s /∂Ċ v i j = √ C v ri √ C v js . Now using the partial results of (B.1), the first term in (29) is readily written as Gathering the last two results provides the GSM equation Multiplying next this last equation from the left side with C v and taking the trace, such that trĊ v C v-1 =Ċ v · C v-1 , one gets Contracting again Eq. (B.5) from both sides with C v leads to Substituting the intermediate result (B.6) into this last equation, leads to After some straightforward algebra, this last relation may be recast to the compressible form given in (30). The incompressible form is then readily obtained in Eq. (31).

Appendix C. Comparison with a simplistic torque model
In some recent publications concerned with h-MREs [45,47,88,89], the magnetization of the structure under consideration has been assumed to be known a priori and more importantly be piecewise (per domain) constant. Also, the modeling is done for small magnetic fields such that the amplitude of the magnetization remains "almost" constant thus inducing no dissipative effects. The range of validity of such an approach obviously needs to be assessed to avoid entering in actual hysteretic regimes of response. However, that is only possible via the full framework including dissipation.
Even so, the assumptions above greatly simplify the coupled PDE problem which becomes purely energetic in magnetism. In fact, it reduces the magnetomechanical problem to a purely mechanical one with body torque magnetic loads, whereas the air domain V a is no more needed.
We shall compare this approximate approach with our solutions of the full problem. For that purpose, one may simply assume, even though it might be inconsistent with the Maxwell field equations, the magnetization in the particles to be uniform and equal to 15 m = RM. Therein, R is the rotational component of F and M is the given and fixed magnetization density of a rigid particle. Specifically, M = (m s , 0, 0) for a particle premagnetized in X 1 -direction. The body torque then automatically results from the energy contribution We refer to Mukherjee et al. [42] for a discussion of torques in the more complex case of deformable magnetic bodies. Below we summarize the energy densities employed in the torque-model BVP. As no air domain and no magnetic solution variable are needed, the energy densities of the elastomer matrix (W m , see (25)) and the quasi-rigid magnetically hard particles (W p , see (37)) reduce to W m = ρ 0 Ψ eq (I 1 , J ) + ρ 0 Ψ neq (I e 1 , J e ) (C.2) W p = ρ 0 Ψ mech (I 1 , J ) + Ψ torque (F). (C.3) As this torque approach does not take care of ferromagnetic evolution, only viscous evolution (30) needs to be considered. The corresponding results are shown in Fig. 9, where both the full and the torque-model employ the same value for the magnetic saturation parameter m s = 0.67 MAm −1 . Note that the choice of using the saturation magnetization value in the torque approach appeared most plausible to us. Other choices are possible but require more involved analysis which undermines the simplicity of the magnetic torque approach.
Remark 8. The torque contribution Ψ torque (F) is objective despite the explicit dependence on F through R because not only R but also the prescribed magnetization vector M transforms in a change of material frame.
Remark 9. The rotational component R can simply be approximated by F in the case of quasi-rigid particles, which greatly facilitates computer implementations.