Curing spurious magneto‐mechanical coupling in soft non‐magnetic materials

Abstract The present work is concerned with the issue of spurious coupling effects that are pervasive in fully coupled magneto‐mechanical finite element simulations involving very soft non‐magnetic or air‐like media. We first address the characterization of the spurious magneto‐mechanical effects and their intuitive interpretation based on energy considerations. Then, as main contribution, we propose two new cures for the issue under consideration that completely prune the undesired spurious magneto‐mechanical coupling in non‐magnetic media. The proposed methods are compared with established methods in the context of magnetic bodies embedded in (i) air or vacuum and (ii) very soft elastic non‐magnetic media. The comparison shows that the proposed approaches are accurate and effective. They, furthermore, allow for a consistent linearization of the coupled boundary value problems, which is crucial for the simulation of compliant structures. For reproducibility and accessibility of the proposed methods, we provide our implementations with Netgen/NGSolve as well as all codes necessary for the reproduction of our results as supplementary material.


Introduction
In this work we address the fundamental numerical issue of spurious magnetic coupling in full-field simulations of extremely compliant magnetic bodies and structures.The bodies that we have in mind are composites consisting of an elastomer matrix with almost rigid ferromagnetic inclusions, i.e., magnetorheological elastomers (MREs).For experimental data on MREs we refer to Jolly et al. (1996), Bednarek (1999), Ginder et al. (2002), Martin et al. (2006), Böse and Röder (2009), Danas et al. (2012), Bodelot et al. (2018) and Garcia-Gonzalez et al. (2021) as well as to the extensive review by Bastola and Hossain (2020).These findings have inspired many interesting applications of MREs such magnetically actuated valves (Böse et al., 2012), magnetically tunable dampers (Stewart et al., 1998) and vibration isolators (Ginder et al., 2001;Hitchcock et al., 2005Hitchcock et al., , 2006) ) as well as micropumps (Cesmeci et al., 2022).With matrix materials becoming softer and bodies and structures becoming more compliant, new opportunities for applications open up in medical and biological contexts Gonzalez-Rico et al. (2021); Moreno-Mateos et al. (2022a).However, wellestablished experimental, numerical and theoretical approaches face new challenges.For experimental data in this regard we refer to Nikitin et al. (2004), Stepanov et al. (2007Stepanov et al. ( , 2008) ) and rather recently Moreno et al. (2021); Moreno-Mateos et al. (2022c), which set several opportunities for future research.
A particular breed of MREs based on remanently magnetized particles dates back at least to the exploratory experiments of Stepanov et al. (2012).Such "hard" MREs (h-MREs) have recently received increased thanks to advances in additive manufacturing of h-MRE structures (Kim et al., 2018).One prominent field of application of this new technology is soft robotics (Kim et al., 2019).Numerical modeling of h-MREs under assumption of a prescribed remanent magnetization has been addressed by Garcia-Gonzalez (2019), Zhao et al. (2019), Garcia-Gonzalez and Hossain (2021), Kadapa and Hossain (2022), Dadgar-Rad and Hossain (2022).The resulting models reduce to finite elasticity with particular loading conditions.In contrast to that, the important aspect of ferromagnetic evolution in the particles, e.g. during the magnetization process of h-MREs, has been accounted in the methodological contributions of Kalina et al. (2017) and Rambausek et al. (2022).Macroscopic material models that are able to reproduce microstructural simulations of h-MREs have been developed by Mukherjee et al. (2021) and Mukherjee and Danas (2022).For comprehensive reviews on h-MREs we refer to Lucarini et al. (2022) but also to Wu et al. (2020) where both (magnetically) "soft" MREs and h-MREs were considered.
No matter whether "soft" or "hard" MREs are considered in simulations, one almost always encounters the pervasive issue that will be the addressed by this contribution: spurious magneto-mechanical coupling in sufficiently soft or compliant magnetoactive bodies as depicted in Figure 1 (Rambausek, 2020).There one magnetized non-magnetic Figure 1: A spatially fixed, rigid and magnetized body (grey) embedded in a soft non-magnetic material (yellow).The rigid, magnetized body does not deform.However, spurious deformation (red vs. grey mesh) is observed in the non-magnetic domain in response to magnetic fields (Rambausek, 2020, Chap. 9.3.1;Figure 9.17a).
can observe excessive deformation in the non-magnetic material (yellow) embedding a clamped, rigid magnetized body (grey) where no deformation at all is expected according to the laws of physics.These effects are ubiquitous but only become relevant in extremely soft materials.This can be for example the case when a sufficiently soft non-magnetic medium is employed to model the air surrounding some magnetic specimen (Keip and Rambausek, 2016).However, it may also be observed when increasingly soft2 matrix materials are employed as is the case for example in the experiments of (Moreno et al., 2021).While both situations have the commonality of an extremely soft non-magnetic material being present in the boundary value problem, in the first case this is purely for convenient modeling of a medium with no significant stiffness at all, in the second case the soft material is physical reality.Due to that, in the former case, one has more modeling and methodological freedom than in the latter.In fact, Pelteret et al. (2016) proposed an elegant staggered approach to account for air or vacuum surrounding a magnetic body that does not make use of an auxiliary soft material and hence does not suffer from spurious coupling in the present sense.However, their approach is not applicable when a soft elastic material is indeed physically present.Also, due to its partially decoupled nature and the inherent inconsistent linearization of the discretized boundary value problem, it suffers from convergence issues in the case of highly compliant structures.As a second method that successfully circumvents spurious coupling effects we mention the staggered approach of Liu et al. (2020).Their approach is based on surface currents computed in a purely magnetic boundary value problem.In a second step, these currents are employed to compute magneto-mechanical tractions on boundaries of magnetic bodies for the mechanical subproblem.Due to the decoupled magnetostatic and mechanical subproblems, the scheme requires a certain number of interations between both of them to account for the strong coupling and attain convergence.While the procedure is effective for the analysis of slender magneto-active structures, its staggered nature and the explicit use of electric currents for force computation render the scheme less suited for adoption in existing fully coupled (monolithic) methods.A third approach for the modeling of surrounding air in finite element simulation are non-local constraints that bind the motion in the air domain to the motion of the boundary of the magnetic bodies under consideration (Psarra et al., 2019;Mukherjee et al., 2021).While this method does not introduce artificial stiffness and linearizes properly it leads to a much denser linear system.Furthermore, it is difficult to implement for general geometries and arrangements of multiple magnetic bodies.For some additional insights on existing monolithic and staggered schemes we refer to Rambausek et al. (2022).
The methods proposed in this contribution exploit the non-magnetic nature of the soft material (surrounding air or soft carrier/matrix material) to obtain formulations that integrate well with existing magnetoelasticity frameworks, i.e. it is not much more difficult to implement than the naive monolithic schemes, for both the surrounding air and the soft matrix case.Furthermore, the schemes obtained linearize properly and at the same time do not suffer from spurious coupling and are thus promising candidates for widespread adoption.To facilitate this process we provide our implementation with Netgen/NGSolve (Schöberl, 1997(Schöberl, , 2014(Schöberl, , 2022) ) as well as codes for our numerical examples as supplementary material.
The remainder of this contribution is structured as follows.In Section 2 we outline the fundamental equations of finite magneto-elasticity as well as important considerations on magneto-mechanical interactions.Section 3 then presents the spurious magneto-mechanical coupling that is the motivation for the present work and discusses the underlying general pattern.Then, in Section 4 we arrive at the heart of our contribution where we propose two approaches that eliminate spurious magneto-mechanical interactions in non-magnetic domains.We sketch numerical implementations of the schemes and assess their accuracy by a comparison with existing methods as well as their convergence under mesh refinement.After that we showcase both proposed approaches in a challenging example involving a magnetic solid and a very soft non-magnetic solid under gravitational and magnetic loads embedded in an air-like medium.We close our contribution with the conclusion in Section 5.

Theoretical background
For the problem of spurious magneto-mechanical coupling investigated in this work, we restrict ourselves to quasi-static finite elasticity coupled to magnetostatics in absence of free currents.We first recall the fundamental theory for the mechanical part, then the magnetostatic and finally bring both fields together.Figure 2 provides a overview of the setting considered.There we depict magnetic bodies (grey;  , ) and non-magnetic The overall setting considered in this work.Magnetic bodies (grey;  , ) and non-magnetic bodies (yellow;  ) in "free" space (light cyan; vacuum, air-like medium) exposed to a uniform external magnetic field ∞ .We employ the symbol  ∞ for the full domain including all sub-bodies and  a for the "remaining" empty or air-like domain.
bodies (yellow;  ) embedded in "empty" space  a (light cyan; vacuum, air-like medium).The bodies are exposed to a uniform external magnetic field ∞ and may have mechanical support.While not indicated, we may of course consider mechanical body forces like gravity as well.The "full" domain consisting of  ∞ and all (solid) bodies  is denoted as  ∞ .The (computational) outer boundary  ∞ is indicated by a dashed line.Physically, there is no boundary but in finite element simulations the infinite space domain is truncated3 at a larger distance from the bodies under consideration (Bustamante et al., 2011;Miehe and Ethiraj, 2012;Keip and Rambausek, 2016;Schröder et al., 2022).For definiteness, we mechanically fix the outer boundary  ∞ .

Quasi-static finite hyperelasticity
In the context of finite strains we have to distinguish between the Eulerian or current configuration and the Lagrangian or reference4 configuration.The former simply refers to the (deformation) state of media as an external observed perceives them, whereas the latter takes the perspective of material points.Following a common convention, we denote quantities in the Eulerian configuration by lowercase symbols, whereas uppercase symbols are employed for the Lagrangian configuration.In some instances, we employ subscript "0" for Lagrangian and subscript "t" for Eulerian quantities for definiteness.The Eulerian and Lagrangian configurations are connected through the map , that assigns to each Lagrangian position (a material point) a Eulerian position = ( ).The corresponding tangent map or Jacobian = .In Cartesian coordinates is often written as and commonly referred to as deformation gradient.Its determinant = Det is a measure for the change of volume elements where d is the Eulerian and d the Lagrangian volume element.On the kinetic side, the balance of momentum in a quasi-static setting reads in the Eulerian and Lagrangian setting, respectively, where σ denotes the Cauchy stress, the "first Piola-Kirchhoff stress" and b denotes the body force density per Eulerian volume.The surface-integral-preserving Piola-type transform relating σ and is The Eulerian body force density b transforms as The corresponding boundary conditions are given as where is a given Eulerian traction, i.e. a surface force density per Eulerian area and is the outward-pointing unit normal.The complicated factor on the right is the ratio of area elements d ∕ d which transforms the Eulerian to a Lagrangian traction.
In the context of hyperelasticity, the stresses are obtained as partial derivatives of a strain energy density per mass ( ) where and 0 are the mass densities per current (Eulerian) and initial (Lagrangian) volume.The latter are related through As the final modeling ingredient we introduce the strain energy density per Lagrangian Ψ 0 ( ) = 0 ( ) volume, where 0 is the respective mass density.Correspondingly, the strain energy density per Eulerian volume is Ψ ( ) = ( ).The balance of momentum is accompanied the balance of moment of momentum expressed as i.e. the symmetry of the Cauchy stress.In the present setting, ( 9) is ensured by an objective (material frame indifferent) strain energy density.Material frame indifference is achieve, for example, by parameterizing the strain energy density in terms of the "right Cauchy-Green" tensor The last term clearly is symmetric and thus σ is symmetric too.
Considering the case where the external volume and surface force densities are conservative, we introduce the respective loading potential densities (per Eulerian volume) and via With these at hand, we may formulate a variational principle for finite hyperelasticity where we remark that the integrals can be individually transformed to the Lagrangian configuration for computational convenience with the help of (2), ( 8) and the area transform factor in (6).For example, the fully transformed version reads

Magnetostatics
The magnetostatic field equations reduced to the scope of the present work read where is the magnetic −field (magnetic flux), is magnetic ℎ−field (magnetic intensity).From (15) it is clear that and can be expressed in terms of potentials as with being the (magnetic) vector potential and the auxiliary scalar magnetic potential.While not of real physical significance, is particularly convenient for computations.It is possible to show that the equations introduced so far can be obtained from two equivalent variational principles.The first, to which we refer as magnetic energy principle reads where Ψ ( , ) is the magnetic energy density per Eulerian volume, a quantity to be specified depending on the medium occupying the respective point in space.At a solution ̂ we have Dual to that we have the magnetic co-energy principle where ( ) is the magnetic co-energy density, which is obtained from ( ) through a Legendre-Fenchel transform At a solution ̂ we have

Coupled magnetoelasticity
In order to couple magnetostatics to finite elasticity we establish the following transforms We note that the transforms involving − preserve line integrals, those involving −1 preserve surface integrals and those involving −1 preserve volume integrals.
In what follows, we exploit that any material-frame indifferent Ψ ℎ 0 ( = T ⋅ , ) can be recast to the form Ψ 0 ( , , ) and Ψ 0 ( = −1 ⋅ , ) to Ψ 0 ( , , ) (Mukherjee et al., 2021;Mukherjee and Danas, 2022).The parameterizations of Ψ 0 and Ψ 0 have the advantage that they directly generalize to the magnetomechanical case since they depend on both (mechanical) kinematic and magnetic quantities.With these at hand, we combine the mechanical ( 14) and magnetic principles ( 17) and ( 19) to the coupled variational principles (Bustamante et al., 2008) and whereby we note that, as a consequence of fixing the outer boundary  ∞ , the loading terms on that boundary in ( 23) and ( 24) indeed take this form with a slight abuse of notation.Moreover, outside  the energy densities Ψ 0 ( , ) and Ψ 0 ( , ) reduce to the respective vacuum energy densities discussed in Subsection 2.4.

Specific relations for non-magnetic media
In empty space (vacuum) and, to very good approximation, in all (technically) non-magnetic5 media the fields and are related through with 0 being the vacuum permeability.The corresponding (co-)energy densities are given as and from which (25) can be easily computed via (18) 2 or (21) 2 , respectively.An important consideration in this context is that non-magnetic media do not experience any "magnetic forces".This means that they will not move or deform under magnetic field.Also, in general they remain non-magnetic under deformation.This is reflected by the following observation: Consider any of the energy densities (26) or ( 27) and employ them in the respective magnetomechanical variational principle (23) or (24).For this purpose, the densities are expressed in their Lagrangian form and Then, irrespective whether magnetostatic energy or co-energy is used, the resulting Cauchy-type stress field in this particular case is which is commonly known as the Maxwell stress6 .
The key properties of the Maxwell stress are (i) that its Eulerian form only depends on magnetic quantities but no mechanical (kinematic) ones and (ii) that it is divergence-free for any magnetostatic solution in nonmagnetic media.Combining both reveals that there cannot be any magnetic force densities in such a medium irrespective of the deformation state.An equivalent viewpoint is that Eulerian magnetostatic energy an coenergy are invariant with respect to deformation in the interior of a non-magnetic medium.
Remark 1.The simplicity of the vacuum constitutive relation (25), in particular the scalar appearance of 0 , hides the fact that and are of different geometric nature.By "pulling-back" this equation to the Lagrangian configuration, we obtain which reveals the tensorial object represented by 0 .In terms of differential geometry, corresponds to a one-form whereas corresponds to a two-form, which roughly means that the former is a "gradient-type" field and the latter is of "curl-type".Both can be expressed as vectors, but they transform differently.A tensorial object that relates fields of different character in the above sense is an instance of a Hodge operator.

Magnetic media and magnetic forces
In accordance with the above "definition" of a non-magnetic medium we regard any medium for which (25) does not hold in general as magneticor magnetized.The mismatch is connected to the magnetization which appears as such that non-magnetic media have ≡ 0.
It is important to note here that the generic constitutive relations (18) 2 and ( 21) 2 still hold.In particular (Bustamante et al., 2008) and where Ψ ( , ) and Ψ ℎ ( , ) are total energy densities including the vacuum energy densities whereas Ψ ( , ) and Ψ ℎ ( , ) are contributions from matter.
Concerning stresses and force densities due to the presence of magnetic fields, we consider (Kankanala and Triantafyllidis, 2004, Table 1) and where we note that results for the divergence were obtained for curl = 0.The fact that the stress contributions are obviously not divergence-free can be interpreted as the presence of magnetic volume force densities.
The corresponding surface force densities at an interface to a non-magnetic medium are (Kankanala and Triantafyllidis, 2004, Table 1) and where is the outward pointing normal vector.
Remark 2. In magnetically anisotropic media, i.e. media where is not necessarily parallel to , the magnetomechanical Cauchy-type stresses ( 35) and ( 37) are not symmetric but only the total Cauchy-type stress ( 11) is.This corresponds to the presence of magnetic torques.Examples for such materials are MREs with anisotropic particle distributions (Boczkowska and Awietjan, 2012;Danas et al., 2012) and h-MREs when remanently magnetized.We refer to Mukherjee et al. (2021) for a related discussion of the latter case.
When considering two rigid magnetic bodies that are separated by a non-magnetic medium, there will be some kind of force transmitted between these two bodies.The transmission indeed happens through the non-magnetic medium.In fact, the net magnetic force  mag on a magnetized body can be computed by the surface integral of the normal component of the Maxwell stress over any closed surface  around the body of interest (Brown, 1966) See Figure 3 for an illustration.Orthogonal to that is the important case of a magnetized body that does not experience a net magnetic force.This can be a body with remanent magnetization that is remote from other magnetic media or a uniformly magnetized body in a uniform external field.Then, even in absence of a net force, there will still be jumps in the normal component of the magneto-mechanical stress across the boundary of the magnetic body that depend on the magnetization (see equations ( 39) and ( 39)).In other words, the body still experiences tractions that are equilibrated in the sense that the net force vanishes.However, they may still cause deformation.The overall effect of these traction on a body depends on its shape, which is the fundamental mechanism behind the shape-dependent magneto-mechanical response of MREs (Keip and Rambausek, 2017).A corresponding energetic point of view has been put forward by Rambausek and Keip (2018).Such settings are also formidable test cases for the present study because of the strong coupling of boundary deformation and magnetic response even in geometrically "simple" settings.

Spurious magneto-mechanical coupling
While the theory clearly states that non-magnetic media do not deform in reaction to magnetic fields, one may observe the contrary in numerical simulations.For the purpose of a brief demonstration, we consider the magneto-mechanical (quarter) boundary value problem (BVP) of a circular quasi-rigid magnetic body surrounded by a rather soft non-magnetic medium such as a very soft elastomer or some extremely soft auxiliary material representing air.The geometrical setting is illustrated by lated in terms of both coupled variational principles ( 23) and ( 24), which are discretized with second-order isoparametric Lagrange elements for both the deformation (displacements) and the out-of plane vector potential component or the scalar magnetic potential, respectively, as commonly done in finite magneto-and electro-elasticity (Vu et al., 2007;Javili et al., 2013;Keip et al., 2014;Miehe et al., 2016;Keip and Rambausek, 2016;Pelteret et al., 2016;Danas, 2017).The resulting fully coupled monolithic nonlinear system obtained is solved by a Newton-Raphson procedure.Before we continue, we need to clarify some terms.With "fully coupled" we simply refer to schemes which do not neglect any coupling that comes out of the variational principle to be discretized, i.e. no physical effects present in the theory are neglected or dropped for computational considerations or convenience.Second, "monolithic" means that one solves the equations representing the fully coupled discretized system for all primary variables at the same time in the whole discretized spatial domain.For example, in the case of ( 24), one solves for the pair { , } in all magnetic and non-magnetic media under consideration.
Below we consider both the energy ({ , }) and energy-co-energy formulations ({ , }).For both the magnetic and non-magnetic domain we consider the energy densities and where vacuum permeability is 0 = 0.4 µT m A −1 .For either formulation, the material parameters employed for the practically rigid magnetic domain are { , ′ , } = {1 × 10 4 kPa, 5 × 10 5 kPa, 10}.For the non-magnetic medium we choose { , ′ , } = {1 × 10 −5 kPa, 5 × 10 −4 kPa, 0} in when using the energy formulation but { , ′ , } = {1 × 10 −3 kPa, 5 × 10 −2 kPa, 0} in the energy-co-energy formulation.For the non-magnetic domain, we for the purpose of demonstration employ different mechanical parameters because the spurious coupling is more pronounced for the energy-co-energy formulation in this example.Figure 4(b) and (c) show color-contours of the magnetic field magnitude in the deformed configuration obtained with the magneto-mechanical energy formulation (23) and the energy-co-energy formulation (24), respectively.The blue mesh corresponds to the deformed configuration while the light grey mesh is the undeformed one.Both meshes coincide in the practically rigid magnetic particle.From the physical theory one would not expect any deformation or displacements in the non-magnetic domain, since the magnetic body is spatially fixed and rigid whereas the deformable material is non-magnetic and clamped at the outer boundary  a .However, the dashed-red triangles marking undeformed triangle cells which undergo significant deformation (dashed red versus solid blue) show that the opposite is the case.This clearly is against the physical theory and, as we argue below, a universal trait of fully coupled monolithic discretizations.The fundamental problem with such approaches is that the solution fields differ for each discretization, which is an inherent property of any numerical method in general.This means that, the values of the governing energies or potentials evaluated for numerical solution states change with the discretization.The crucial step is now to consider the magnetostatic (co-)potentials and respectively, where is considered as a "given" parameter in Π and Π ℎ .Recall from Sections 2.2 and Sections 2.3 that -for a given -the magnetostatic problem is solved for by minimizing Π and for by maximizing Π ℎ , respectively.However, we also solve for via the minimization of the total (co-)energy7 Π | .When we were not in possession of numerical but exact magnetostatic solutions, Π | would be invariant with respect to in the nonmagnetic domain (the magnetic domain is considered to be rigid in the present example).But as our numerical solutions are not exact, some deformation has to happen8 as we have seen in Figure 4b,c.The energetic considerations are illustrated by Figure 5, where subplot (a) shows Figure 5: Magnetostatic (co-)energy Π |ℎ for a rigid magnetic particle embedded in a non-magnetic solid depending on the applied field magnitude ∞ (a).The respective difference to the (co-)energy Π |ℎ rigid obtained for a rigid non-magnetic domain (b).
the magnetostatic potentials Π |ℎ at a solution state and (b) depicts the difference to the rigid case.The graphs in Figure 5a are practically mirrored (neglecting FE errors) as is expected from the duality of the energy and the co-energy formulations in the case of linear magnetic materials.In Figure 5b we observe that both potentials Π and Π ℎ are smaller for the deformable non-magnetic medium than for the rigid case, i.e. there was a minimization of the magnetostatic energy through deformation.In the energy formulation this is formally favorable as we interprete a smaller computed total energy as being closer to the (true) solution than a higher total energy.But of course, this minimization is only possible because of numerical errors and the resulting deformation is entirely unphysical and undesirable.In the energy-co-energy case, we do not have any "formally" favorable outcome.In particular, the deformation pattern is to our experience in general quite different from the one obtained for the energy formulation.Moreover, it is typically more pronounced -recall that we had = 1 × 10 −3 kPa in the energy-co-energy formulation but = 1 × 10 −5 kPa in the energy formulation.In the present example, the energy-co-energy would crash for applied field magnitudes somewhere between around ∞ = 0.5 T and ∞ = 0.6 T. In contrast, the energy formulation runs through until at least ∞ = 1.0 T. While the latter appears to be more robust, it might also crash eventually.However, if the soft non-magnetic medium is indeed a physical material and not some auxiliary model for an air-like environment, any spurious deformation must be avoided.In the next section we present two approaches that achieve this goal.
Remark 3. The above considerations have parallels to those behind variational r-adaptivity for finite elasticity (Thoutireddy and Ortiz, 2004;Kuhl et al., 2004;Askes et al., 2004).In that context, the "spatial motion" (the "usual") finite strain elasticity problem on the Lagrangian configuration plays a role comparable to the magnetostatic problem in the Eulerian configuration.The Lagrangian mesh optimization, on the other hand, is governed by the "material motion" problem which is parallel to the "spatial" motion problem -in a sense Eulerian mesh optimization -in the present context.We furthermore point out that the "material motion" problem corresponds to configurational mechanics and the respective forces are configurational forces (Gurtin, 1994(Gurtin, , 1999;;Eshelby, 1999).They are expected to vanish everywhere except at material defects such that one may attain the viewpoint that the discretization of a body indeed introduces defects.
4. Two novel approaches to cure spurious magneto-mechanical coupling in non-magnetic media Without loss of generality, we from now on restrict ourselves to the use of the magnetostatic co-energy and the variational principle (24).Nevertheless, we emphasize that all derivations and observations presented below have their equivalent counterpart when ( 23) is employed instead.In particular, we recall that the Maxwell stress has a unique notion in non-magnetic media.As a result, the methods introduced in this Section are directly applicable to both the energy and energy-co-energy formulation.
Below we shall distinguish two cases.First, we consider magnetic bodies surrounded by vacuum or air, i.e. a non-magnetic medium of practically negligible stiffness.Then we do not have to care too much about the actual deformation outside the bodies.What matters is that the deformation acceptable from a numerical perspective.The second case is concerned with magnetic bodies that are embedded in a very soft carrier medium such as an extremely soft elastomer, a gel or some biological tissue.In that case one expects a specific mechanical response of non-magnetic material such that we are not free to choose the mechanical parameters.Because of this constraint, we regard the latter case as more difficult.
To our best knowledge, up to now, all remedies to spurious magneto-mechanical coupling in non-magnetic media are concerned with spurious coupling in vacuum or air-like domains.Therefore, we start with this important case, where we first briefly review existing approaches and then propose two novel cures to this issue.They will be assessed by a comparison with each other and the staggered scheme of Pelteret et al. (2016).After that demonstration, in Subsection 4.5 we will turn to the extension of the proposed methods to extremely soft, non-magnetic solids.
Before jumping into the details, we remark that we only consider the case when spurious deformations are an issue only in the non-magnetic domains but not in the magnetic ones.In other words, throughout this work we assume magnetic domains to be sufficiently stiff such that spurious deformations are negligible.We anticipate though, that this might not always be the case.

Eliminating spurious coupling in vacuum and air-like domains
In case of vacuum or air surrounding the body of interest, one might look into schemes which just do not solve the coupled problem in a monolithic but in a staggered way along the lines of Pelteret et al. (2016).In their scheme, one first solves the fully coupled problem in a monolithic may but with deformation blocked in the interior of the non-magnetic domain.By that one ends up with a thin deformable layer between the boundary of the solid bodies and the first FE nodes in the vacuum domain with vanishing stiffness.As a consequence, this first part of problem is free of any artificial elastic model for the air or vacuum domain.The second step consists of smoothing the displacement field in the vacuum domain by an auxiliary problem.While this adaption does not have any effect on the stiffness experienced by the solid bodies from their "empty" surrounding, there is an effect on the magnetostatic solution due to the change of the Eulerian configuration.And since the problem is physically coupled, the mesh adaption step perturbs the equilibrium computed in the first, coupled step.This can be a disadvantage when very soft bodies or compliant (slender) structures are considered.The situation can be improved by iterating between the coupled and the displacement-smoothing step.
A second highly successful scheme for the treatment of vacuum or air in this context is the use of nonlocal algebraic constraints that bind the displacement in the vacuum domain to the motion of the boundaries of the embedded bodies (Psarra et al., 2019).This scheme is in fact monolithic and thus does not suffer from the convergence issues of the staggered approach discussed before.However, the non-local constraints can be difficult to setup for complex geometries (Dorn et al., 2021) and furthermore lead to an unfavorable sparsity structure of the system matrix.Both, the staggered and the non-local constraint approaches -if applicable -lead to quantitatively better results than the probably most straightforward method from an implementation perspective: assigning a small but sufficient artificial stiffness to the vacuum domain (Keip and Rambausek, 2016;Dorn et al., 2021).By that one ends up with a perturbed boundary value problem similar to the one employed in Section 3 for the demonstration of spurious magnetic forces.We will refer to this approach as "naive monolithic scheme".As we have seen in Section 3, spurious coupling may eventually lead to crashes.Thus, for the "naive monolithic scheme" we face the problem that the spurious coupling sets a lower bound on the (auxiliary) stiffness added.Since the accuracy of the "vacuum" approximation depends on the ratio of "solid" to "vacuum" stiffness, the lower bound on the auxiliary stiffness will at some point lead to significant deviations from the unperturbed problem and its solution.The "naive" approach can be significantly improved by adapting the stiffness added to the magnitude and the "profile" of the magnetic fields in the vacuum domain as demonstrated by Rambausek et al. (2022).However, this requires a certain knowledge of the problem or complicated heuristics to automatically assign optimal stiffness parameters.Moreover, such improvements only shift the limit but do not solve the underlying problem.
A somewhat distinct method is that of Liu et al. (2020) which is essentially based on magnetic forces and tractions derived from numerically computed fields and (localized and free) electric currents.By that, no forces emerge in non-magnetic domains.In order to adapt the deformation of the vacuum domain to that of the slender structure under consideration they thus may assign a much smaller -indeed negligibleartificial stiffness.However, their specific scheme for the computation of magnetic volume and surface force densities requires to compute derivatives of the magnetization, which are not directly available in standard FE simulations and thus prevents a monolithic approach.Instead, their solution procedure is a staggered one based on the alternate solution of the magnetostatic and the elasticity problem.Hence, one needs a certain number of iterations between the two subproblems in order to solve the coupled problem.
What we propose below has partially great similarities with the method of Liu et al. ( 2020) on a higher level in the sense that we focus on what is going on at the interface between magnetic and non-magnetic bodies.However, in contrast to (Liu et al., 2020), our schemes essentially rely on the same building blocks as the naive monolithic scheme.By that the proposed methods can be implemented with minimal adaptions to existing code.In particular, there is no need to compute forces emerging from currents nor anything else that would require a staggered algorithm.Thus, the new schemes can be easily linearized and will be implemented in a monolithic way.When starting from an implementation of the naive monolithic scheme, only minor adaptions of code are required.

The Maxwell traction approach
For the derivation of the first approach we consider the "mechanical" weak form corresponding to the first variation of the governing potentials with respect to deformation restricted to vacuum domains  a .We emphasize that at this stage, the only (co-)energy density considered in vacuum is the magnetostatic one.Thus, the only stress present is the vacuum Maxwell stress (30) and the "mechanical" weak forms for the energy and the co-energy versions coincide As discussed in the context of (30), σ MW , the vacuum Maxwell stress is a divergence-free field such that we may safely9 transform the equation above by means of the divergence theorem where MW is the first-Piola-Kirchhoff-type instance of the Maxwell stress.The important observation here is that, in the case the FE discretizations under consideration10 , the boundary integrals in (47) do contribute any discrete "interior" forces11 .As a consequence, no spurious coupling within  a can emerge from such boundary integrals whereas the volume integral in (47) indeed suffers from this numerical artifact.12This is somehow parallel to the use of magnetic surface forces by Liu et al. (2020) but approached from a quite different angle.Indeed, the jump of the normal components of the "magnetic" stresses corresponds to the magnetic surface tractions (see ( 35) and ( 35)), but here we only compute the respective contribution on the vacuum side.
In order to obtain a non-singular system matrix in FE simulations we equip the vacuum domain with a very soft elastic material.However, in contrast to the "naive" monolithic scheme, the artificial stiffness does not have to counterbalance the effects of spurious forces and is thus only limited from below by requirements of the linear solver and numerical stability, respectively.Summarizingly, the above considerations ask us to implement the right-hand-side of instead of the left-hand-side, where Ψ a is the auxiliary strain energy density.Finite element implementations of the boundary integral in (48) approach are briefly outline below.The actual implementations in Netgen/NGSolve are provided as supplementary material.
Direct implementation of the boundary integral.The most obvious finite element implementation of the "Maxwell traction" approach is the implementation of the boundary integral in ( 48) just as what it is -a boundary integral.The essential additional algorithmic ingredients are • a routine the gathers all element faces that belong to the boundary of the non-magnetic domain and of which the parent (volume) element belongs to the vacuum domain and • the evaluation of the outward-pointing unit normal on these faces.
For our implementation we choose Netgen/NGSolve which directly provides the outward-point unit normal vector and a symbolic language for the integration over element boundaries.However, as will be shown in the numerical examples in Section 4.4, the convergence for refined meshes of such an implementation is not satisfactory.We also point out that while the boundary integral can readily be linearized to obtain its contribution to the system matrix, this contribution is non-symmetric.
Implementation via discrete force omission.Alternative to the boundary integral implementation, one can exploit the fact that MW is expected to be divergence free in a non-magnetic domain, such that all resulting discrete "interior" forces can just be ignored.In contrast, the discrete forces corresponding to boundary deformation degrees of freedom balance discrete forces stemming from applied tractions or the body on the other side of an interface.Hence, instead of actually computing the boundary integral in strict sense, we may simply compute the volume integral in (48).This is what one would also do in the "naive" monolithic approach but ignoring all discrete forces resulting from Ψ vac, 0 = MW ∶ that correspond to deformation degrees of freedom (DoFs) in the interior of  a .The actual implementation thus only needs to be able to distinguish between discrete forces in the interior of the vacuum domain and its boundary.Depending on the finite element software, the deletion of forces might be done by simply omitting the respective contribution during assembly.Alternatively, one can assemble the residual vectors resulting from Ψ a and Ψ vac, 0 separately in a first step.In a second step, the "interior" forces from the latter residual vector are deleted such that only the boundary forces remain.The final residual vector is the obtained as the sum of the auxiliary residual vector and the fleshed-out vacuum residual.It is important to note that the omission or deletion of discrete force contribution entails the deletion of the corresponding contributions to the system matrix, rendering the latter non-symmetric.
Our actual implementation employs the latter method as this can be done without any low-level changes to the assembly routines of NGSolve.

The traction compensation approach for air-like media
The second approach proposed also relies on tractions but, besides that, builds upon a converse line of thinking.The basic idea is to assign an auxiliary finite strain energy Ψ a to the vacuum domain then compensate the effect of the stiffness added by subtracting the corresponding mechanical traction from the boundary  a .This means adding an "integral zero" to the functional under the condition that the corresponding stress field is divergence-free.We start from the weak form ( 46) and proceed as outlined to obtain with While the right-hand-side of ( 49) is the basis for our implementations, further manipulations provide additional insight into the method: From this equation one can see that, irrespective of the actual value of σ a , only the Maxwell stress is exerted on embedded bodies ("across" the boundaries of the vacuum domain).Thus, as the actual deformation in the interior of the vacuum domain does not affect its neighboring bodies, which allows for great freedom in the choice (form and magnitude) of Ψ a or σ a , respectively.A second aspect is that one effectively solves for div σ = div(σ MW + σ a ) = 0 which allows for σ a ≠ 0 such that the auxiliary stiffness can counterbalance any spurious forces corresponding to σ MW ≠ 0. In that case, the terms added in the derivation do not exactly add up to zero.The crucial aspect from an implementation point of view is the treatment of the boundary integral in ( 49).In the remainder of this subsection we present two possible finite element approaches for which the actual implementation in Netgen/NGSolve is provided as supplementary material.
Direct computation of the traction compensation integral.Similar to the Maxwell traction approach, we start with the direct implementation of (49).The volume integrals are exactly the same as for the "naive" monolithic implementation.The only addition in this sense is the boundary integral in that equation, of which the implementation is technically parallel to that for the Maxwell tractions.It also has the same drawback of bad convergence under mesh refinement as will be shown in the numerical examples in Section 4.4.Similarly, the resulting system matrix is non-symmetric.

Implementation via discrete force omission.
Reusing the viewpoint of force omission approach for Maxwell tractions, one can implement traction compensation via manipulations of the discrete force vector resulting from standard evaluation of volume integrals and the corresponding system matrix contributions.However, in contrast to the Maxwell tractions, now one leaves the discrete forces corresponding to displacements in the interior as they are but removes the auxiliary discrete forces that correspond to displacement DoFs at the boundary of the vacuum domain.This means the we may think of implementing (49) as or Thus, for any implementation the essential algorithmic ingredient is a mean to distinguish "interior" and boundary" displacement DoFs.Our own implemention first assembles the volume integrals in (53) separately, removes the boundary forces from the auxiliary residual vector and finally adds both vectors.Accounting for the manipulations of discrete forces in the system matrix leads to a loss of symmetry of the latter, as is the case for all approaches proposed in this contribution.

Assessment of the "Maxwell traction" and the "traction compensation approach for air-like nonmagnetic domains
We critically assess the proposed schemes by their performance in the boundary value problem presented in Section 3, in particular Figure 4a.The energy densities employed are of the form (43).The material parameters for the magnetic domain are { , ′ , } = {1 kPa, 50 kPa, 10}.We consider the non-magnetic domain to be air-like, where = 0 and = ′ = 0.However, depending on the numerical scheme applied, we employ different values for the auxiliary hyperelastic parameters a and ′ a .While the effect of their absolute value is to be studied, we keep the ratio of a ∕ ′ a = 1∕50 throughout all examples.Concerning the discretization we employ triangular meshes and second-order finite element spaces for each component of the deformation and for the scalar magnetic potential.This appears to be a common choice in literature such that the results reported in this work are directly relevant to research applications.
From a theoretical perspective, the "traction compensation" approach is exact whereas the "Maxwell traction" approach presented is not.This is because in the latter case the effect of the auxiliary stiffness is very small but nevertheless finite.In contrast, all effects from the auxiliary stiffness and body force contributions are compensated in the former case.However, despite being exact in the above sense, the results obtained with the traction compensation approach still depend on the magnitude of added stiffness and forces.This is particular observable for rather small added stiffness such that the spurious coupling may still affect the deformation of the non-magnetic domains.Therefore, in what follows, we investigate the effect of the absolute value of the auxiliary shear modules a and the discrete resolution of the domain, i.e. the mesh size.
We start the comparison with individual parameter studies of the "naive" monolithic approach, continue with the "Maxwell traction" and "traction compensation" implementations and finally present a comparative convergence study.All results presented below refer to the vertical displacement 2 of the point with initial position = (0, ), i.e. the intersection of the boundary of the magnetic domain with the -axis, at an applied magnetic field magnitude ∞ = 0.7 T.
The coarsest mesh generated with Netgen is depicted in Figure 6a, a comparison of the coarsest (refinement level zero) and the finest (refinement level four) mesh is provided in Figure 6b.Level zero corresponds to a total of 222 cells and 1473 DoFs.After four refinement steps one has 56 832 cells and 343 203 DoFs.

The "naive" fully-coupled monolithic auxiliary stiffness approach
For some further insight on the underlying problem of spurious magnetic forces we have a closer look at the capabilities of the "naive" fully-coupled monolithic auxiliary stiffness approach.Interestingly, as can be seen from Figure 7, the approach works quite well for very fine meshes (refinement level four (RL 4) corresponds to more than 1 × 10 5 DoFs) but not so for coarse meshes.This underlines the fact that spurious magnetic forces result from inaccurate solutions of the magnetostatic part of the problem, which has already been put forward in Section 3. Figure 7 furthermore shows that -if results can be obtained at all -convergence with respect to mesh refinement is only a minor issue compared to the convergence with respect to the auxiliary stiffness parameter a .From a practical point of view, one might say "acceptable" results were obtained for already for a ≤ 1 × 10 −2 kPa and one refinement step."Accurate" results were only obtained a ≤ 1 × 10 −3 kPa requiring at least two mesh refinement steps, i.e. going from roughly 1000 to 20 000 DoFs.A direct comparison with the methods proposed in this work is presented in Subsection 4.4.6.

Maxwell traction approach via direct boundary integral computation
Figure 8a shows the effect of the auxiliary stiffness parameter a for different mesh refinement levels whereas subplot (b) shows the effect of the refinement level for given a .From the plots one can see that accurate results are only obtained for a ≤ 1e−3 kPa as already observed for the "naive monolithic" approach.Concerning the mesh resolution, three refinement steps (86 355 DoFs) yield a high-quality solution.The displacement is plotted over (a) the auxiliary free-space stiffness a and (b) over the mesh refinement level.The results for a = 1 × 10 −4 and 1 × 10 −6 kPa practically coincide such the graph for the latter is omitted in subplot (b).

Maxwell traction implementation via discrete force omission
The effect of the auxiliary stiffness parameter a is essentially the same as in Figure 8a.However, as shown in the mesh convergence plot depicted in Figure 9, the range of displacements is much more narrow than for the boundary integral implementation.This indicates a much higher accuracy of the volume-integral-based "discrete force omission" implementation compared to the "direct boundary integral" version.Also, convergence with respect to mesh refinement seems to be attained earlier with the force omission implementation.A study of actual convergence rates is deferred to Section 4.4.6.

Traction compensation implementation via direct boundary integral computation
In this case the value of the auxiliary stiffness parameter a is responsible for preventing excessive deformation due to spurious magnetic forces.Therefore, we employ considerably higher values than for the "Maxwell traction" implementations.Figure 10a shows a pronounced effect of a that clearly depends on the mesh resolution.Figure 10b indicated that the value of the auxiliary stiffness a is of great relevance for coarse meshes but looses its influence for increasingly fine discretizations.For example, with a = 0.1 kPa two refinement steps yield approximately equal results as four refinement steps and a = 1.0 kPa.Thus, Figure 10: Scaled displacement 2 (0, )∕ at ∞ = 0.7 T obtained with the "Traction compensation" approach implemented via direct boundary integral computation.Subplot (a) depicts the displacement plotter over the auxiliary stiffness parameter a whereas subplot (b) shows the displacement plotter over the mesh refinement level.Note that the parameter a plays different roles in the "Traction compensation" and "Maxwell traction" such that a comparison of its effect cannot be made directly.
lower auxiliary stiffness tends to lead to a significantly higher accuracy which is somewhat unexpected from the theory.

Traction compensation implementation via discrete force omission
As can be seen from Figure 11a, the force omission implementation of the traction compensation approach shows almost no sensitivity with respect to the auxiliary stiffness parameter a as expected theoretically.Only for the coarsest mesh (refinement level zero), on can see a small effect in Figure 11b.This is due to finite deformations caused by spurious magnetic forces in the case of a = 0.1 kPa.For finer meshes, the spurious forces decline such that the solutions for a = 0.1 kPa and a = 1 kPa practically coincide.
Also, from the tick-labels of the -axis in both plots one can see that the volume integral implementation overall is much more accurate than the implementation based direct boundary integral computation.

Comparative convergence study
For the convergence study below, we consider as reference the solution obtained with a staggered scheme along the lines of Pelteret et al. (2016).In order to enhance the accuracy of their scheme, we keep alternating between the coupled and the free-space adaption sub-steps until the overall solution is converged.To avoid visual clutter, we only consider one representative set of parameters for each method.Figure 12: Comparative convergence study in terms of (a) the scaled displacement 2 (0, )∕ plotted over the refinement level and (b) the relative displacement error over the number of DoFs in logarithmic axes.The data points in (b) correspond to refinement levels zero to three.The reference solution for (b) is that obtained with a staggered scheme for mesh refinement level four, i.e. 343 203 degrees of freedom.The shorthand legend labels expand to: "NV" -"naive" with a = 1e−3 kPa, "MW (B|V)" -"Maxwell traction" via direct Boundary integrals or Volume integrals and force omission with a = 1e−6 kPa, "TC (B|V)" -"Traction compensation" via direct Boundary integral or Volume integral and force omission with a = 1 kPa and "ST" -staggered scheme.
plementations are much more accurate than their boundary integral counterparts.Indeed, the former practically coincide with the solution obtained with the staggered scheme.The plot in Figure 12b depicts the convergence of the methods under consideration in terms of the deviation from the solution of the staggered scheme at refinement level four (343 203 DoFs) in dependence of the number of DoFs.The solid gray line corresponds to linear convergence and the solid black line to quadratic convergence in mesh resolution.One observes that the volume integral implementations are much more accurate and furthermore converge at super-linear rate.In contrast, the boundary integral implementations converge at most at a linear rate.Different from all other graphs, the one for the "naive" approach does not exhibit any significant convergence because its solution (or error, respectively) is dominated by the perturbation of the actual problem through the auxiliary stiffness.As a small, but nevertheless noteworthy detail we point out that only for very fine meshes, one can see a small deviation of the "Maxwell traction" implementation via force omission ("MW (B)") from the staggered and the "traction compensation" via force omission ("TC (B)").This stems from the auxiliary stiffness a = 1e−6 kPa that perturbs the solution in the case of the "Maxwell traction" approach.Due to the superior performance of the force omission implementations over their boundary integral counterparts, the latter will be discarded from further considerations.

Eliminating spurious coupling in non-magnetic solid domains
In the case of non-magnetic solid domains that are sufficiently soft to suffer from spurious magnetic forces, neither the "staggered" (Pelteret et al., 2016) nor the non-local constraint Psarra et al. (2019) approach can be applied, because both of the methods rely on the absence of physical stiffness.In contrast, both the "Maxwell traction" and "Traction compensation" approach can be extended to non-magnetic solids.
The adaption of the "Maxwell traction" approach and its implementations to this scenario is straight forward.The auxiliary strain energy density from before is simply substituted by the actual mechanical material model.For the "traction compensation" approach, the modification for solid domains is more delicate and requires a generalization of the fundamental idea: increase mechanical stiffness and compensate this increase by an appropriate increase of loads.Consider the weak balance of momentum for prescribed tractions and body force densities where the (total) Cauchy stress σ consists of the Maxwell stress, the mechanical stress resulting from the actual mechanical material model.In order to add a mechanical zero holding down spurious magnetomechanical interactions, we need to be more careful than before in order to not change the solution of the original problem.First recall that, within a non-magnetic material, the presence of magnetic field does not affect the mechanics.Thus, we may factor out anything related to the Maxwell stress, i.e. and whereby ( 56) is automatically fulfilled for any solution of the magnetostatic problem, since the medium under consideration does not experience any magnetic forces.Next, as the purely mechanical part (55) can be treated separately, we may scale it, i.e. multiply by some constant factor (1 + ), without affecting the solution to this equation.For the purpose of a more compact notation we denote any quantity multiplied by the compensation factor with a bar, i.e. σmech = σ mech .Then, and any (deformation or displacement) solution to (55) is a solution to (57) and vice versa.Doing the usual integration by parts followed by application of the divergence theorem yields The first set of terms are contained in the original boundary value problem whereas the second set of terms shall be employed to suppress spurious interactions.In order to proceed towards an implementation that is again based on omission of discrete forces obtained by volume integrals, we have a close look at the second "zero" in (58).The volume integrals can be implemented as they are, but instead of equating them with the boundary integral exploiting σ ⋅ = ̄ we again just omit the resulting discrete forces corresponding to boundary (interface) deformation, i.e. the additional terms to be considered are or In contrast to the traction compensation via discrete force omission for air-like media, here we employ a multiple of the mechanical strain energy density and a multiple of the body forces which must be scaled by the same factor.Applied tractions, however, do not need any special treatment.
4.6.Assessment of the "Maxwell traction" and the "traction compensation approach for non-magnetic solid domains In this subsection we build upon the trust in the force omission implementations of the "Maxwell traction" and the "traction compensation" approach gained in Subsection 4.4.The direct boundary integral implementations will not be considered due to their inferior accuracy.Now, in the case of non-magnetic solids, we are not any longer in possession of "reference" solutions obtained with the staggered scheme.Instead, we compare the solutions of the proposed methods to results obtained with the "naive" monolithic scheme, which admits direct extension to non-magnetic solids and is quite accurate and robust as long as the non-magnetic solid is sufficiently stiff.In fact, what we need to demonstrate in this section is not the accuracy of the new methods for very soft non-magnetic media as this has already been covered in the previous section.Instead, we want to demonstrate the correctness of the extension to actual solids.Therefore, we may choose material parameters such that a comparison with the "naive" scheme is indeed reasonable.

A magnetic disk in a non-magnetic carrier under gravitation-type and magnetic loading
Here we in essence reuse the boundary value problem from Sections 3 and 4.4 but with a specific energy density assigned to the non-magnetic domain.For simplicity, we again choose an energy density of the form (43) with parameters { , ′ , } = {1 kPa, 50 kPa, 10} for the magnetic domain and { , ′ , } = {0.5 kPa, 25 kPa, 0} for the non-magnetic domain.In addition, we assign mass densities13 to both the magnetic and the non-magnetic solid leading to gravitation-like forces.To be specific, the mass density in the magnetic domain is where denotes the gravitational loading parameter. 14In Figure 13a we show the scaled vertical displacement 2 (0, )∕ at = 1 m s −1 and ∞ = 0, i.e. purely gravitational loading.The corresponding mesh convergence plot is shown in Figure 13b.In both subplots one can observe the perfect agreement of all three methods.The picture is similar in Figure 14, where the scaled vertical displacement 2 (0, )∕ is depicted "naive" "Maxwell traction" (vol.int.) "Traction compensation" (vol.int.) at the magneto-mechanical loading state = 1 m s −1 and ∞ = 1 T. The deformed magnetic bodies at purely gravitational and combined loading are depicted in subplots (a) and (b), respectively, of Figure 15.There one can see that displacement of point "A" is downwards under gravitational load (subplot a) and upwards under combined loading (subplot b).

A bilayer beam under gravity air magnetic field embedded in an air-like domain
As final example that demonstrates the capabilities of the proposed schemes we consider a beam consisting of a magnetic and non-magnetic layer that is embedded in an air-like domain as depicted in Figure 16.It depicts a bilayer beam consisting of a magnetic  magn and a non-mangetic layer  nonm that is clamped at its vertical symmetry axis and surrounded by "empty" space  a .This setting is inspired by prior examples of MRE beams (Zhao et al., 2019;Mukherjee et al., 2021;Rambausek et al., 2022;Moreno-Mateos et al., (a)  2022b) and the prospective application of magnetoactive materials as mechanically active substrate for biological experiments (Gonzalez-Rico et al., 2021).In the simulations, the beam is first exposed to gravity = 9.81 m s −2 and in a second step to a combined loading through gravity and a uniform external magnetic field ∞ = 1 T. The generic energy density function employed is given as It has the same purely mechanical neo-Hookean type contribution as in examples.The magnetic part, however, now models saturation ‖ ‖ → s for ‖ ‖ → ∞.
The material parameters for this example are collected in Table 1.They have been chosen with great  magn 2 × 10 3 100 × 10 3 2 × 10 3 10 1  nonm 1 50 1 × 10 3 0 care to render a physically reasonable example that challenges the numerical methods under evaluation.Furthermore, when the "Maxwell traction" approach is applied, the air-like domain  a is equipped with a = 1e−6 kPa.In the case of "traction compensation", both the non-magnetic and air-like domains have a = 1e2 kPa, which amounts to a compensation factor of = 100 in the solid domain  nonm .This means Figure 16: Boundary value problem of a beam that consists of a magnetic  magn and a non-magnetic layer  nonm .It is clamped at its vertical symmetry axis and surrounded by "empty" space  a .The bilayer beam is exposed to gravity and a uniform external magnetic field ∞ .
that the auxiliary stiffness added to the solid is much higher that the actual stiffness.This deliberate choice helps us to confirm that all effects from the additional stiffness are compensated properly.
In contrast to all foregoing examples, the geometry of the bilayer beam has sharp corners.This complicates convergence studies because the "naive" scheme has severe problems in such cases.The reason are magnetic field concentrations near the corners leading to pronounced spurious magnetic forces that cannot be simply alleviated with mesh refinement (Rambausek, 2020, Chapter 9, Section 3.1.3).Therefore, we omit a rigorous convergence study and instead show contours of solution states in observe a gravity-driven bending-dominated deformation of the magnetic layer  magn and a slightly more general deformation pattern of the non-magnetic layer  nonm .The color contours in (a) correspond to the deformation magnitude.Subplot (b) shows the deformed configurations under combined loading.The deformation modes are very similar to those under purely gravitational loading, the displacements magnitudes are further increased through the applied magnetic field.The color contours in (b) correspond to the magnetization magnitude.They first of all confirm that only the magnetic layer exhibits magnetization as expected and, second, that the magnetization distribution is far from being perfectly uniform or trivial.This underlines the importance of full-field simulations of such boundary value problems.We highlight that in both subplots of Figure 4.6.2one can observe severe deformations of the "empty" domain surrounding the beam in the vicinity of the tip of the beam.While this was not of a great concerns here, in practice one could "smoothen" the deformation field by carefully increasing the auxiliary stiffness in such regions relative to that in other areas of the "empty" domain.Please also note that the post-processed deformation is only of first order but the actual deformation is the simulation is of second order, such that severely distorted elements in the vicinity of the beam tip are not displayed exactly how they actually appear in the simulations.In any case, the solutions obtained for the Maxwell traction and the traction compensation scheme coincide almost perfectly.
In this context we highlight that we did not need to fine-tune the auxiliary stiffness parameters to achieve excellent agreement between both approaches.This underlines not only their computational robustness but also their robustness with respect to the methods parameters.In general, we for the Maxwell traction scheme recommend to choose the auxiliary stiffness parameter in a range of 10 −6 of the softest solid parameter.For the "Traction compensation" we recommend stiffness parameters roughly in the range of those employed for the magnetic solids.

Conclusion
This work is centered around the issue of spurious magneto-mechanical interactions that is pervasive in "naive" fully-coupled numerical simulations of deformable magnetic bodies such as MREs.The key contributions are, first, a thorough characterization of the underlying issue.Second, we present two novel approaches that effectively eliminate or suppress spurious coupling in both vacuum-or air-like media and non-magnetic solid.The first scheme relies on the so-called "Maxwell tractions" that removes all unwanted magneto-mechanical interaction from the interior of the non-magnetic domain.For definiteness of the solution in "empty" (air-like or vacuum domains) the method, only needs a negligibly small auxiliary stiffness even in comparison with very soft solids.The second scheme is somewhat dual to the first in that it employs a sufficiently large auxiliary mechanical stiffness to suppress unphysical magneto-mechanically driven deformation in non-magnetic domains.The additional stiffness is balanced by additional body force contributions and a removal or a compensation of the resulting tractions on the interface to neighboring bodies.The common advantages of our proposed methods over existing successful schemes, in particular in comparison with the staggered schemes (Pelteret et al., 2016;Liu et al., 2020) and monolithic schemes based on non-local constraints (Psarra et al., 2019;Rambausek et al., 2022) are twofold.The first is the ease of implementation atop of "naive" monolithic FE simulations because, as has been shown, the successful volume-integral-based versions rely on only slight modifications of the weak form.There is no need for modifications to the overall solution procedure, nor cumbersome preprocessing and modifications to the sparsity structure of the linearized system matrix.Second, as demonstrated successfully, both of our approaches are not only applicable for air-like environments but also for actual, very soft non-magnetic solids.Another advantage over the staggered scheme is that both of the proposed approaches directly allow for the consistent linearization of the fully coupled system, which positively affects the convergence of nonlinear solvers.A mild disadvantage is that the resulting linear systems are non-symmetric, which increases the effort required for linear solves.To our experience, this is usually outweighed by the improved convergence and robustness of the methods.
We have assessed implementations of both proposed methods regarding their accuracy and their effectiveness regarding the elimination of spurious coupling in comparison with existing approaches.In particular, we have shown the convergence of the "Maxwell traction" approach for decreasing auxiliary stiffness, while the "traction compensation" method converges for increasing auxiliary stiffness.In this context is impor-tant to note that while there remains a parameter to be set, namely the auxiliary stiffness in one form or another, there is typically no need for tuning this parameter with great care.For the "Maxwell traction" one may go right in the range of 10 −6 of the softest solid parameter, for the "Traction compensation" one may simply employ some stiffness in the range of the magnetic solids under consideration.The critical comparison with existing schemes has demonstrated the competitiveness of the force-omission-based variants of the approaches proposed, with small advantages for the "traction compensation" approach.Thus, in view of the ease of computer implementation and simple choice of parameters both of the proposed methods can hope for wide adoption in future numerical investigations on MREs and related materials or problems, e.g. in electromechanics.

Figure 3 :
Figure3: A magnetized body  surrounded by a non-magnetic medium.The closed surface  lies in the non-magnetic domain and encloses the magnetic body.

Figure 4 :
Figure4: Test problem demonstrating a mild case of spurious deformation in a non-magnetic deformable elastic domain.Subplot (a) depicts the (quarter) boundary value problem of a rigid magnetic particle  magn embedded in a non-magnetic elastic medium  nonm exposed to a uniform external magnetic field ∞ .Subplot (b) show the result obtained with the magneto-mechanical energy formulation (23) with color-contours indicating to the magnetic field magnitude.The blue mesh corresponds to the deformed configurations, whereas the grey lines indicated the undeformed mesh.The dashed red triangles highlight parts of the undeformed configuration which are significantly deformed under applied field.The corresponding plot for the energy-co-energy formulation (24) is depicted in subplot (c).

Figure 6 :
Figure 6: Coarsest and finest mesh of the circular magnetic domain embedded in a non-magnetic domain.Subplot (a) depicts the "full" (actually quarter domain) mesh.Details of the coarsest and the finest mesh employed are shown in (b).Note that actual geometry representation is of second order, which has only been lost during postprocessing.

Figure 7 :Figure 8 :
Figure7: Scaled displacement 2 (0, )∕ at ∞ = 0.7 T obtained with the "naive" monolithic auxiliary stiffness approach.The displacement is plotted over (a) the auxiliary free-space stiffness a and (b) over the mesh refinement level.The "incomplete" graphs in (a) and (b) reflect the fact that the required auxiliary stiffness decreases with mesh refinement.

Figure 11 :
Figure11: Scaled displacement 2 (0, )∕ at ∞ = 0.7 T obtained with the "Traction compensation" approach implemented via volume integrals.Subplot (a) depicts the displacement plotter over the auxiliary stiffness parameter a whereas subplot (b) shows the displacement plotter over the mesh refinement level.The results for the displacement are shown in Figure12a.The plot confirms that the volume integral im- m −3 whereas in the non-magnetic domain we employ nonm 0 = 100 kg m −3 .Gravity points in negative vertical ( ) direction such that the resulting body forces are given as

Figure 13 :
Figure13: Comparative convergence study for gravitational load ( = 9.81 m s −2 , ∞ = 0) in terms of (a) the scaled displacement 2 (0, )∕ plotted over the refinement level and (b) the relative displacement error over the number of DoFs in logarithmic axes.The data points in (b) correspond to refinement levels zero to three.The reference solution for (b) is that obtained with the "Traction compensation" scheme at mesh refinement level four, i.e. 343 203 DoFs.

Figure 15 :.
Figure 15: Deformed configurations under (a) purely gravitational load and (b) combined gravitational and magnetic loads.The solid red line indicates the deformed magnetic domain  magn whereas the dashed red line indicates the boundary of the undeformed magnetic domain  magn 0 .The contours refer to the displacement magnitude, whereby the colors are adjusted to displacement range in (b).Please note that while the mesh nodes are connected by straight lines in the figure, simulations actually have been carried out with second-order geometry descriptions.

Table 1 :
Material parameters for the solid domains in the bilayer beam example Domain