# A continuum magneto-viscoelastic model for isotropic soft magnetorheological elastomers: experiments, theory and numerical implementation

Kostas Danas, Masami Nakano, Gaël Sebald

► To cite this version:

Kostas Danas, Masami Nakano, Gaël Sebald. A continuum magneto-viscoelastic model for isotropic soft magnetorheological elastomers: experiments, theory and numerical implementation. Mechanics of Materials, 2025, 200, pp.105187. 10.1016/j.mechmat.2024.105187. hal-04781789

HAL Id: hal-04781789

https://hal.science/hal-04781789v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Editor Invited Article

# A continuum magneto-viscoelastic model for isotropic soft magnetorheological elastomers: experiments, theory and numerical implementation

Kostas Danas<sup>a,b,*</sup>, Masami Nakano<sup>c</sup>, Gaël Sebald<sup>d</sup>

ARTICLE INFO

Keywords:

Viscoelasticity

Magneto-mechanics

Magnetorheological elastomer (MRE)

Coupled Finite Elements

Experiments

Finite strains

Constitutive modeling

ABSTRACT

This study deals with the experimental, theoretical and numerical investigation of the nonlinear viscoelastic response of magnetically soft magnetorheological elastomers (commonly known as s-MREs and denoted here simply as MREs) subjected to combined magnetic and simple shear loads. We consider a fairly soft mechanically MRE. The experiments show a strong effect of the magnetic field on the resulting viscosity and hence dissipated energy expanded by the material during a simple shear cycle. Moreover, the effect of frequency on the response is weak indicating strongly nonlinear viscous effects similar to non-Newtonian fluids. An analytical magneto-viscoelastic model is proposed exhibiting magneto-mechanical coupling at both equilibrium and non-equilibrium energies as well as on the dissipation potential. The model is calibrated by solving in a semi-analytical way a simplified boundary value problem (BVP) of an infinite thin MRE strip embedded in an infinite air domain. These simplified solutions are cross-validated by full-field finite element simulations of the experimental setup showing very good agreement between the experimental data and model estimates. This illustrates the validity of the simplified material model for the proposed experimental setup and sets the ground for a more universal experimental protocol to characterize properly the finite strain response of MREs more generally.

1. Introduction

The present work deals with the modeling and experiments of the large strain viscoelastic response of isotropic, magnetically soft magnetorheological elastomers (s-MREs) subjected to combined magnetic and simple shear loads. s-MREs or simply MREs here more generally have gained a lot of thrust the last years owing to their interesting coupled response at finite strains (Ginder et al., 1999; Diguet et al., 2010; Danas et al., 2012) as well as their fine control of surface patterning via instabilities (Psarra et al., 2017, 2019) or actuation in biomedical devices (Dorn et al., 2021; Moreno-Mateos et al., 2022). More recently, novel MREs using permanently magnetizable particles (such as NdFeB) have been proposed (Stepanov et al., 2014, 2017) to trigger higher torques at smaller applied fields (Kim et al., 2018; Mukherjee et al., 2021; Mukherjee and Danas, 2022; Yan et al., 2023; Stewart and Anand, 2023).

The literature in these materials is growing with an extremely fast pace and it is beyond the goal of this work to carry out such a review

but it is important to note that there are still numerous aspects of these materials that are yet to be understood or modeled in connection with experiments. One such aspect pertains to their viscoelastic response at finite strains and large magnetic fields. Finite strain viscoelasticity in the absence of magnetic fields is a rather old subject in mechanics with seminal works found already in the 50–60’s (Green and Rivlin, 1957; Pipkin and Rogers, 1968) and much later using an internal variable formalism (such as the one that will be used in the present study) (Le Tallec et al., 1993; Bergström and Boyce, 1998; Reese and Govindjee, 1998; Kumar and Lopez-Pamies, 2016). In the context of composite materials and homogenization, one may refer to the works of Francfort and Suquet (1986), Lahellec and Suquet (2007), Diani et al. (2013), Chatzigeorgiou et al. (2016) and Idiart et al. (2020) for small strain viscoelasticity and Ghosh et al. (2021) for large strain viscoelasticity.

In turn, very little and much more recent work has been done in the context of magneto-viscoelasticity from the modeling point of

<footnote><sup>a</sup> LMS, CNRS, École Polytechnique, Institut Polytechnique de Paris, Palaiseau, 91128, France</footnote>

<footnote><sup>b</sup> ELyTMAx, CNRS, Tohoku University, Sendai, 980-8577, Japan</footnote>

<footnote><sup>c</sup> SmartTECH Lab. Inc., Material Solutions Center, Tohoku University, 2-1-1, Katahira, Aoba-ku, Sendai, Miyagi, 980-8577, Japan</footnote>

<footnote><sup>d</sup> ELyTMAx, CNRS, INSA Lyon, Centrale Lyon, Université Claude Bernard Lyon 1, Tohoku University, Sendai, 980-8577, Japan</footnote>

<footnote>* Corresponding author at: LMS, CNRS, École Polytechnique, Institut Polytechnique de Paris, Palaiseau, 91128, France.</footnote>

<footnote>E-mail addresses: konstantinos.danas@polytechnique.edu (K. Danas), mnakano@smarttech-lab.com (M. Nakano), gael.sebald@insa-lyon.fr (G. Sebald).</footnote>
view (Jolly et al., 1996; Saxena et al., 2014; Hossain et al., 2015; Garcia-Gonzalez and Hossain, 2020; Haldar, 2021; Rambausek et al., 2022; Lucarini et al., 2022; Gonzalez-Saiz and Garcia-Gonzalez, 2023; Stewart and Anand, 2023). Some of these studies are theoretical or numerical and very few of them attempt to compare with corresponding experiments. In turn, there exist some experiments that study the viscoelastic response of MREs under combined magneto-mechanical loads most of them in the context of s-MREs (Stepanov et al., 2007, 2014; Hiptmair et al., 2015; Sebald et al., 2017; Lucarini et al., 2022; Pierce et al., 2024) and only very few for h-MREs (Yan et al., 2023) (but not extensively).

In most of these studies, the experimental magneto-mechanical setup involves an air gap between the magnetic source and the MRE and as a consequence the fields in the MRE specimen are highly heterogeneous (Lefèvre et al., 2017; Moreno-Mateos et al., 2023) while the response exhibits strong structural effects or shape effects (Diguet et al., 2010). Those two combined make the modeling analytically impossible since the boundary value problem (BVP) at hand has no explicit solution. In turn, use of a full field finite element (FE) simulator coupled with an optimization method that would allow for the calibration of material parameters is currently also out of reach. This is mainly a consequence of the extremely long computation times involved in magneto-mechanical problems that require – especially in this case of complex experimental boundary conditions – the modeling of all parts of the setup as well as a large portion of the surrounding air (Moreno-Mateos et al., 2023). And this is even the case in two-dimensional calculations. It is therefore understandable that still this problem remains up to date largely open.

Motivated by those observations and lack of a robust viscoelastic model for soft MREs in the literature, we propose in the present work a combined experimental-theoretical-numerical analysis that allows to model several features present in viscoelastic soft MREs subjected to large simple shear strains and large magnetic fields. The main difference of the present study lies in the use of the experimental setup of Sebald et al. (2017), which we will show is capable under careful analysis to deliver the actual material response free from boundary and specimen shape effects. The experimental data are then confronted with the semi-analytical solution of a simplified but rigorous BVP allowing us to calibrate efficiently a general MRE material model with viscoelastic effects. In order to cross-validate the assumptions involved in the simplified BVP, we carry out additional full-field FE simulations of the actual experimental setup showing the consistency of the model but also revealing several heterogeneities of the local magneto-mechanical fields that however do not affect the average measured response. This last step – omitted in most of the studies in the literature – allows to close the loop of the modeling approach and lift any possible uncertainties related to magnetic edge and corner effects as well as non-uniformities that may develop due to large strains and specimen shape effects.

The manuscript is organized as follows. In Section 2, we introduce the principal mechanical and magnetic quantities needed for the analysis of the magneto-mechanical problem. In Section 3, we describe the experimental setup together with material fabrication and present a small but representative set of experimental results allowing to reveal the main features of the material response. Section 4 describes the mathematical formulation of the model in the context of finite magneto-viscoelasticity, whereby the energy density and dissipation potential are defined. With this at hand, in Section 5, we propose a simplified BVP that combined with the previous model definitions allows for the semi-analytical incremental solution of the problem under simple shear loads and normal magnetic fields. An optimization process is set in place to allow the calibration of the various material parameters with a selection of the experimental data. Section 6 describes a full-field finite element (FE) BVP simulation of the experimental setup using the previously identified magneto-viscoelastic model thus allowing to probe the assumptions made in the simplified BVP. Finally, we conclude the manuscript discussing potential extensions of the present work and future directions in the context of magneto-viscoelasticity.

## 2. Kinematics and magnetostatics

For simplicity in presenting the main field quantities of the problem, we consider a deformable, magneto-active solid with volume $V_0$ ($V$) in its reference (current) configuration. The magneto-active solid is considered in the present work as a homogeneous material described by coupling energy densities and a dissipation potential to be defined in detail in Section 4. The boundary of the solid is assumed to be smooth and is designated by $\partial V_0$ ($\partial V$), while $\mathcal{N}(\mathbf{n})$ denotes the unit normal on $\partial V_0$ ($\partial V$) in the reference (current) configuration.

In the present formulation, we consider as primary variables

(i) the mechanical displacement field $\mathbf{u}(\mathbf{X})$, which relates the current, $\mathbf{x}$, and reference, $\mathbf{X}$, position vectors via $\mathbf{x} = \mathbf{X} + \mathbf{u}(\mathbf{X})$,

(ii) the magnetic vector potential field $\mathbf{A}(\mathbf{X})$ defined in the reference configuration.

This allows to write the deformation gradient $\mathbf{F}$ and Lagrangian magnetic flux $\mathbf{B}$ as

$$ \mathbf{F} = \mathbf{I} + \mathrm{Gradu}(\mathbf{X}), \quad \mathbf{B}(\mathbf{X}) = \mathrm{CurlA}(\mathbf{X}). \qquad (1) $$

Here, $\mathbf{I}$ is the second-order identity tensor and both Grad and Curl operators are defined with respect to the reference position $\mathbf{X}$. The deformation gradient $\mathbf{F}$ must satisfy the impenetrability condition $J = \det \mathbf{F} > 0$ at all times and positions. The above definitions satisfy automatically compatibility and the absence of magnetic monopole condition, i.e.,

$$ \mathrm{Curl} \mathbf{F} = \mathbf{0}, \quad \mathrm{Div} \mathbf{B} = \mathbf{0}, \qquad (2) $$

as well as jump conditions across perfect interfaces (or boundaries) such that

$$ \llbracket \mathbf{F} \rrbracket \times \mathcal{N} = 0, \quad \llbracket \mathbf{B} \rrbracket \cdot \mathcal{N} = 0. \qquad (3) $$

We introduce next the first Piola or engineering stress $\mathbf{S}$ and the magnetic field strength (or H-field) $\mathbf{H}$ as conjugate fields to $\mathbf{F}$ and $\mathbf{B}$, respectively. Conservation of linear and angular momentum (and in absence of mechanical body forces and inertial effects) as well as Ampère's law (in the absence of electric currents or charges) lead to the pointwise differential equations

$$ \mathrm{Div} \mathbf{S} = \mathbf{0}, \quad \mathrm{Curl} \mathbf{H} = \mathbf{0}, \qquad (4) $$

and jump conditions across perfect interfaces (or boundaries)

$$ \llbracket \mathbf{S} \rrbracket \mathcal{N} - \mathbf{T} = \mathbf{0}, \quad \llbracket \mathbf{H} \rrbracket \times \mathcal{N} = \mathbf{0}. \qquad (5) $$

Here, $\mathbf{T}$ denotes the mechanical traction in the reference configuration.

Use of standard push forward transformations allows to write the stress and magnetic fields in the current configuration, i.e., their Eulerian counterparts read (Ogden, 1997; Dorfmann and Ogden, 2004; Kankanala and Triantafyllidis, 2004)

$$ \sigma = \frac{1}{J} \mathbf{S} \mathbf{F}^T, \quad \mathbf{b} = \frac{1}{J} \mathbf{F} \mathbf{B}, \quad \mathbf{h} = \mathbf{F}^{-T} \mathbf{H}. \qquad (6) $$

For later use in presenting the results, one may also define the current magnetization vector $\mathbf{m}$ via the constitutive relation

$$ \mathbf{b} = \mu_0(\mathbf{h} + \mathbf{m}) \quad \text{or} \quad \mathbf{m} = \frac{1}{\mu_0}(\mathbf{h} - \mathbf{h}), \qquad (7) $$

where $\mu_0$ is the magnetic permeability of vacuum, air or non-magnetic solids.

This last expression is a definition of the magnetization vector in the current volume $V$, which however is not defined on its boundary $\partial V$. By definition $\mathbf{m} = \mathbf{0}$ in a non-magnetic body. This implies that as a quantity is insufficient to describe the presence of magnetic lines (in the sense of Maxwell) in the surrounding air or in a non-magnetic solid. Henceforth, we will focus on the original Maxwell fields $\mathbf{B}$ and $\mathbf{H}$ that are related via linear and/or nonlinear constitutive laws, while $\mathbf{m}$ will serve as a quantity that can be readily computed by Eq. (7). In fact,
<figure><img src="image_3.png" /><figcaption>Fig. 1. Scanning electron microscope (SEM) image of the as fabricated isotropic MRE sample.</figcaption></figure>

in Mukherjee et al. (2021), Mukherjee and Danas (2022) and Danas (2024), it was shown that **m** may be directly related to an internal state variable in the general case of dissipative magnetic solids, which is beyond the scope of the present work that focuses only on magnetically soft (i.e., non-dissipative) MREs.

# 3. Material fabrication and experiments

## 3.1. Fabrication of soft isotropic MREs

We fabricate the mechanically and magnetically soft magnetorheological elastomers (MREs) by dispersing carbonyl iron particles (CIP CS, BASF) of 6-8µm average diameter in an elastomer matrix composed of silicone rubber and silicone oil (378 364, Sigma-Aldrich Co. Ltd.). The composition of the fabricated isotropic MRE tested for the characterization is as follows; carbonyl iron (70 wt%, 24.07 vol%), silicone rubber (15 wt%, 33.75 vol%), and silicone oil (15 wt%, 42.18 vol%). The silicone rubber is made from the base material and curing agent (KE-1241 and CLA-9, Shin-Etsu Chemical Co. Ltd.) mixed with a weight ratio of 10:1, respectively. The addition of silicone oil decreases the mechanical modulus (modulus of transverse elasticity) of the pure elastomer matrix down to approximately 13 kPa. During the fabrication process, we first mix and stir in a beaker the carbonyl iron particles, silicone rubber and silicone oil at defined concentrations. Subsequently, we place the mixture in a vacuum chamber to eliminate air bubbles. We pour the vacuumed mixture into a non-magnetic mold of internal dimensions of 65 mm × 12 mm × 2 mm allowing to obtain a rectangular sheet MRE sample. During the curing process, we continuously rotate the mold by use of a servo motor at 60 rpm to avoid the settlement of the iron particles. In addition, the sample is heated with a heat gun set at 80 °C for 60 min to accelerate the curing. A typical microstructural image of the fabricated isotropic MRE sample obtained by a scanning electron microscope (SEM, VE-9800, Keyence) is displayed in Fig. 1, where we observe a fairly uniform and isotropic dispersion of the iron particles in the elastomer matrix.

## 3.2. Experimental setup

The characterization bench (shown in Fig. 2a) consists of a rectangular-shaped electromagnet comprising a yoke made of laminated electrical steel with a cross section of 12 mm × 50 mm and an excitation coil (wire diam.= 0.6 mm, N = 1730 turns, R = 22.9 Ω), as shown in Fig. 2b. The yoke has a 5.8 mm air gap opening allowing to sandwich two MRE samples in-between separated by a mechanically stiff magnetic steel plate of 1 mm thickness. We bond the MRE samples

to the upper and lower yoke and steel plate, respectively using adhesive layers of 0.2 mm thickness, as sketched in Fig. 2c. The horizontal displacement of the magnetic steel plate generates a shear straining of the MREs. We actuate the magnetic steel plate using a slider-crank system with a speed-controlled motor, employing an imposed sinusoidal displacement which can reach a maximum value of 2 mm (leading to a maximum of 100% shear strain for the 2 mm thickness samples) at a frequency ranging between 0.3–10 Hz. We measure the horizontal displacement of the steel plate with a laser sensor of 0.2 µm resolution (LC-2440, Keyence) with a laser displacement meter (LC-2400, Keyence). In turn, the force exerted on the two MRE samples is measured as the force acting on the oscillating plate by a load cell of nonlinearity within ±0.02% RO (LUH-50KF, Kyowa, Tokyo, Japan) with a dynamic-strain amplifier (DPM-711B, Kyowa, Japan). We generate the magnetic induction field in the gap by powering the excitation coil with a constant current (of precision ±0.01A). The static magnetic field is monitored by a Hall sensor of a Gauss meter (Model 502, Nihon Denji Sokki Co. Ltd., Tokyo, Japan) inserted between the yoke and the MRE sample under test. The apparent shear strain in MRE is calculated by dividing the measured plate displacement with the MRE sample height *h* = 2 mm. The corresponding engineering shear stress is evaluated by dividing the measured force by the combined contact area between the two MRE samples and the steel plate, i.e., by 2A, where A = 12×50 mm².

## 3.3. Experimental data

Fig. 3 presents a reduced set of representative experimental data with the aim to reveal the main effects of the loading triplet, i.e., (a) the applied magnetic field $b^0$, (b) the shear strain amplitude $\gamma^0$ and (c) the frequency f. Those observations will directly guide the modeling effort of the following sections. Specifically, we observe in Fig. 3a that the magnetic field leads to a substantial increase of the dissipation and thus of the viscosity at a fixed shear strain amplitude $\gamma^0$ = 0.4 and frequency f = 1 Hz. In addition, we observe a significant increase of the stress–strain slope at $\gamma$ = 0 and during unloading at $\gamma$ = $\pm\gamma^0$. Those two slopes together with the viscosity increase are evidence of the intrinsic coupling of the equilibrium, non-equilibrium and dissipative parts of the material and thus need to be addressed during the modeling effort. Fig. 3b examines the effect of the shear strain amplitude $\gamma^0$ at a fixed magnetic field $b^0$ = 0.5T and frequency f = 1 Hz. We observe, in particular, a notable change of the stress–strain slope at $\gamma$ = 0, with that slope increasing inversely with $\gamma^0$, known as Payne effect (Payne, 1961; Lion, 1996; Höfer and Lion, 2009). This effect has already been observed in various experimental data for MREs (Jolly et al., 1996; Sorokin et al., 2014; Sebald et al., 2017; Vatandoost et al., 2020; Wang et al., 2023) and was found to be a strong function of the magnetic field. This slope change, in practice, implies a form of a memory effect (which however is reversible in long times) of the MRE with respect to the peak strain. This feature is obviously of extreme complexity and in the present work a simplified approach will be considered in the following sections to provide model predictions valid for symmetric cyclic loads only. Finally, Fig. 3c shows the effect of the frequency f for a fixed magnetic field $b^0$ = 0.5T and shear strain amplitude $\gamma^0$ = 0.4. We observe that for the range between f = 0.5–5 Hz there is practically a negligible effect of the frequency on the shear stress–strain curves, which implies a form of non-Newtonian fluid response. This specific effect, which is observed even in the absence of the applied magnetic field (as shown later in Fig. 9a), will require special attention and departure from the classical quadratic viscoelasticity potentials proposed in the literature over the past forty years.

# 4. Thermodynamics and constitutive modeling

This section provides the kinematic and constitutive relations necessary for the description of the nonlinear magneto-viscoelastic MRE response. The only source of dissipation in the present problem is
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Electromagnet consisting of rectangular-shaped yoke with one air gap to place two sandwiched MRE samples on both sides of a magnetic steel plate. (a) Picture of the experimental setup. (b) Sketch of the electromagnet with the characterization bench and (c) 2D sketch of the bench and MRE sample dimensions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Experimental shear stress–strain data. (a) Effect of the applied magnetic field $b^0 = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8$T for frequency $f = 1$ Hz and shear strain amplitude $\gamma^0 = 0.4$. (b) Effect of the shear strain amplitude $\gamma^0 = 0.1, 0.2, 0.3, 0.4$ for frequency $f = 1$ Hz and applied magnetic field $b^0 = 0.5$T. (c) Effect of the frequency $f = 0.5, 1, 5$ Hz for shear strain amplitude $\gamma^0 = 0.4$ and applied magnetic field $b^0 = 0.5$T.</figcaption></figure>
the viscoelasticity of the MRE, as we ignore any magnetic dissipative effects, which are extremely small in the present materials containing CIP particles. Following standard literature in the context of mechanical viscoelasticity (Le Tallec et al., 1993; Bergström and Boyce, 1998; Reese and Govindjee, 1998; Kumar and Lopez-Pamies, 2016), we first introduce as an internal variable the second-order viscous deformation gradient tensor $F^v$, such that the total deformation tensor $F$ may be written as

$$F = F^e F^v \quad \text{or} \quad F^e = FF^{v-1}. \qquad (8)$$

Following this definition and for later use, we define the corresponding right Cauchy-Green total, elastic and viscous tensors, respectively, as

$$C = F^T F, \quad C^e = F^e T F^e, \quad C^v = F^{vT} F^v, \quad F^{vT} = F^v. \qquad (9)$$

It is important to note at this point that we assume that all rotations are associated only with the elastic part of the deformation gradient $F^e$ following the seminal works of Le Tallec et al. (1993) and Bergström and Boyce (1998).<fnref n="1" /> This readily leads to the symmetry of the viscous part $F^v$ and thus to $F^v = \sqrt{C^v}$. Note that $C \neq C^e C^v$ but instead $C = F^v C^e F^v = \sqrt{C^v} C^e \sqrt{C^v}$.

Next, we introduce the spatial velocity gradient $L$ as

$$L = \dot{F}F^{-1} = L^e + L^v = \dot{F}^e F^{e-1} + FF^{v-1} \dot{F}^v F^{-1}. \qquad (10)$$

where $L^e$ and $L^v$ are the elastic and viscous velocity gradients, respectively. We assume following that all rotations are associated only with the elastic part $L^e$ such that $L^v$ may be directly related to a viscous strain-rate $d^v$ thus leading to

$$d^v = FF^{v-1} \dot{F}^v F^{-1} \quad \text{or} \quad \dot{F}^v = F^v F^{-1} d^v F, \qquad (11)$$

which may be shown to be symmetric as a consequence of the assumption of rotation-free $L^v$. By direct algebraic manipulations and use of the symmetry of $F^v$, we have that $C^v = 2\dot{F}^v F^v = 2F^v \dot{F}^v$ and thus by substitution in (11), one gets

$$C^v = 2 C^v F^{-1} d^v F. \qquad (12)$$

This last expression provides an alternative option to work with the internal variable $C^v$ instead of its original form $F^v$. Moreover, $d^v$ will be defined later in Section 4.5 via a dissipation potential and its evaluation requires the complete constitutive viscoelastic model.

## 4.1. Power balance and governing equations

Considering isothermal conditions, we define an energy density function that depends on the deformation gradient $F$, the magnetic flux $B$, and the internal variable via $F^v$, denoted as $W$ ($F$, $B$, $F^v$). The local dissipation density $D$ is given as the difference of the external power $P$ minus the rate of change of the internal energy $\dot{W}$, i.e., (Coleman and Noll, 1974; Rambausek and Danas, 2021; Lucarini et al., 2022)

$$D = P - \dot{W} = \left( S - \frac{\partial W}{\partial F} \right) \cdot \dot{F} + \left( \dot{H} - \frac{\partial W}{\partial B} \right) \cdot \dot{B} - \frac{\partial W}{\partial F^v} \cdot \dot{F}^v \geq 0. \qquad (13)$$

Owing to the arbitrariness of $\dot{F}$ and $\dot{H}$, we employ the standard Coleman-Noll-Gurtin argument to arrive at the well-established constitutive relations

$$S = \frac{\partial W}{\partial F}, \quad H = \frac{\partial W}{\partial B}. \qquad (14)$$

The last term in Eq. (13), establishes the consistency conditions allowing to define the evolution of $F^v$ (or equivalently $C^v$) by use of the generalized standard material (GSM) (Halphen and Son Nguyen, 1975) formalism, such that

$$\frac{\partial D}{\partial F^v} + \frac{\partial W}{\partial F^v} = 0. \qquad (15)$$

The thermodynamic conditions presented in this section together with the field equations and interface/boundary conditions introduced in Section 2 define the entire set of equations that need to be considered to solve a general boundary value problem in magneto-viscoelasticity. An incremental variational formulation of the problem may also be constructed based on these local set of equations. For brevity this is not demonstrated in the present manuscript, but the reader is referred to the general magneto-viscoelastic formulation presented in Rambausek et al. (2022) for more details.

Neglecting inertia and mechanical body forces, the set of governing equations for the unknown displacement field $u(X, t)$ and magnetic vector potential $A(X, t)$ in a general boundary value problem in magneto-viscoelasticity is summarized as:

$$ \left\{ \begin{array}{ll} \text{Div } S = 0, & \text{Curl } H = 0, \qquad (X, t) \in V_0 \times [0, T] \\ \|S\| \mathcal{N} - T = 0, & \|H\| \times \mathcal{N} = 0 \qquad (X, t) \in \partial V_0^{T,H} \times [0, T] \\ S = \frac{\partial W}{\partial F}, & H = \frac{\partial W}{\partial B}, \qquad (X, t) \in V_0 \times [0, T] \\ \frac{\partial D}{\partial F^v} + \frac{\partial W}{\partial F^v} = 0, & (X, t) \in V_0 \times [0, T] \\ u(X, t) = u^0(t) & (X, t) \in \partial V_0^u \times [0, T] \\ A(X, t) = A^0(t) & (X, t) \in \partial V_0^A \times [0, T] \\ F^v(X, t) = I & (X, t) \in V_0 \times [0] \end{array} \right. \qquad (16) $$

## 4.2. Constitutive relations

The materials analyzed in this work are isotropic and exhibit even magneto-mechanical coupling, i.e., the response is independent of the direction of the applied magnetic field. In addition, both $W$ and $D$ need to satisfy frame indifference (Kumar and Lopez-Pamies, 2016; Mukherjee et al., 2021). A straightforward and elegant way to satisfy these conditions is to express the energy density and dissipation in terms of properly chosen isotropic invariants. First, we define the set of invariants that will be used in the present work (this set not being exhaustive). This set comprises the minimum number of invariants that are deemed necessary to model the present MRE material and is inspired directly from the previous homogenization-guided model of Mukherjee et al. (2020).

We write first with the mechanical invariants

$$ \begin{aligned} I_1 &= F \cdot F = \text{tr}C, & J &= \sqrt{\det C} = \det F, & I_1^e &= \text{tr}C^e, & J^e &= \sqrt{\det C^e} \\ &= \det F^e, & I_1^v &= \text{tr}C^v. \end{aligned} \qquad (17) $$

The above definitions readily imply that $J^v = \det F^v = J/J^e$. The corresponding magnetic and magneto-mechanical invariants read

$$I_4 = \mathbf{B} \cdot \mathbf{B}, \quad I_5 = \mathbf{B} \cdot \mathbf{C} \mathbf{B}, \quad I_6 = \mathbf{B} \cdot \mathbf{C}^2 \mathbf{B}, \quad I_5^e = \mathbf{B} \cdot \mathbf{C}^e \mathbf{B}. \qquad (18)$$

These invariants are the standard ones employed usually in the context of s-MREs (Mukherjee et al., 2020), while $I_5^e$ is proposed here additionally in order to couple the viscoelastic strains with the magnetic fields.

Following earlier work in the context of magneto-elasticity (Lefèvre et al., 2017, 2020; Mukherjee et al., 2020), we assume an additive decomposition of $W$ of the form

$$W(C, B, C^v) = \rho_0 \Psi^{\text{eq}}(C, B) + \rho_0 \Psi^{\text{neq}}(C^e, B) + \frac{1}{2\mu_0 J} \mathbf{B} \cdot \mathbf{C} \mathbf{B}. \qquad (19)$$

Here, $\mu_0$ is the magnetic permeability of vacuum, $\rho_0 \Psi^{\text{eq}}$ and $\rho_0 \Psi^{\text{neq}}$ correspond to the equilibrium and non-equilibrium parts, while the last term serves to describe the background magnetic energy present in all domains including the air domain, as well as magnetic and non-magnetic solids. It allows for the magnetic field lines to travel in the background space. In turn, the dissipation potential is assumed to be, in general, a function of $\dot{F}^v$ and the current value of $B$ but will be further specialized for convenience to be a direct function of $d^v$ defined in (11), such that

$$D(\dot{F}^v; F, F^v, B) = D(d^v; B). \qquad (20)$$

<footnote n="1">¹ This assumption is not absolutely necessary and one may show that the final energy and dissipation potential is effectively a function only of $C^v$ as discussed in Kumar and Lopez-Pamies (2016).</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Diagram of the rheological model of magneto-viscoelasticity proposed in this study.</figcaption></figure>

In the following, we propose specific choices for these energy densities and dissipation potential. Fig. 4 shows the diagram of the different energy and dissipation terms that follow, including the relevant coupling between mechanics and magnetism.

### 4.3. Equilibrium energy

We begin by prescribing the equilibrium part of the energy, which is directly inspired from the work of Mukherjee et al. (2020) and is further decomposed into a purely mechanical, purely magnetic and a coupling magneto-mechanical part as

$$ \rho_0 \Psi^{\text{eq}}(\mathbf{F}, \mathbf{B}) = \rho_0 \Psi_{\text{mech}}^{\text{eq}}(I_1, J) + \rho_0 \Psi_{\text{mag}}^{\text{eq}}(I_5) + \rho_0 \Psi_{\text{couple}}^{\text{eq}}(I_4, I_5, I_6). \quad (21) $$

Specifically, one may choose any form for the mechanical energy that is available in the literature. In the present work, we choose to work with the Lopez-Pamies hyperelastic model (Lopez-Pamies, 2010) keeping only one term, i.e.,

$$ \rho_0 \Psi_{\text{mech}}^{\text{eq}}(I_1, J) = \frac{3^{1-\alpha}}{2\alpha} G(I_1^\alpha - 3^\alpha) - G \ln J + \frac{G'}{2} (J-1)^2. \quad (22) $$

Here, G and G' are the Lamé shear and compressibility moduli, and $\alpha > 0.5$ in order to satisfy automatically polycnvexity. This last restriction is not absolutely necessary in the strict sense but is a useful property for materials that are rank-one convex as is the case of isotropic soft MREs.

The purely magnetic energy is given in terms of the hypergeometric function, i.e.,

$$ \rho_0 \Psi_{\text{mag}}^{\text{eq}}(I_5) = -\frac{\chi}{2\mu_0(1+\chi)} I_5 \mathcal{R}_1 \left[ \frac{1}{k}, \frac{2}{k}, 1+\frac{2}{k}, -\left(\frac{\chi\sqrt{I_5}}{(1+\chi)\mu_0 m^s}\right)^k \right], \quad (23) $$

with $\chi$ and $m^s$ denoting the magnetic susceptibility and magnetization saturation of the MRE, respectively.

The coupling energy is chosen of the form

$$ \rho_0 \Psi_{\text{couple}}^{\text{eq}}(I_4, I_5, I_6) = \rho_0 \delta(\bar{I}_1^v) [(1+\theta)\hat{\psi}(I_4) + (1-\theta)\hat{\psi}(I_6) - 2\hat{\Psi}(I_5)], \quad (24) $$

with $\theta \in [-1, 1]$ and

$$ \rho_0 \hat{\psi}(I_i) = \frac{\mu_0 \beta_1 (m^s)^2}{2\beta_2} \ln \left[ 1 + \beta_2 \frac{I_i}{(\mu_0 m^s)^2} \right], \quad i = 4, 5, 6. \quad (25) $$

In these expressions, $\beta_1 > 0$ and $\beta_2 > 0$ are calibration parameters. The form of the coupling function (25) is the same with that in Mukherjee et al. (2020) except that we keep only one term inside the log function.

In turn, $\delta(\bar{I}_1^v) \in [0, 1]$ is a bounded function of an additional internal variable $\bar{I}_1^v$, which describes the history of the viscoelastic strains and in particular of the first viscous invariant $I_1^v = \text{tr}\mathbb{C}^v$ defined over a symmetric cycle as

$$ \delta(\bar{I}_1^v) = 1 - q_s^v \tanh \left[ 10^3 q_0^v (\bar{I}_1^v - 3) \right], \quad \bar{I}_1^v(t) = \max_{t>0} I_1^v(t), \\ q_s^v \in [0, 1], \quad q_0^v > 0. \quad (26) $$

In this last equation, $q_0^v$ are calibration parameters satisfying the above inequalities. This function allows to model the observed Payne effect, i.e., change of the apparent slope of the stress-strain experimental curves with respect to the overall applied strain. The bounded character of $\delta(\bar{I}_1^v)$ ranging between zero and unity together with the bounded range of $\theta$ interpolation parameter in (24) guarantees that the coupling energy $\rho_0 \Psi_{\text{couple}}^{\text{eq}}(I_4, I_5, I_6)$ remains bounded for finite values of $\mathbb{C}$ and $\mathbf{B}$. We note here that the function (26) is a mere choice sufficient for the present study but may become inaccurate for more complex loading states and thus this part of the energy needs to be further improved in the future. In particular, this choice is considered to be insufficient to describe the expected decrease of the Payne effect at large times. Nevertheless, this function may be easily modified in future studies if additional experimental data along this direction become available.

**Remark.** We would like to remark at this point that the split of the energy into a magnetic part, defined via $\rho_0 \Psi_{\text{mag}}^{\text{eq}}$ in (23), and a coupling part, defined via $\rho_0 \Psi_{\text{couple}}^{\text{eq}}$ in (24), allows to effectively split the energy contributions to the magnetization and the magnetic stress, respectively. This provides ample flexibility during the modeling approach allowing to calibrate magnetostrictive/Maxwell-stress effects independently without affecting the magnetization amplitude response. This is a requirement in the case of (quasi-)incompressible MREs as has been discussed in several experimental observations (see for instance the pre-stress independent magnetization response in Danas et al. (2012), which is also true for hard MREs as discussed recently in Yan et al. (2023) and Danas and Reis (2024)).

### 4.4. Non-equilibrium energy

The non-equilibrium part of the energy is decomposed into a purely mechanical part and a coupling part, such that

$$ \rho_0 \Psi^{\text{neq}}(I_1^e, J^e, I_5^e) = \rho_0 \Psi^{\text{neq}}_{\text{mech}}(I_1^e, J^e) + \rho_0 \Psi^{\text{neq}}_{\text{couple}}(I_1^e, J^e, I_5^e). \quad (27) $$

Following the work of Rambausek et al. (2022), we assume that the non-equilibrium part of the energy does not contribute (or only weakly) on the magnetic response of the MRE. Hence, a non-equilibrium magnetic energy is not considered at this point.

The mechanical non-equilibrium part is chosen to take the same form as the energetic counterpart and reads

$$ \rho_0 \Psi_{\text{mech}}^{\text{neq}}(I_1^e, J^e) = \frac{3^{1-\alpha_e^{(1)}}}{2\alpha_e} G_e^{(1)} [(I_1^e)^{\alpha_e^{(1)}} - 3^{\alpha_e^{(1)}}] - G_e^{(1)} \ln J^e + \frac{G'_e}{2} (J^e - 1)^2. \quad (28) $$

Here, $\alpha_e > 0.5$ is a power exponent and $G_e$ and $G'_e$ correspond to the mechanical, non-equilibrium Lamé shear and compressibility moduli. It is noted that more terms may be considered in the above definition if deemed required similar to the works of Kumar and Lopez-Pamies (2016) and Rambausek et al. (2022). In the present case, we will show that one modulus $G_e$ and one power exponent $\alpha_e$ are enough to describe the experimental data.

We choose the non-equilibrium coupling energy to be of the form

$$ \rho_0 \Psi_{\text{couple}}^{\text{neq}}(I_1^e, J^e, I_5^e) = g^{\text{neq}}(I_5^e) G_e^{(2)} \left( \frac{3^{1-\alpha_e^{(2)}}}{2\alpha_e^{(2)}} [(I_1^e)^{\alpha_e^{(2)}} - 3^{\alpha_e^{(2)}}] - \ln J^e \right), \\ \text{with } g^{\text{neq}}(I_5^e) = \frac{q_0^e I_5^e}{1 + q_0^e (q_s^e)^{-1} I_5^e}. \quad (29) $$

Here, $G_e^{(2)}$ has the units of the shear modulus and is introduced for dimensional consistency, whereby $\alpha_e^{(2)} > 0.5$, $q_0^e > 0$ and $g^{\text{neq}}(I_5^e) \in [0, q_s^e]$ is a bounded non-dimensional function of a sigmoid character. Specifically, the parameter $q_0^e$ (with dimensions $T^{-2}$) serves to describe the initial response of the function $g^{\text{neq}}(I_5^e) = q_0^e I_5^e + O((I_5^e)^2) \sim |\mathbf{B}|^2 + O(|\mathbf{B}|^4)$ as $|\mathbf{B}| \to 0$. It affects, albeit only weakly (as it should), the initial magnetization response, while it influences mainly the viscous unloading stress-strain response. In turn, $q_s^e$ sets an upper saturation value to $g^{\text{neq}}(I_5^e) = q_s^e$ as $\mathbf{B} \to \infty$.
**Remark.** We insist here on the fact that the function $g^{\text{neq}}(I_5^e)$ needs to exhibit a saturating response at large magnetic fields $\mathbf{B} \to \infty$ as well as a quadratic dependence on $\mathbf{B}$ as $\mathbf{B} \to 0$. Otherwise, discontinuous or non-physical solutions may appear at the origin $\mathbf{B} = 0$, which would lead to fictitious non-zero magnetic field strength $\mathbf{H}$ and magnetization $\mathbf{m}$ fields even in the absence of an applied magnetic field $\mathbf{B}$; a response that is incorrect for s-MREs, which can only magnetize when $\mathbf{B} \neq 0$.$^2$

## 4.5. Dissipation potential

The experimental data presented in Section 3 indicate a very weak effect of the frequency and by extension of the strain-rate in the range of values that have been considered, i.e., $f = 0.3 - 5$ Hz. Given the sinusoidal character of the mechanical excitation and the 40% shear strain amplitude applied, this frequency involves strain-rates spanning the range $[0 - 12] \text{ s}^{-1}$.

This directly implies that the response resembles more that of a non-newtonian fluid rather than that of a viscoelastic solid. Motivated by this observation, we propose a dissipation potential of a power-law form (van Dommelen et al., 2003; Danas et al., 2008b; Danas and Ponte Castañeda, 2009; Garcia-Gonzalez and Jerusalem, 2019)

$$D(\mathbf{d}^v, I_5) = \frac{\eta(I_5)}{m(I_5) + 1} |\mathbf{d}^v|^{m(I_5)+1}, \quad \text{trd}^v = 0. \qquad (30)$$

From this definition, $\mathbf{d}^v$ is a purely deviatoric symmetric tensor, whereas $\eta$ is the viscosity coefficient and $m \in (0, 1]$ the power-law exponent. When $m = 1$ and $\eta$ is not a function of the viscous strains, the material becomes linear viscoelastic similar to the works of Le Tallec et al. (1993), Bergström and Boyce (1998) and Reese and Govindjee (1998). Instead, in the limit of $m = 0$, the material becomes rate independent (Idiart et al., 2006; Danas and Ponte Castañeda, 2009). For any value between zero and one, the dissipation potential (30) leads to a nonlinear with respect to $\mathbf{d}^v$ viscoelastic response, which is that of a non-Newtonian fluid. Both $\eta$ and $m$ are also functions of the magnetic field and in particular are assumed to be functions of $I_5$. It is remarked that in addition, one could allow $\eta$ to be a function of the viscous strains or viscous stress similar to the work of Reese and Govindjee (1998) and Kumar and Lopez-Pamies (2016). In the present context, such a dependence is found to be unnecessary but larger strains or different loads may reveal such dependence and in that case, one may add it in the above potential. In turn, the non-linear dependence of $\mathbf{d}^v$ via $m$ is found to be of critical importance to capture the very weak frequency dependence.

Using the properties of the derivative of the hypergeometric function introduced in (23), we propose the following functional dependence on $I_5$ for $\eta$

$$\eta(I_5) = \eta_0 \eta_B, \quad \eta_B = 1 + \frac{s_0^\eta \sqrt{I_5}}{[1 + (s_0^\eta (s_s^\eta)^{-1} \sqrt{I_5})^{p_\eta}]^{1/p_\eta}}, \quad \eta_B \in [1, s_s^\eta] \quad (31)$$

with $s_0^{0,s} > 0$, $p_\eta \ge 1$ and $\eta_0$ denoting the viscosity in the absence of a magnetic field.

Similarly, using the same function (but with different coefficients), we write the power-law exponent as

$$m(I_5) = m_0 m_B, \quad m_B = \left\{ 1 + \frac{s_0^m \sqrt{I_5}}{\left[ 1 + \left( s_0^m (s_s^m)^{-1} \sqrt{I_5} \right)^{p_m} \right]^{1/p_m}} \right\}^{-1}, \\ m_B \in \left[ \frac{1}{s_s^m}, 1 \right]. \qquad (32)$$

Again, $s_{0,s}^m > 0$, $p_m \ge 1$ and $m_0$ denotes the power-law exponent in the absence of a magnetic field.

In both cases, the calibration parameters $s_0^{\eta,m}$ control the increase (decrease) of $\eta_B$ ($m_B$) at small values of $\mathbf{B}$, whereas the parameters $s_s^{\eta,m}$ impose an upper (lower) saturation value on $\eta_B$ ($m_B$). The powers $p_\eta$ and $p_m$ control how fast the corresponding functions reach the saturation values. The calibrated response of these functions is shown later in Fig. 10.

Using the constitutive relation (15) together with the definition for $D$ in (30), we can write by means of the chain rule

$$\frac{\partial D}{\partial \mathbf{d}^v} \cdot \frac{\partial \mathbf{d}^v}{\partial \mathbf{F}^v} + \frac{\partial W}{\partial \mathbf{F}^v} = \mathbf{0}. \qquad (33)$$

After some straightforward algebra, one can show that

$$\eta |\mathbf{d}^v|^{m-1} \mathbf{d}^v = \boldsymbol{\sigma}^{\mathrm{vd}} \quad \text{with} \quad \boldsymbol{\sigma}^{\mathrm{vd}} = \mathrm{dev} \left( \frac{1}{J} \frac{\partial \rho_0 \Psi^{\mathrm{neq}}}{\partial \mathbf{F}} \mathbf{F}^T \right) \\ = \mathrm{dev} \left( \frac{1}{J} \frac{\partial \rho_0 \Psi^{\mathrm{neq}}}{\partial \mathbf{F}^\mathrm{e}} \mathbf{F}^{\mathrm{v}-T} \mathbf{F}^T \right). \qquad (34)$$

Here, “dev” denotes the deviatoric operator whereby the viscous stress $\boldsymbol{\sigma}^{\mathrm{vd}}$ is a deviatoric symmetric second-order tensor. This operation is carried out in order to impose the incompressibility constraint $\text{trd}^v = 0$ of the dissipation potential (30). One may readily solve for $\mathbf{d}^v$ in relation (34) to obtain

$$\mathbf{d}^v = \eta^{-\frac{1}{m}} |\boldsymbol{\sigma}^{\mathrm{vd}}|^{\frac{1}{m}-1} \boldsymbol{\sigma}^{\mathrm{vd}}. \qquad (35)$$

Substitution of this last result in (12) provides an explicit evolution equation for $\dot{\mathbf{F}}^v$, i.e.,

$$\dot{\mathbf{F}}^v = \eta^{-\frac{1}{m}} |\boldsymbol{\sigma}^{\mathrm{vd}}|^{\frac{1}{m}-1} \mathbf{F}^v \mathbf{F}^{-1} \boldsymbol{\sigma}^{\mathrm{vd}} \mathbf{F}, \qquad (36)$$

or alternatively for $\dot{\mathcal{C}}^v$, i.e.,

$$\dot{\mathcal{C}}^v = 2\eta^{-\frac{1}{m}} |\boldsymbol{\sigma}^{\mathrm{vd}}|^{\frac{1}{m}-1} \mathcal{C}^v \mathbf{F}^{-1} \boldsymbol{\sigma}^{\mathrm{vd}} \mathbf{F}, \qquad (37)$$

with $\boldsymbol{\sigma}^{\mathrm{vd}}$ defined in (34)$_2$. In the present work, following the earlier implementations by Rambausek et al. (2022), we have used the evolution Eq. (37).

## 4.6. Maxwell stress

For later use in the solution of the simplified BVP in Section 5, we recall here the definition of the Maxwell stress (Dorfmann and Ogden, 2003). Following the works of Kankanala and Triantafyllidis (2004) and Danas (2017), we define first the energy function equivalence $W(\mathbf{F}, \mathbf{B}, \mathbf{F}^v) = W(\mathbf{F}, J^{-1}\mathbf{F}^{-1}\mathbf{b}, \mathbf{F}^v) = w(\mathbf{F}, \mathbf{b}, \mathbf{F}^v)$, such that

$$w(\mathbf{F}, \mathbf{B}, \mathbf{F}^v) = \rho_0 \Psi^{\mathrm{eq}}(\mathbf{F}, \mathbf{b}) + \rho_0 \Psi^{\mathrm{neq}}(\mathbf{F}^e, \mathbf{b}) + \frac{J}{2\mu_0} \mathbf{b} \cdot \mathbf{b}. \qquad (38)$$

Under a prescribed Eulerian magnetic field $\mathbf{b}$ (as is the case in our experimental setup), the total Cauchy stress may be decoupled in a mechanical and a Maxwell (magnetic) part, such that $\boldsymbol{\sigma} = \boldsymbol{\sigma}^{\mathrm{mech}} + \boldsymbol{\sigma}^{\mathrm{maxw}}$. Recalling that in the present formulation, $\mathbf{b} = \mathbf{b}(\mathbf{F}, \mathbf{B})$ from relation (6), one may write

$$\boldsymbol{\sigma} = \underbrace{\frac{\partial w}{\partial \mathbf{F}}}_{\boldsymbol{\sigma}^{\mathrm{mech}}} \mathbf{F}^T + \underbrace{\left( \frac{\partial w}{\partial \mathbf{b}} \cdot \frac{\partial \mathbf{b}}{\partial \mathbf{F}} \right) \mathbf{F}^T}_{\boldsymbol{\sigma}^{\mathrm{maxw}}}. \qquad (39)$$

After some straightforward algebraic calculations, we get (Danas, 2017)

$$\boldsymbol{\sigma}^{\mathrm{maxw}} = \mathbf{h} \otimes \mathbf{b} - \frac{\mu_0}{2} (|\mathbf{h}|^2 - |\mathbf{m}|^2) \mathbf{I} = \mathbf{h} \otimes \mathbf{b} - \frac{1}{2\mu_0} (|\mathbf{b}|^2 - 2\mu_0 \mathbf{b} \cdot \mathbf{h}) \mathbf{I} \\ = \frac{1}{J} \mathbf{S}^{\mathrm{maxw}} \mathbf{F}^T. \qquad (40)$$

<footnote>$^2$ In the recent literature, coupling functions of a saturation form such as (but not only) $\tanh(\sqrt{I_5})$ or $\tanh(\sqrt{I_5})$ have been proposed to couple the mechanical equilibrium or non-equilibrium parts with the magnetic fields. This, however, may lead to fictitious magnetization fields at $\mathbf{B} = 0$ or strong dependence of the magnetization saturation with straining, a result that is in direct contradiction with available experiments such as those in Danas et al. (2012) and multiple numerical RVE simulations such as Danas (2017), Mukherjee et al. (2020, 2021), Mukherjee and Danas (2022) and Danas (2024).</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Sketch of the simplified BVP consisting of an semi-infinite (along e₁) thin MRE layer subjected to a uniform Eulerian magnetic field **b** = b°e₂ across its thickness and a simple shear loading along direction e₁ with maximum amplitude F₁₂ = γ°. The continuity conditions of the magnetic fields together with the fact that the layer is film leads to the shown known magnetic fields b₂ and h₁ in the MRE.</figcaption></figure>

This last relation will be used in Section 5 to estimate the jump in the total stress between the MRE layer and the surrounding materials.

## 5. Simplified boundary value problem (BVP) and model calibration

In this section, we discuss the calibration of the previously proposed constitutive model with a selected set of the experiments reported in Section 3. The proposed model is used subsequently to predict a wider range of experimental data not used during calibration. Specifically, the experimental setup described in Fig. 2 requires a direct and refined modeling approach. This is necessary in order to claim that the proposed model is able to describe the underlying material response. Unlike simpler mechanical experiments, in a magneto-mechanical experiment, it is extremely difficult to maintain uniform magnetic and mechanical fields throughout the loading process, especially due to corner effects that localize the magnetic fields and thus lead to heterogeneities of both the mechanical and magnetic fields in the material under study. This makes the direct modeling of MREs non-analytical and requires in most of the cases a full-field simulation of the boundary value problem (BVP) at hand.

In this work, however, we show that the present experimental setup allows, under careful assumptions, to extract the MRE material response from the BVP experimental data. In order to validate this property of the present experimental setup, in Section 6, we cross-validate the model estimates by analyzing the full BVP using 2D FE simulations.

### 5.1. Boundary conditions and assumptions

We simplify the experimental BVP by considering that a semi-infinite (along X₁) thin MRE layer lies between two semi-infinite blocks of rigid non-deformable magnetic material, as sketched in Fig. 5. The third dimension is assumed to satisfy plane-strain conditions with all magnetic fields equal to zero along e₃. As a consequence, both the magnetic and mechanical fields are uniform in the MRE specimen. This assumption is rather strong since we will see in Section 6 that both the magnetic and mechanical fields exhibit a certain level of heterogeneity, especially near the free ends of the specimen in the actual experimental setup (see Fig. 2). Even so, we will show that the force measured by the experimental load cell corresponds to the average response of the sample and in particular to the MRE material response as if the fields were uniform.

With this assumption in mind, effectively, we replace the complex BVP with a thin slab of MRE that is infinite in two out of the three directions, i.e., X₁, X₃ ∈ (-∞, ∞). In turn, the MRE is thin in direction X₂ and together with the tangential continuity of the magnetic field

strength **h** (and **H**) introduced in (5), one may assume h₁ = 0 in the MRE sample throughout the entire deformation process. This is due to the fact that h₁ = 0 in the upper and lower parts of the device (e.g., the electrical steel blocks here, neglecting in this analytical treatment of the BVP additional layers used between the MRE and the electrical steel parts as shown in Fig. 2c). This assumption becomes increasingly more accurate with the reduction of the MRE sample thickness. In addition, a thin sample allows for a fairly uniform simple shear strain state in the pure mechanical problem. This does not always mean, however, that the fields will remain uniform at large strains in the coupling problem (see for instance the work of Lefèvre et al. (2017)).

In turn, the current applied magnetic field in the electrical steel blocks is **b** = b°e₂ and owing to the normal continuity of the current magnetic flux, [[**b**]] · e₂ = 0 across the MRE sample, b₂ = b° also in the MRE sample in this simplified BVP. Considering further a uniformly applied shear deformation with maximum amplitude |F₁₂| = γ° as imposed by the experimental device, one may write the current magnetic and deformation gradient fields in the MRE material as

$$ \mathbf{b} = b_1 \mathbf{e}_1 + b^0 \mathbf{e}_2, \quad \mathbf{h} = h_2 \mathbf{e}_2, \quad \mathbf{F} = \gamma^0 \mathbf{e}_1 \otimes \mathbf{e}_2 + \mathbf{I}. \qquad (41) $$

Use of the transformation rules (6), leads to

$$ \mathbf{B} = B_1 \mathbf{e}_1 + B_2 \mathbf{e}_2 = (b_1 - b^0 \gamma^0) \mathbf{e}_1 + b^0 \mathbf{e}_2, \quad \mathbf{H} = H_2 \mathbf{e}_2 = h_2 \mathbf{e}_2. \qquad (42) $$

Following the definition for the Maxwell stress in (40) together with the previous forms of the fields **b**, **h** and **F**, we get

$$ \sigma^{\text{maxw}} = b_1 h_2 \mathbf{e}_2 \otimes \mathbf{e}_1 + b^0 h_2 \mathbf{e}_2 \otimes \mathbf{e}_2 + p^{\text{maxw}} \mathbf{I}, \quad p^{\text{maxw}} = \frac{1}{\mu_0} (b_1^2 + (b^0)^2) - b^0 h_2, \qquad (43) $$

which leads to the corresponding 1st Piola Maxwell stress

$$ S^{\text{maxw}} = B_1 H_2 \mathbf{e}_2 \otimes \mathbf{e}_1 + B_2 H_2 \mathbf{e}_2 \otimes \mathbf{e}_2 + p^{\text{maxw}}(-\gamma^0 \mathbf{e}_2 \otimes \mathbf{e}_1 + \mathbf{I}). \qquad (44) $$

An important observation in expression (43) is that the Maxwell shear stress component

$$ \sigma_{12}^{\text{maxw}} = S_{12}^{\text{maxw}} = 0. \qquad (45) $$

This is a direct consequence of the slender dimension along e₂, which in turn leads to h₁ = 0 and thus to (45). This further implies that the measured shear stress S₁₂ in the BVP is free of additional external magnetic stresses and is only a function of magnetic stresses induced by particle interactions. This is not true for the remaining components which have Maxwell contributions.

From the previous analysis and relations (41) or (42), it is clear that the only unknown in the present problem is B₁ (or equivalently
<table><caption>**Table 1**

Triplets of input loading parameters for model calibration.</caption><thead><tr><th>{f (Hz), b<sup>0</sup>(T)} \ γ<sup>0</sup></th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th></tr></thead><tbody><tr><th>{1, 0}</th><td></td><td>✓</td><td></td><td>✓</td></tr><tr><th>{1, 0.3}</th><td></td><td>✓</td><td></td><td>✓</td></tr><tr><th>{1, 0.5}</th><td></td><td>✓</td><td>✓</td><td>✓</td></tr><tr><th>{1, 0.8}</th><td>✓</td><td>✓</td><td></td><td>✓</td></tr><tr><th>{5, 0}</th><td></td><td></td><td></td><td>✓</td></tr><tr><th>{5, 0.5}</th><td></td><td></td><td></td><td>✓</td></tr></tbody></table>

$b_1$), which may be solved for incrementally by use of the non-linear algebraic equation

$$h_1 = H_1 = \frac{\partial W}{\partial B_1} = 0. \qquad (46)$$

The solution of the BVP is first done for $B_1$ using this last equation and then the viscous strains are updated from the evolution Eq. (37) using a 5th-order Runge-Kutta method (see Kumar and Lopez-Pamies (2016) and Wijaya et al. (2023) for more details). Once $B_1$ is evaluated, one may readily compute the entire stress tensor and magnetic field strength from (14).

## 5.2. Model calibration and comparison with experiments

The aforementioned simplified BVP is analyzed numerically using a custom made Fortran code that solves the nonlinear algebraic Eq. (46) together with the evolution Eq. (37).

Following closely the experimental loading conditions, we first apply a magnetic field $b_2 = b^0$ with a prescribed rate $\dot{b}^0$ and subsequently a frequency, f, dependent shear strains, summarized as follows

$$ \left\{ \begin{array}{l} b_2(t) = \dot{b}^0 t, \quad \dot{b}^0 = 0.5T/s, \quad t_{mag} = \frac{b^0}{\dot{b}^0}, \quad t \in [0, t_{mag}] \\ \gamma(t) = \gamma^0 \sin [2\pi f(t - t_{mag})], \quad t_{mech} = \frac{5}{4f}, \quad t \in [t_{mag}, t_{mech} + t_{mag}] \end{array} \right. \qquad (47) $$

It is noted here that the applied magnetic field rate is applied fairly fast in the experiment but is not exactly controlled and thus is chosen somehow arbitrarily here. We observe, however, that for the model considered in the present work (especially the non-equilibrium energy function in (29), the value of $\dot{b}^0$ has no particular effect but one should be cautious since for very soft polymers effects may be present in the normal stress components as shown in Moreno et al. (2021).

In order to achieve proper convergence for the explicit Runge-Kutta scheme, we use a time step of $\Delta t = 0.0005$. Note further that the numerical data that are used to calibrate the model parameters are those obtained at the time interval $t - t_{mag} \in [3/4f, 5/4f]$. This corresponds to the interval $\gamma \in [-\gamma^0, \gamma^0]$ after 3/4 of a cycle has been carried out allowing the stabilization of the response.

The corresponding experimental data used for calibration are also pre-processed to have a smoother character as well as in order to make sure we have a regular interval of the stress measures at pre-defined shear strains which are exactly the same with that of the numerical ones. For this, we only keep the upper part of the stress-strain curve corresponding to $\gamma \in [-\gamma^0, \gamma^0]$. Those data are then imported in Mathematica and are fitted with a higher-order polynomial (usually of degree 10 or higher) owing to their extremely nonlinear character. The polynomial is then used to export equally spaced in $\gamma \in [-\gamma^0, \gamma^0]$ data. Thus, the experimental and numerical data are extracted at exactly the same shear strain interval and the same exact shear strain increments.

Subsequently, we form an objective error function that reads

$$ \text{Obj Fcn} = \sum_{N_{irpl}} \frac{|S_{12}^{num}(1:N_p) - S_{12}^{exp}(1:N_p)|}{|S_{12}^{exp}|}, \qquad (48) $$

where $|\cdot|$ denotes the Frobenius norm of order $N_p$ corresponding to the number of time data extracted for a single loading triplet $(\gamma^0, f, b^0)$. We use a total of $N_{irpl} = 12$ triplet loading states as reported in Table 1.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Magnetization response of the MRE material using the exponent $k$, magnetic susceptibility $\chi$ and magnetization saturation reported in Table 2.</figcaption></figure>

The optimization objective function (48) is minimized by use of the Nelder-Mead and SLSQP (Sequential Least Squares Programming) optimizers in Python SciPy library to obtain the material parameters introduced in the previous sections. Both optimizers do not require the prescription of the Hessian matrix, which is a necessary condition in our problem since it involves time dependent evolution viscous strains. In order to reach an acceptable optimized solution, a first estimate of the material parameters is obtained using the Nelder-Mead simplex algorithm, which is extremely stable but at the same time exhibits a very slow convergence especially as we approach to an optimal solution. Then, the nearly optimal vector obtained from the Nelder-Mead approach is used as an initial estimate in the SLSQP optimizers, which converges to a solution in a few iterations (usually less than fifteen).

Table 2 reports sixteen material parameters as calibrated by the experiments. The parameters shown in bold are either adapted for the relevant particle volume fraction from previous works such as $\chi$ and $m^s$ reported in Psarra et al. (2017) and Diguet et al. (2019), respectively, and $k$ in Mukherjee et al. (2020), or simply considering that the MRE is quasi-incompressible for $G'$.

Using the parameters reported in Table 2 together with $G_e^{(2)} = 1$ MPa (which was used for dimensional consistency in (29)), we trace the corresponding magnetization response defined via Eq. (23) in Fig. 6 in the absence of mechanical strains (i.e. $F = I$). It is easily observed that the present data are in the range of applied magnetic flux $b^0 \in [0, 0.8]$T and thus lie mostly in the linear part of the magnetization curve, whereas saturation effects kick in for $b^0 \ge 0.9$T. It is however well-known by now that approaching magnetic saturation coupling effects also saturate and thus the response becomes of less interest in such problems.

Fig. 7 compares calibrated and predicted model curves with corresponding experimental data for frequency $f = 1$ Hz and various applied magnetic fields $b^0$ and shear strain amplitudes $\gamma^0$. Overall, the agreement between the model and the experiments is very good. The non-elliptical form of the shear stress-strain response is partially a consequence of the power-law dissipation potential (30) as well as the non-equilibrium energy function (27). Several observations are in order. First, we observe that the magnetic field has a pronounced effect on the dissipation of the MRE, which increases substantially with $b^0$. (A more quantitative analysis is discussed later in Fig. 10.) The increase of the stress with $b^0$ is a combined effect of viscoelasticity and pure energetic (equilibrium) coupling, a feature that is usually ignored in the literature where the increase of the stress response with magnetic field is simply referred to as magnetorheological effect and is usually
<table><caption>**Table 2**
Material parameters.</caption><thead><tr><th></th><th>α (-)</th><th>G'/G (-)</th><th>χ (-)</th><th>μ₀ m² (T)</th><th>k (-)</th><th>β₁ (°)</th><th>β₂ (-)</th></tr></thead><tbody><tr><td>G (kPa)</td><td>26.93</td><td>1.41</td><td>500</td><td>0.93</td><td>0.556</td><td>6</td><td>0.46</td><td>0.695</td></tr><tr><td>G<sub>e</sub><sup>(1)</sup> (kPa)</td><td>29.21</td><td>0.96</td><td>0.97</td><td>500</td><td>1.693</td><td>0.251</td><td>0.013</td><td>0.58</td></tr><tr><td>q<sub>0</sub><sup>e</sup> (T<sup>-2</sup>)</td><td>1.42</td><td>0.422</td><td>10.8</td><td>12.15</td><td>4.18</td><td>1.18</td><td>0.3558</td><td>1</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Experiments (markers) versus model (continuous lines) predictions for frequency f = 1 Hz and various applied magnetic fields $b^0 = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.8$T for different shear strain amplitudes: (a) $γ^0 = 0.1$, (b) $γ^0 = 0.2$, (c) $γ^0 = 0.3$ and (d) $γ^0 = 0.4$. Note the different range in the y-axes set independently for each plot for better visualization.</figcaption></figure>

connected to the pure energetic modulus of the MRE. In fact, the pure energetic modulus of the MRE is identified in Fig. 7 with the slope of the curve at $γ = 0$ and is obtained in the model by introducing the equilibrium coupling energy $ρ\Psi_{couple}^{eq}$ in (24). If the MRE were not viscoelastic then the sole energetic effect on the stress would have been much weaker. We observe finally that the unloading stress state after reaching the maximum applied shear strain amplitude $γ = γ^0$ depends strongly on $b^0$. This otherwise non-trivial dependence is captured by introducing a magneto-mechanical coupling via the non-equilibrium coupling energy $ρ\Psi_{couple}^{neq}$ in (29). This energy also affects the fast increase of the viscous and thus total shear stress with applied shear strain. This non-equilibrium slope increases with $b^0$ as described by the nonlinear saturation-type $g^{neq}$ function introduced in (29)$_2$.

Fig. 8 reports the effect of the applied shear strain amplitude $γ^0$ on the shear stress response for two applied magnetic fields $b^0 = 0.3, 0.5$. We observe that as $γ^0$ increases the shear stress-strain slope at $γ = 0$ decreases. This has been related to the well-known Payne effect in elastomers and purely mechanical loads. Physically, it corresponds to a reversible, albeit with limited memory, effect owing to significant reconfiguration of the polymer chains. The interesting observation in the

present context is that this effect becomes substantially more significant with the application of the magnetic field $b^0$ as can be illustrated by the relative comparison of the curves for $b^0 = 0.3$T and $0.5$T. This effect is captured by the proposed model via the use of the internal variable $\bar{\bar{L}}_1$ introduced in the definition of the *equilibrium* coupling energy (24). This implies that in the present model this effect is only present when a magnetic field is applied since it is found to be fairly negligible in the pure mechanical case. Nevertheless, a similar term with additional parameters could have been introduced if deemed necessary in the mechanical equilibrium part of the energy defined in (22). Finally, we recall that the proposition (26) is only valid for symmetric cyclic loads, whereas a more complex functional form with potentially more history/memory variables may be needed for non-symmetric cyclic loads (see a similar case in pure magnetism discussed in Mukherjee and Danas (2019)).

Fig. 9 shows the effect of the shear loading frequency upon the shear stress-strain response for different applied magnetic fields and shear strain amplitude $γ^0 = 0.4$. As already discussed in Section 3.3, the present MRE material exhibits a very weak frequency effect, which
<figure><img src="image_10.png" /><figcaption>Fig. 8. Experiments (markers) versus model (continuous lines) predictions for frequency f = 1 Hz and shear strain amplitudes $\gamma^0 = 0.1, 0.2, 0.3, 0.4$ for two applied magnetic field amplitudes : (a) $b^0 = 0.3T$ and (b) $b^0 = 0.5T$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Experiments (markers) versus model (continuous lines) predictions for three representative frequencies f = 0.5, 1.5 Hz and shear strain amplitude $\gamma^0 = 0.4$ for three applied magnetic field amplitudes : (a) $b^0 = 0T$, (b) $b^0 = 0.3T$ and (c) $b^0 = 0.5T$. Note the different range in the y-axes set independently for each plot for better visualization.</figcaption></figure>

may also be translated to a low strain-rate effect. This directly implies that the response resembles more that of a non-Newtonian fluid or viscoplastic metal than that of a usual rubber. In order to capture that effect, we have used the power-law dissipation potential in (13) with a power exponent $m < 1$. After calibration, we obtain a fairly low power-law exponent, $m_0 = 0.251$ (see Table 2) in the absence of the magnetic field, as depicted in Fig. 9a. This value implies a very low rate dependency. As one increases the magnetic field the rate effect

also becomes more important and thus in order to maintain the rate-independency observed in the experiments, the power exponent needs to decrease even more as implied by the positive values of the constants $s_0^m$ and $s_\infty^m$ in Table 2.

In order to quantify and visualize clearly the increase of the viscosity $\eta$ and decrease of the power-law exponent $m$ with the applied magnetic field $b^0$, we show Fig. 10. In part (a) of this figure, we observe that the calibrated normalized viscosity $\eta/\eta_0$ increases almost ten-fold for
<figure><img src="image_12.png" /><figcaption>Fig. 10. (a) Normalized viscosity $\eta/\eta_0$ and (b) power-law exponent $m$ shown as a function of the applied magnetic field $b^0$.</figcaption></figure>

$b^0 = 0.8$T revealing the strong influence of the magnetic field on the dissipative response of the MRE material. This viscosity increase explains directly the increase of the area enclosed between the shear stress–strain responses in the previous figures. In turn, in Fig. 10b, we observe conversely that the calibrated power-law exponent decreases to even lower values with the increase of $b^0$ thus leading to a more rate-independent response. A rather non-trivial observation in the context of these two curves is that the viscosity $\eta$ has a larger linear range as compared with the power-law exponent $m$. This is attributed to the fact that the exponents $p_\eta = 4$ and $p_m = 1$ (see Table 2) introduced in the saturation functions $\eta_B$ in (31) and $m_B$ in (32), respectively, take very different values. This further reveals a highly non-trivial influence of the magnetic field in different parts of the dissipation potential.

## 6. Finite element (FE) simulations and model assessment

The previous model predictions and calibration were done using the simplified BVP shown in Fig. 5. While this allows for a semi-analytical resolution of the problem and an efficient calibration procedure, it is based on the assumption of uniformity of the mechanical and magnetic fields throughout the MRE specimen and the fact that tangential field $h_1 = 0$, which is true for vanishingly thin films. Nevertheless, the real experimental setup involves additional materials between the electrical steel blocks where the magnetic field is present, whereas the MRE layer has a finite thickness of 2 mm. Moreover, the MRE layer extends all the way to the boundary of the electrical steel blocks and thus there exist corner effects amplifying the magnetic field in those regions and thus affecting locally the magnetic field distributions in the MRE specimens. In order to show unambiguously that the simplified BVP presented in the previous section constitutes a valid assumption and delivers accurate predictions for the actual experimental material response, we carry out two-dimensional finite element (FE) simulations of the actual experimental setup, called henceforth as finite element boundary value problem (FE-BVP). For this, we have implemented the previous magneto-viscoelastic model in a user-element subroutine (UEL) in Abaqus (ABAQUS, 2023) and compared the response obtained by the FE with that of the simplified BVP and the experiments.

### 6.1. Geometry, meshing and material choice

Starting from the experimental setup, we assume that the out-of-plane ($X_3$ direction) is under plane-strain conditions (i.e., $F_{33} = 1$, $F_{13} = F_{23} = F_{31} = F_{32} = 0$), while the magnetic fields are also zero $b_3 = h_3 = m_3 = 0$. This is a reasonable assumption except of course near the corners of the device, where magnetic field gradients exist. Nevertheless, in order to keep the analysis tractable, we focus on two-dimensions but include all the parts used in the experimental device

including the corners and the surrounding air. Specifically, following the actual setup presented in Fig. 2c, in Fig. 11, we show the geometry and the mesh used to analyze the problem at hand. The domain is subdivided into several parts such as the two MRE layers, the stainless steel (non-magnetic), the magnetic steel plate (magnetic) used to impose the mechanical shear load and the electrical steel blocks (also magnetic). In addition, we introduce a coil type domain (non-magnetic) at the top and bottom part of the BVP, which in the present case is used to apply a gradient of the vector potential $A_3$ along $X_1$. This is achieved by simply setting a non-zero value of $A_3 = -B_2X_1$ at the right part of the coil domain and a zero value at the left parts of the coil domains. The amplitude of $B_2$ is then raised to the desired value with a ramp and pre-defined rate 0.5 T/s as in the simplified BVP. What is important in the present case is that the electrical steel blocks are as long as the actual experimental device allowing a fairly uniform distribution of the b field throughout the yoke region. In order to have an even more faithful description of the experimental setup, we also model the adhesive layers (non-magnetic), shown in yellow color in Fig. 11. The FE-BVP is completed by including the surrounding air which is subdivided into a deformable and an undeformable part. The deformable part, due to the complexity of the FE-BVP, is assumed to have a very low mechanical stiffness, while both in-plane displacements in the undeformable part of the air are entirely blocked. The geometrical parameters used to for the FE-BVP are summarized in Table 3. We note further that in Fig. 11, the $L_1$ dimensions of the electrical steel part, MRE, adhesive and stainless steel layers are all equal. The same holds for the $L_1$ dimension of the undeformed and deformed air domains.

Finally, a graded mesh is used for the FE-BVP with a smaller element size as one approaches the MRE layers where the strains and magnetic fields are concentrated and the corresponding field gradients are higher. We use 4-node isoparametric bi-linear elements with four Gauss points, whereas the hydrostatic part of the energy is under-integrated at a single Gauss point at the mid of the element. More elaborate schemes may be applied to deal with the quasi-incompressibility of the energy but this has been deemed sufficient for the present analysis. Due to the two-dimensional character of the FE-BVP, we have a total of three degrees of freedom (d.o.f.) per node, i.e., the two in-plane displacements $u_1$ and $u_2$ and the scalar magnetic vector potential $A_3$. Finally, the mesh comprises a total of 34 675 nodes, which corresponds to 104 025 d.o.f.

We define next the material laws used to simulate the remaining parts that make up the FE-BVP. Specifically, we consider that the electrical steel blocks, stainless steel layers and magnetic steel plate are substantially stiffer mechanically than the MRE and adhesive layers. Moreover, the steel plate and electrical steel blocks are ferromagnetic with negligible hysteresis, while the stainless steel and the adhesive
<figure><img src="image_13.png" /><figcaption>Fig. 11. Two-dimensional geometry of the experimental setup. The constituting parts are shown with different colors. The dimensions are also denoted. The mesh near the central area is shown using four node isoparametric bilinear elements. The magnetic load is applied by specifying the vector potential $A_3 = -B_2X_1$ for all $X_1$ that belong to the “coil” domain. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<table><caption>Table 3
FE-BVP geometrical parameters (mm).</caption><thead><tr><th>L<sup>es</sup><sub>1</sub></th><th>L<sup>es</sup><sub>2</sub></th><th>L<sup>are</sup><sub>2</sub></th><th>L<sup>ss</sup><sub>2</sub></th><th>L<sup>adh</sup><sub>2</sub></th><th>L<sup>plate</sup><sub>1</sub></th><th>L<sup>plate</sup><sub>2</sub></th><th>L<sup>c</sup><sub>1</sub></th><th>L<sup>c</sup><sub>2</sub></th><th>L<sup>undef.a</sup><sub>1</sub></th><th>L<sup>undef.a</sup><sub>2</sub></th><th>L<sup>def.a</sup><sub>2</sub></th></tr></thead><tbody><tr><td>50</td><td>50</td><td>2</td><td>1</td><td>0.2</td><td>90</td><td>1</td><td>10</td><td>10</td><td>110</td><td>51</td><td>5.8</td></tr></tbody></table>

layers are non-magnetic. All of these components are assumed to be purely elastic with no mechanical or magnetic dissipation.

In order to minimize the definition of various materials, we use the same energy functions with those introduced for the MRE material but with different properties. Specifically, for all remaining parts we consider an energy density of the form

$$W(\mathbf{F}, \mathbf{B}) = \rho_0 \mathcal{P}_{\text{mech}}^{\text{eq}}(I_1, J) + \rho_0 \mathcal{P}_{\text{mag}}^{\text{eq}}(I_5) + \frac{1}{2\mu_0 J} \mathbf{B} \cdot \mathbf{C} \mathbf{B}. \quad (49)$$

where $\rho_0 \mathcal{P}_{\text{mech}}^{\text{eq}}(I_1, J)$ and $\rho_0 \mathcal{P}_{\text{mag}}^{\text{eq}}(I_5)$ are given by the same exact expressions defined in (22) and (23), respectively. Nonetheless, the material properties are chosen differently for the various parts. Those involve the corresponding Lamé moduli $G$ and $G'$, the susceptibility $\chi$ and magnetization saturation $m^3$. In all materials the powers are set equal to $\alpha = 1$ thus leading to a pure neo-Hookean response. The specific values used in the subsequent calculations are summarized in Table 4. It is important to note that the air domain is further subdivided in more domains which are smaller in size as we approach closer to the MRE. This allows to decrease gradually the shear modulus of the air as we move further away from the MRE specimen. This approach has been used successfully in Rambausek et al. (2022) to mitigate the artificial stiffness that may be introduced in the system by considering a non-zero air shear modulus. Alternative meshfree methods (Kumar et al., 2019) could allow setting the air stiffness exactly equal to zero in complex boundary domains.

We also note in passing that the shear modulus of the adhesive layers is chosen to describe a double face adhesive tape whose modulus ranges between 1000 kPa -2000 kPa. In the present case, we considered a value of ~ 1400 kPa, which corresponds to approximately 50 times that of the MRE modulus but even 1000 kPa does not alter the results. For the rest of the metal parts, a value that is 500 that of the MRE is found to be sufficient to describe the corresponding large stiffness contrast avoiding at the same time an ill-posed stiffness matrix of the numerical system. A calculation where all displacements are blocked in those stiff parts gave no difference. Finally, the shear loading is applied directly

on the entire steel plate in the same exact manner as described in (47). The simulations are conducted in a cluster using 40 cores per simulation and their duration ranges between 2 to 5 h each depending on the amplitude of the applied magnetic field and shear strain amplitude.

## 6.2. FE results versus experiments and model BVP predictions

Fig. 12 compares the shear stress-strain experimental data (markers termed as Exps.) with the simplified BVP (continuous lines termed as Mod.) and the full field FE-BVP (dashed lines termed as FE) predictions. The resulting shear stress and strain in the FE-BVP are extracted by considering the reaction force in the steel plate similar to the experimental measurement and divided by the contact area of both MRE samples. It is worth noting at this point that by considering the average 1st Piola shear stress in the MRE specimens (instead of the reaction force in the plate), we get the exact same result indicating that the surrounding materials (especially the adhesive layer domains) exhibit a negligible deformation. It is clearly illustrated in that figure that the FE-BVP delivers almost the same response as the simplified BVP which has already been shown in the previous section to be in very good agreement with the experimental data. Yet one can observe some small differences between the model and the FE data. Those differences are directly attributed to the heterogeneities in the magnetic and mechanical fields present naturally in the full field simulations but being absent in the simplified BVP model problem.

## 6.3. FE contours and local fields

In order to gain a better understanding of these heterogeneities, we show next a series of contours and local data of the relevant magnetic and mechanical fields. All subsequent data and contours consider the triplet of $b^0 = 0.5T$, $\gamma^0 = 0.4$ and $f = 1$ Hz loading conditions.

We begin by showing in Fig. 13 the local (a) $b_2$ and (b) $h_2$ fields at a central region that involves all the MRE specimens and surrounding
<table><caption>**Table 4**
Material parameters for remaining components in the FE-BVP.</caption><thead><tr><th></th><th>G/G<sup>MRE</sup> (-)</th><th>G'/G (-)</th><th>χ (-)</th><th>μ<sub>0</sub>m<sup>2</sup> (T)</th><th>k (-)</th></tr></thead><tbody><tr><td>Electrical steel blocks/steel plate</td><td>500</td><td>500</td><td>500</td><td>2.5</td><td>4</td></tr><tr><td>Stainless steel</td><td>500</td><td>500</td><td>0</td><td>-</td><td>-</td></tr><tr><td>Adhesives</td><td>50</td><td>500</td><td>0</td><td>-</td><td>-</td></tr><tr><td>Air</td><td>0.005</td><td>100</td><td>0</td><td>-</td><td>-</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Experiments (markers) versus model (continuous lines) and FE (dashed lines) predictions for shear strain amplitude $\gamma^0 = 0.4$ for applied magnetic field amplitudes $b^0 = 0, 0.3, 0.5$T and: (a) $f = 1$ Hz and (b) $f = 5$ Hz.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** (a) $b_2$ (T) and (b) $h_2$ (MA/m) contours at $\gamma^0 = 0.4$, $b^0 = 0.5$ (T) and $f = 1$ Hz. The contour comprises all the domains including the air showing the magnetic fields extending outside of the testing region. The black lines are sketched to show the position of the MRE layers and the shearing plate.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Shear strain $\gamma$ contours at applied magnetic field $b^0 = 0.5$ (T) and $f = 1$ Hz. (Top) $t_{mech} = 0$ or applied $\gamma = 0$ and (Bottom) $t_{mech} = 1.25$ s or $\gamma = 0.4$. The contour comprises the adhesive layers, MRE and steel plate domains. The yellow lines indicate the boundaries between the MRE, steel plate and the adhesive parts. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

materials. It is demonstrated that the $b_2$ is fairly uniform in the main gauge area where the MREs are being tested except near the ends of the specimen where magnetic field gradients are present due to edge and corner effects. The steel plate, being longer, seems to maintain a good magnetic flux continuity even outside the yoke domain. A similar observation can be made for the $h_2$ field which, contrary to the continuous $b_2$ field exhibits as expected jumps across the MRE/adhesive/plate directions. It is also interesting to observe that the magnetic fields extend well beyond the main testing area, partially due to the magnetic character of the steel plate.

In the following, we show more quantitative data from the FE-BVP focusing at the MRE parts only. Figs. 14 and 15 show FE results for the shear strain $F_{12} = \gamma$ and the total 1st Piola shear stress $S_{12}$ (MPa). It is interesting to observe in the data and the contours that the shear strain is more uniform than the shear stress. Focusing on the shear strain results, we observe a very good reproduction of simple shear state almost everywhere in the MRE specimen, with very weak dispersion only visible in the data (Fig. 15a). In turn, the shear stress contour and data (part (b) of both figures) exhibit a much larger heterogeneity especially near the corners of the specimen. Moreover, we observe that this dispersion of $S_{12}$ is already present in the first step of the simulation, that of the application of the magnetic field, where the strains are negligible. Those non-zero $S_{12}$ values are clearly attributed to the Maxwell stresses due to the magnetic fields and in particular in the presence of a small but non-zero $h_1$ component. Nevertheless, as already discussed in the context of the simplified BVP in Section 5, those stresses are expected to be zero when $h_1 = 0$ in the MRE layer. This assumption is further analyzed below.

In order to have a more quantitative aspect of the magnetic fields in the MRE specimens, we show in Fig. 16 all in plane components (a) $b_1$, (b) $b_2$, (c) $h_1$, (d) $h_2$, (e) $m_1$ and (f) $m_2$. It is clear from all these figures that all those fields exhibit a certain level of heterogeneity right from the beginning of the loading history. Starting from the **b** field components, it is clearly observed that the $b_1$ field exhibits a rather notable dispersion attaining $\pm 0.1$T, for an applied average $b_2 = 0.5$T. The corresponding dispersion of the $b_2$ field is in the range of $(-0.1, +0.05)$, but with most of the data lying around the average value. Turning now to the **h** field, which is of paramount importance in assessing the simplifying assumption in the simplified BVP discussed in Section 5, we observe that despite the thin dimension of the MRE layers ($L_{2}^{mre} = 2$ mm), several regions therein exhibit values of the tangential component $h_1$ larger than 0.01 MA/m. These non-zero values are evidently due to the fact that despite its slender geometry the MRE specimen still has a finite dimension along the $e_2$ direction and this

leads to non-zero shear stresses during the purely magnetic loading regime and subsequent dispersion during the mechanical loading. Yet, those $h_1$ values remain reasonably low, thus allowing to have an average shear stress that is impressively close to that of the simplified BVP and therefore to the experimental data (see Fig. 12). The $h_2$ field exhibits similar trends and dispersion as the $b_2$ field. Finally, the magnetization fields **m** follow very similar behavior with that of the **b** field as expected.

This numerical exercise reveals clearly the complexity of the magneto-mechanical experiments and the inherent heterogeneities present in such coupled problems. Yet, this cross-validation of the simplified model with a full field simulation of the experimental setup is more often than not ignored in the literature. Most studies to date assume uniformity of the magnetic and mechanical fields even in problems that this uniformity was shown to be very far from the truth as is the case in uniaxial experiments where there is a significant air gap between the magnetic source and the MRE specimens. These important points were raised already in Brown (1966) and more recently in Lefèvre et al. (2017) and Moreno-Mateos et al. (2022). In the latter study, a series of standard experimental setups where analyzed showing substantial differences in the measured magnetostrictions. As a result, several calibrated models in the literature do not correspond directly to the real material response but to the specific experimental geometry and material combined together. In passing, we note that a side study of the magnetostriction and magnetization response of the pure material is detailed in the Appendix A.

## 7. Concluding remarks and open questions

In the present study, we have presented a combined experimental-theoretical-numerical analysis for the constitutive characterization of the finite strain and large magnetic response of viscoelastic s-MREs. We have shown that the simple shear experimental setup proposed originally in Sebald et al. (2017) is capable of producing sufficiently uniform mechanical and magnetic fields in the MRE specimen, whose average corresponds closely to the actual material response. This allows to use these experimental data to calibrate in a semi-analytical fashion a magneto-viscoelastic model.

The proposed model includes a set of equilibrium and non-equilibrium mechanical and magnetic energy densities as well as a dissipation potential in terms of the viscous strain-rate or equivalently of the corresponding viscous right Cauchy-Green tensor. The explicit model is used to solve a simplified BVP which is introduced to approximate the actual experimental setup. This leads to an efficient and
<figure><img src="image_17.png" /><figcaption>Fig. 15. Mechanical shear field data inside the MRE layers in the FE simulations as a function of time (s). (a) Shear strain $\gamma$ and (b) engineering shear stress $S_{12}$ (kPa).</figcaption></figure>

very fast incremental resolution of the magneto-viscoelastic problem at hand. Use of an optimizer then allows to calibrate the underlying material parameters by use of a small but carefully chosen set of experimental data. The calibrated model is able to predict fairly accurately additional experimental data for various frequencies, applied magnetic fields and simple shear strains.

Specifically, we have found that the magnetic field leads to a significant increase of the viscosity in the MRE thus increasing substantially the resulting dissipation under a cyclic load. In addition, for the present MRE material analyzed, we show that its response resembles closely a non-Newtonian fluid and thus a power-law dissipation potential with a rather low value ($<$ 0.25) is required to model viscoelastic response. Such a low power exponent leads to negligible frequency effects on the shear stress–strain response. The power exponent is also found (by use of the calibrated model) to decrease even further with the applied magnetic field. Furthermore, the experiments exhibit a Payne effect, which becomes marked with the increase of the magnetic field. An additional history variable is introduced in the model to capture this otherwise very complex effect.

In order to validate the modeling approach and in particular the assumptions made in the context of the simplified BVP, we have implemented the MRE model in a user element (UEL) routine in Abaqus. This allows for the full-field simulation of the actual experimental BVP, and thus the cross-validation of the semi-analytical approach. It has been found that the FE-BVP is in extremely close agreement with the simplified BVP and by extension with the experimental data. The FE simulations allow as a by-product to analyze the local mechanical and magnetic fields in the MRE specimens and the surrounding components of the device and air. We have found that while the shear strains are rather uniform throughout the MRE specimens, the magnetic fields and by consequence the stress fields exhibit a measurable non-uniformity despite the rather thin width of the MRE specimen. This heterogeneity, however, does not affect significantly the average shear stress–strain response, which is the experimentally measured one. This non-trivial outcome reveals the complexity of carrying out experiments and even more of interpreting them in the context of magneto-mechanics. More often than not those non-uniformities are neglected in the literature thus leading to models that describe the device and not the actual material response.

The present experimental setup admits a number of improvements and extensions that are beyond the scope of the present study. Briefly, one could add normal mechanical loads to probe the normal Maxwell stresses and thus provide additional data that could further aid in the characterization of the material especially if that one is compressible (which is not the case in the present study). Moreover, additional sensors may be placed in the front surface of the MRE to measure the normal $h_2$ fields and thus provide additional information on the

evolution of the magnetic response throughout the deformation process. Finally, from the point of view of the modeling, one can further fine-tune the Payne effect predictions following perhaps the approach proposed by Höfer and Lion (2009).

In closing, it is perhaps relevant to remark that the majority of a large number of recent studies using MREs in the active control of metamaterials (Pierce et al., 2020; Sim and Zhao, 2023; Leanza et al., 2024) and mechanobiology (Shou et al., 2023; Gomez-Cruz et al., 2024) neglect the effect of viscoelasticity. This is perhaps due to its inherent complexity both experimentally and theoretically but also due to the lack of available experimental devices and models that will allow the more precise study of such more complex systems. As has been shown viscoelasticity may have extremely important effects on the magneto-mechanical response and change entirely the way we think the material behaves in time. In this regard, we hope that the present study will yield a better understanding of the viscoelastic effects in MREs via a properly designed experimental setup and convey the important message that mechanical and magnetic heterogeneities need to be taken into account in order to reach a sound material characterization free from boundary effects and magnetic field concentrations.

### CRediT authorship contribution statement

**Kostas Danas:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Masami Nakano:** Writing – review & editing, Validation, Resources, Methodology, Investigation, Data curation, Conceptualization. **Gaël Sebald:** Writing – review & editing, Investigation, Conceptualization.

### 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.

### Acknowledgments

K.D would like to acknowledge support from C.N.R.S. Ingénierie for providing the opportunity to spend one year at the ELyTMax laboratory at Tohoku University, Sendai, Japan. K.D. would also like to acknowledge the extensive use of the Cholesky cluster at Ecole Polytechnique for this study.
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Magnetic field data inside the MRE layers in the FE simulations as a function of time (s). (a) $b_1$ (T), (b) $b_2$ (T), (c) $h_1$ (MA/m), (d) $h_2$ (MA/m), (e) $m_1$ (MA/m), (f) $m_2$ (MA/m).</figcaption></figure>

## Appendix A. Material magnetostriction after eliminating periodic Maxwell stress

In this section, we present magnetostriction and magnetization results obtained by employing the elimination of the material Maxwell stress similar to the work of Danas (2017). In that work, the goal was to extract the purely homogenized response of an MRE free from boundary effects of the specimen. For this, one may extract from the total energy the part that is created by the Maxwell stresses of the surrounding material points (considering that the medium is infinite). Fig. A.17 shows magnetostriction (a,c) and magnetization (d) as a function of time as the applied magnetic field $b_2(T)$ shown for clarity also in part

(b). The applied magnetic field is applied using a sinusoidal signal, i.e., $b_2 = b^0 \sin 2\pi ft$ with $f = 0.5, 1, 5$ Hz, $b^0 = 2$T and $t > 0$. It is interesting to observe that for higher frequencies, a significant increase of the magnetostrictive stretches (Fig. A.17a,c) is observed at early times, while as time grows the stretches stabilize to a value that is almost insensitive to the frequency. The slower frequency $f = 0.5$ Hz stabilizes almost immediately after the first cycle, whereas for $f = 1$ Hz and 5 Hz multiple cycles are necessary in this theoretical exercise. In turn, we observe that in Fig. A.17d that the magnetization is insensitive to the frequency as well as on viscoelasticity in general. This seems to be a critical component of the present model which is in agreement with earlier experiments (Danas et al., 2012; Moreno et al., 2021)
<figure><img src="image_19.png" /><figcaption>Fig. A.17. (a) Model predictions for the magnetostriction $F_{22} > 1 > F_{11,33}$ and (b) the applied magnetic field $b_2$ as a function of time (in sec). (c) Magnetostriction $F_{22} > 1 > F_{11,33}$ and (d) Eulerian magnetization $m_2$ (MA.m) as a function of the applied magnetic field $b_2$ (T) using three frequencies $f = 0.5, 1, 5$.</figcaption></figure>

as well as full-field numerical simulations (Rambausek et al., 2022; Lucarini et al., 2022).

## Appendix B. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.mechmat.2024.105187.

### Data availability

Data will be made available on request.

### References

* ABAQUS, 2023. CAE/2023. Simulia, RI, USA.

* Bergström, J.S., Boyce, M.C., 1998. Constitutive modeling of the large strain time-dependent behavior of elastomers. J. Mech. Phys. Solids 46, 931–954. http://dx.doi.org/10.1016/S0022-5096(97)00075-6.

* Brown, W.F., 1966. Magnetoelastic Interactions. Vol. 9, Springer.

* Chatzigeorgiou, George, Charalambakis, Nicolas, Chemisky, Yves, Meraghni, Fodil, 2016. Periodic homogenization for fully coupled thermomechanical modeling of dissipative generalized standard materials. International Journal of Plasticity 81, 18–39. http://dx.doi.org/10.1016/j.ijplas.2016.01.013.

* Coleman, B.D., Noll, W., 1974. The thermodynamics of elastic materials with heat conduction and viscosity. In: The Foundations of Mechanics and Thermodynamics. Springer Berlin Heidelberg, pp. 145–156. http://dx.doi.org/10.1007/978-3-642-65817-4_9.

* Danas, K., 2017. Effective response of classical, auxetic and chiral magnetoelastic materials by use of a new variational principle. J. Mech. Phys. Solids 105, 25–53. http://dx.doi.org/10.1016/j.jmps.2017.04.016.

* Danas, K., 2024. Electro- and Magneto-Mechanics of Soft Solids. Springer Cham. chapter 3. CISM International Centre for Mechanical Sciences, pp. 65–157.

* Danas, K., Idiart, M., Ponte Castañeda, P., 2008b. A homogenization-based constitutive model for isotropic viscoplastic porous media. Int. J. Solids Struct. 45, 3392–3409. http://dx.doi.org/10.1016/j.ijsolstr.2008.02.007, URL https://www.sciencedirect.com/science/article/pii/S0020768308000656.

* Danas, K., Kankanala, S., Triantafyllidis, N., 2012. Experiments and modeling of iron-particle-filled magnetorheological elastomers. J. Mech. Phys. Solids 60, 120–138. http://dx.doi.org/10.1016/j.jmps.2011.09.006, URL http://www.sciencedirect.com/science/article/pii/S0022509611001736.

* Danas, K., Ponte Castañeda, P., 2009. A finite-strain model for anisotropic viscoplastic porous media: I – theory. Eur. J. Mech. A Solids 28, 387–401. http://dx.doi.org/10.1016/j.euromechsol.2008.11.002, URL http://www.sciencedirect.com/science/article/pii/S0997753808001150.

* Danas, K., Reis, P.M., 2024. Stretch-independent magnetization in incompressible isotropic hard magnetorheological elastomers. J. Mech. Phys. Solids 191, 105764. http://dx.doi.org/10.1016/j.jmps.2024.105764, URL https://www.sciencedirect.com/science/article/pii/S0022509624002308.

* Diani, Julie, Gilormini, Pierre, Merckel, Yannick, Vion-Loisel, Fabien, 2013. Micromechanical modeling of the linear viscoelasticity of carbon-black filled styrene butadiene rubbers: the role of the filler–rubber interphase. Mechanics of Materials 59, 65–72. http://dx.doi.org/10.1016/j.mechmat.2012.12.007.

* Diguet, G., Beaugnon, E., Cavaillé, J., 2010. Shape effect in the magnetostriction of ferromagnetic composite. J. Magn. Magn. Mater. 322, 3337–3341. http://dx.doi.org/10.1016/j.jmmm.2010.06.020, URL http://www.sciencedirect.com/science/article/B6TJJ-5093N8V-2/2/d5878c0abd20adc88cdb94281950f912.
Diguet, G., Sebald, G., Nakano, M., Lallart, M., Cavaillé, J.Y., 2019. Magnetic particle chains embedded in elastic polymer matrix under pure transverse shear and energy conversion. J. Magn. Magn. Mater. 481, 39–49. http://dx.doi.org/10.1016/j.jmmm.2019.02.078, URL https://www.sciencedirect.com/science/article/pii/S030488531832910X.

Dorffmann, A., Ogden, R., 2003. Magnetoelastic modelling of elastomers. Eur. J. Mech. A Solids 22, 497–507. http://dx.doi.org/10.1016/s0997-7538(03)00067-6.

Dorffmann, A., Ogden, R., 2004. Nonlinear magnetoelastic deformations of elastomers. Acta Mech. 167, 13–28. http://dx.doi.org/10.1007/s00707-003-0061-2.

Dorn, C., Bodelot, L., Danas, K., 2021. Experiments and numerical implementation of a boundary value problem involving a magnetorheological elastomer layer subjected to a nonuniform magnetic field. J. Appl. Mech. 88, http://dx.doi.org/10.1115/1.4050534.

Francfort, Gilles A, Suquet, Pierre M., 1986. Homogenization and mechanical dissipation in thermoviscoelasticity. Archive for Rational Mechanics and Analysis 96 (3), 265–293. http://dx.doi.org/10.1007/BF00251909.

Garcia-Gonzalez, D., Hossain, M., 2020. A microstructural-based approach to model magneto-viscoelastic materials at finite strains. Int. J. Solids Struct. http://dx.doi.org/10.1016/j.ijsolstr.2020.10.028, URL http://www.sciencedirect.com/science/article/pii/S0020768320304224.

Garcia-Gonzalez, D., Jerusalem, A., 2019. Energy based mechano-electrophysiological model of cns damage at the tissue scale. J. Mech. Phys. Solids 125, 22–37. http://dx.doi.org/10.1016/j.jmps.2018.12.009, URL https://www.sciencedirect.com/science/article/pii/S0022509618306872.

Ghosh, Kamalendu, Shrimali, Bhavesh, Kumar, Aditya, Lopez-Pamies, Oscar, 2021. The nonlinear viscoelastic response of suspensions of rigid inclusions in rubber: i—gaussian rubber with constant viscosity. Journal of the Mechanics and Physics of Solids 154, 104544. http://dx.doi.org/10.1016/j.jmps.2021.104544.

Ginder, J.M., Nichols, M.E., Elie, L.D., Tardiff, J.L., 1999. Magnetorheological elastomers: properties and applications. In: Smart Structures and Materials 1999: Smart Materials Technologies. International Society for Optics and Photonics, pp. 131–139. http://dx.doi.org/10.1117/12.352787.

Gomez-Cruz, C., Fernandez-de la Torre, M., Lachowski, D., Prados-de Haro, M., del Río Hernández, A.E., Perea, G., Muñoz-Barrutia, A., Garcia-Gonzalez, D., 2024. Mechanical and functional responses in astrocytes under alternating deformation modes using magneto-active substrates. Adv. Mater. 36, 2312497. http://dx.doi.org/10.1002/adma.202312497, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.202312497, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202312497.

Gonzalez-Saiz, E., Garcia-Gonzalez, D., 2023. Model-driven identification framework for optimal constitutive modeling from kinematics and rheological arrangement. Comput. Methods Appl. Mech. Engrg. 415, 116211. http://dx.doi.org/10.1016/j.cma.2023.116211, URL https://www.sciencedirect.com/science/article/pii/S0045782523003353.

Green, A.E., Rivlin, R.S., 1957. The mechanics of non-linear materials with memory. Arch. Ration. Mech. Anal. 1, 1–21. http://dx.doi.org/10.1007/BF00297992.

Haldar, K., 2021. Constitutive modeling of magneto-viscoelastic polymers, demagnetization correction, and field-induced poynting effect. Internat. J. Engrg. Sci. 165, 103488. http://dx.doi.org/10.1016/j.ijengsci.2021.103488, URL https://www.sciencedirect.com/science/article/pii/S00207722521000355.

Halphen, B., Son Nguyen, Q., 1975. Sur les matériaux standard généralisés. J. Méc. 14, 39–63.

Hiptmair, F., Major, Z., Haßlacher, R., Hild, S., 2015. Design and application of permanent magnet flux sources for mechanical testing of magnetoactive elastomers at variable field directions. Rev. Sci. Instrum. 86, 085107. http://dx.doi.org/10.1063/1.4927714, arXiv:https://pubs.aip.org/aip/rsi/article-pdf/doi/10.1063/1.4927714/14056288/085107_1_online.pdf.

Höfer, P., Lion, A., 2009. Modelling of frequency- and amplitude-dependent material properties of filler-reinforced rubber. J. Mech. Phys. Solids 57, 500–520. http://dx.doi.org/10.1016/j.jmps.2008.11.004, URL https://www.sciencedirect.com/science/article/pii/S0022509608001981.

Hossain, M., Saxena, P., Steinmann, P., 2015. Modelling the curing process in magnetosensitive polymers: Rate-dependence and shrinkage. Int. J. Non-Linear Mech. 74, 108–121. http://dx.doi.org/10.1016/j.ijnonlinmec.2015.04.008, URL http://www.sciencedirect.com/science/article/pii/S0020746215000797.

Idiart, Martín I., Danas, Kostas, Ponte Castañeda, Pedro, 2006. Second-order theory for nonlinear composites and application to isotropic constituents. Comptes Rendus Mécanique 334 (10), 575–581. http://dx.doi.org/10.1016/j.crme.2006.06.006.

Idiart, Martín I., Lahellec, Noel, Suquet, Pierre, 2020. Model reduction by mean-field homogenization in viscoelastic composites. i. primal theory. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476 (2242), 20200407. http://dx.doi.org/10.1098/rspa.2020.0407.

Jolly, M.R., Carlson, J.D., Muñoz, B.C., Bullions, T.A., 1996. The magnetoviscoelastic response of elastomer composites consisting of ferrous particles embedded in a polymer matrix. J. Intell. Mater. Syst. Struct. 7, 613–622. http://dx.doi.org/10.1177/1045389x9600700601.

Kankanala, S., Triantafyllidis, N., 2004. On finitely strained magnetorheological elastomers. J. Mech. Phys. Solids 52, 2869–2908. http://dx.doi.org/10.1016/j.jmps.2004.04.007.

Kim, Y., Yuk, H., Zhao, R., Chester, S.A., Zhao, X., 2018. Printing ferromagnetic domains for untethered fast-transforming soft materials. Nature 558, 274–279. http://dx.doi.org/10.1038/s41586-018-0185-0.

Kumar, Siddhant, Danas, Kostas, Kochmann, Dennis M., 2019. Enhanced local maximum-entropy approximation for stable meshfree simulations. Computer Methods in Applied Mechanics and Engineering 344, 858–886. http://dx.doi.org/10.1016/j.cma.2018.10.030.

Kumar, A., Lopez-Pamies, O., 2016. On the two-potential constitutive modeling of rubber viscoelastic materials. C. R. Méc. 344, 102–112. http://dx.doi.org/10.1016/j.crme.2015.11.004.

Lahellec, Noël, Suquet, Pierre, 2007. Effective behavior of linear viscoelastic composites: a time-integration approach. International Journal of Solids and Structures 44 (2), 507–529. http://dx.doi.org/10.1016/j.ijsolstr.2006.04.038.

Le Tallec, P., Rahier, C., Kaiss, A., 1993. Three-dimensional incompressible viscoelasticity in large strains: Formulation and numerical approximation. Comput. Methods Appl. Mech. Engrg. 109, 233–258. http://dx.doi.org/10.1016/0045-7825(93)90080-H.

Leanza, S., Wu, S., Sun, X., Qi, H.J., Zhao, R.R., 2024. Active materials for functional origami. Adv. Mater. 36, 2302066. http://dx.doi.org/10.1002/adma.202302066, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.202302066, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202302066.

Lefèvre, V., Danas, K., Lopez-Pamies, O., 2017. A general result for the magnetoelastic response of isotropic suspensions of iron and ferrofluid particles in rubber, with applications to spherical and cylindrical specimens. J. Mech. Phys. Solids 107, 343–364. http://dx.doi.org/10.1016/j.jmps.2017.06.017.

Lefèvre, V., Danas, K., Lopez-Pamies, O., 2020. Two families of explicit models constructed from a homogenization solution for the magnetoelastic response of MREs containing iron and ferrofluid particles. Int. J. Non-Linear Mech. 119, 103362. http://dx.doi.org/10.1016/j.ijnonlinmec.2019.103362, URL http://www.sciencedirect.com/science/article/pii/S0020746219306237.

Lion, A., 1996. A constitutive model for carbon black filled rubber: Experimental investigations and mathematical representation. Contin. Mech. Thermodyn. 8, 153–169. http://dx.doi.org/10.1007/BF01181583.

Lopez-Pamies, O., 2010. A new I1-based hyperelastic model for rubber elastic materials. C. R. Méc. 338, 3–11. http://dx.doi.org/10.1016/j.crme.2009.12.007.

Lucarini, S., Moreno-Mateos, M., Danas, K., Garcia-Gonzalez, D., 2022. Insights into the viscohyperelastic response of soft magnetorheological elastomers: Competition of macrostructural versus microstructural players. Int. J. Solids Struct. 256, 111981. http://dx.doi.org/10.1016/j.ijsolstr.2022.111981, URL https://www.sciencedirect.com/science/article/pii/S0020768322004346.

Moreno, M., Gonzalez-Rico, J., Lopez-Donaire, M., Arias, A., Garcia-Gonzalez, D., 2021. New experimental insights into magneto-mechanical rate dependences of magnetorheological elastomers. Composites B 224, 109148. http://dx.doi.org/10.1016/j.compositesb.2021.109148, URL https://www.sciencedirect.com/science/article/pii/S1359836821005291.

Moreno-Mateos, M.A., Danas, K., Garcia-Gonzalez, D., 2023. Influence of magnetic boundary conditions on the quantitative modelling of magnetorheological elastomers. Mech. Mater. 184, 104742. http://dx.doi.org/10.1016/j.mechmat.2023.104742, URL https://www.sciencedirect.com/science/article/pii/S0167663623001886.

Moreno-Mateos, M.A., Gonzalez-Rico, J., Nunez-Sardinha, E., Gomez-Cruz, C., Lopez-Donaire, M.L., Lucarini, S., Arias, A., Muñoz-Barrutia, A., Velasco, D., Garcia-Gonzalez, D., 2022. Magneto-mechanical system to reproduce and quantify complex strain patterns in biological materials. Appl. Mater. Today 27, 101437. http://dx.doi.org/10.1016/j.apmt.2022.101437, URL https://www.sciencedirect.com/science/article/pii/S2352940722000762.

Mukherjee, D., Bodelot, L., Danas, K., 2020. Microstructurally-guided explicit continuum models for isotropic magnetorheological elastomers with iron particles. Int. J. Non-Linear Mech. 103380. http://dx.doi.org/10.1016/j.ijnonlinmec.2019.103380.

Mukherjee, D., Danas, K., 2019. An evolving switching surface model for ferromagnetic hysteresis. J. Appl. Phys. 125, 033902. http://dx.doi.org/10.1063/1.5051483.

Mukherjee, D., Danas, K., 2022. A unified dual modeling framework for soft and hard magnetorheological elastomers. Int. J. Solids Struct. 257, 111513. http://dx.doi.org/10.1016/j.ijsolstr.2022.111513, URL https://www.sciencedirect.com/science/article/pii/S0020768322000725, special Issue in the honour Dr Stelios Kyriakides.

Mukherjee, D., Rambausek, M., Danas, K., 2021. An explicit dissipative model for isotropic hard magnetorheological elastomers. J. Mech. Phys. Solids 151, 104361. http://dx.doi.org/10.1016/j.jmps.2021.104361.

Ogden, R.W., 1997. Non-Linear Elastic Deformations. Courier Corporation.

Payne, A.R., 1961. Dynamic mechanical properties of filler loaded vulcanisates. Rubber Plast. Age.

Pierce, C.D., Salim, N.J., Matlack, K.H., 2024. Dynamic magneto-mechanical analysis of isotropic and anisotropic magneto-active elastomers. Experimental Mechanics 64 (9), 1601–1618. http://dx.doi.org/10.1007/s11340-024-01115-4.

Pierce, Connor D, Willey, Carson L, Chen, Vincent W, Hardin, James O, Berrigan, J Daniel, Juhl, Abigail T, Matlack, Kathryn H, 2020. Adaptive elastic metastructures from magneto-active elastomers. Smart Materials and Structures 29 (6), 065004. http://dx.doi.org/10.1088/1361-665X/ab80e4.
Pipkin, A., Rogers, T., 1968. A non-linear integral representation for viscoelastic behaviour. J. Mech. Phys. Solids 16, 59-72. http://dx.doi.org/10.1016/0022-5096(68)90016-1, URL https://www.sciencedirect.com/science/article/pii/0022509668900161.

Psarra, E., Bodelot, L., Danas, K., 2017. Two-field surface pattern control via marginally stable magnetorheological elastomers. Soft Matter 13, 6576-6584. http://dx.doi.org/10.1039/c7sm00996h.

Psarra, E., Bodelot, L., Danas, K., 2019. Wrinkling to crinkling transitions and curvature localization in a magnetoelastic film bonded to a non-magnetic substrate. J. Mech. Phys. Solids 133, 103734. http://dx.doi.org/10.1016/j.jmps.2019.103734.

Rambausek, M., Danas, K., 2021. Bifurcation of magnetorheological film-substrate elastomers subjected to biaxial pre-compression and transverse magnetic fields. Int. J. Non-Linear Mech. 128, 103608. http://dx.doi.org/10.1016/j.ijnonlinmec.2020.103608.

Rambausek, M., Mukherjee, D., Danas, K., 2022. A computational framework for magnetically hard and soft viscoelastic magnetorheological elastomers. Comput. Methods Appl. Mech. Engrg. 391, 114500. http://dx.doi.org/10.1016/j.cma.2021.114500, URL https://www.sciencedirect.com/science/article/pii/S0045782521007064.

Reese, S., Govindjee, S., 1998. A theory of finite viscoelasticity and numerical aspects. Int. J. Solids Struct. 35, 3455-3482. http://dx.doi.org/10.1016/s0020-7683(97)00217-5.

Saxena, P., Hossain, M., Steinmann, P., 2014. Nonlinear magneto-viscoelasticity of transversally isotropic magneto-active polymers. Proc. R. Soc. A: Math. Phys. Eng. Sci. 470, 20140082. http://dx.doi.org/10.1098/rspa.2014.0082.

Sebald, G., Nakano, M., Lallart, M., Tian, T., Diguet, G., Cavaille, J.Y., 2017. Energy conversion in magneto-rheological elastomers. Sci. Technol. Adv. Mater. 18, 766-778. http://dx.doi.org/10.1080/14686996.2017.1377590, arXiv:https://doi.org/10.1080/14686996.2017.1377590, pMID: 29152013.

Shou, Y., Le, Z., Cheng, H.S., Liu, Q., Ng, Y.Z., Becker, D.L., Li, X., Liu, L., Xue, C., Yeo, N.J.Y., Tan, R., Low, J., Kumar, A.R., Wu, K.Z., Li, H., Cheung, C., Lim, C.T., Tan, N.S., Chen, Y., Liu, Z., Tay, A., 2023. Mechano-activated cell therapy for accelerated diabetic wound healing. Adv. Mater. 35, 2304638. http://dx.doi.org/10.1002/adma.202304638, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.202304638, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202304638.

Sim, J., Zhao, R.R., 2023. Magneto-mechanical metamaterials: A perspective. J. Appl. Mech. 91, 031004. http://dx.doi.org/10.1115/1.4063816, arXiv:https://asmedigitalcollection.asme.org/appliedmechanics/article-pdf/91/3/031004/7055728/jam_91_3_031004.pdf.

Sorokin, V.V., Ecker, E., Stepanov, G.V., Shamonin, M., Monkman, G.J., Kramarenko, E.Y., Khokhlov, A.R., 2014. Experimental study of the magnetic field enhanced payne effect in magnetorheological elastomers. Soft Matter 10, 8765-8776. http://dx.doi.org/10.1039/C4SM01738B.

Stepanov, G., Abramchuk, S., Grishin, D., Nikitin, L., Kramarenko, E., Khokhlov, A., 2007. Effect of a homogeneous magnetic field on the viscoelastic behavior of magnetic elastomers. Polymer 48, 488-495. http://dx.doi.org/10.1016/j.polymer.2006.11.044, URL https://www.sciencedirect.com/science/article/pii/S0032386106012912.

Stepanov, G.V., Borin, D.Y., Bakhtiarov, A.V., Storozhenko, P.A., 2017. Magnetic properties of hybrid elastomers with magnetically hard fillers: rotation of particles. Smart Mater. Struct. 26, 035060. http://dx.doi.org/10.1088/1361-665X/aa5d3c.

Stepanov, G.V., Borin, D.Y., Kramarenko, E.Y., Bogdanov, V.V., Semerenko, D.A., Storozhenko, P.A., 2014. Magnetoactive elastomer based on magnetically hard filler: Synthesis and study of viscoelastic and damping properties. Polym. Sci. Ser. A 56, 603-613.

Stewart, E.M., Anand, L., 2023. Magneto-viscoelasticity of hard-magnetic soft-elastomers: Application to modeling the dynamic snap-through behavior of a bistable arch. J. Mech. Phys. Solids 179, 105366. http://dx.doi.org/10.1016/j.jmps.2023.105366, URL https://www.sciencedirect.com/science/article/pii/S0022509623001709.

van Dommelen, J., Parks, D., Boyce, M., Brekelmans, W., Baaijens, F., 2003. Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers. J. Mech. Phys. Solids 51, 519-541. http://dx.doi.org/10.1016/S0022-5096(02)00063-7, URL https://www.sciencedirect.com/science/article/pii/S0022509602000637.

Vatandoost, H., Hemmatian, M., Sedaghati, R., Rakheja, S., 2020. Dynamic characterization of isotropic and anisotropic magnetorheological elastomers in the oscillatory squeeze mode superimposed on large static pre-strain. Composites B 182, 107648. http://dx.doi.org/10.1016/j.compositesb.2019.107648, URL https://www.sciencedirect.com/science/article/pii/S1359836819334924.

Wang, B., Bustamante, R., Kari, L., Pang, H., Gong, X., 2023. Modelling the influence of magnetic fields to the viscoelastic behaviour of soft magnetorheological elastomers under finite strains. Int. J. Plast. 164, 103578. http://dx.doi.org/10.1016/j.ijplas.2023.103578, URL https://www.sciencedirect.com/science/article/pii/S0749641923000645.

Wijaya, I.P., Lopez-Pamies, O., Masud, A., 2023. A unified determinant-preserving formulation for compressible/incompressible finite viscoelasticity. J. Mech. Phys. Solids 177, 105312. http://dx.doi.org/10.1016/j.jmps.2023.105312, URL https://www.sciencedirect.com/science/article/pii/S0022509623001163.

Yan, D., Aymon, B.F., Reis, P.M., 2023. A reduced-order, rotation-based model for thin hard-magnetic plates. J. Mech. Phys. Solids 170, 105095. http://dx.doi.org/10.1016/j.jmps.2022.105095, URL https://www.sciencedirect.com/science/article/pii/S0022509622002721.