Field-induced interactions in magneto-active elastomers - A comparison of experiments and simulations

In this contribution, field-induced interactions of magnetizable particles embedded into a soft elastomer matrix are analyzed with regard to the resulting mechanical deformations. By comparing experiments for two-, three- and four-particle systems with the results of finite element simulations, a fully coupled continuum model for magneto-active elastomers is validated with the help of real data for the first time. The model under consideration permits the investigation of magneto-active elastomers with arbitrary particle distances, shapes and volume fractions as well as magnetic and mechanical properties of the individual constituents. It thus represents a basis for future studies on more complex, realistic systems. Our results show a very good agreement between experiments and numerical simulations—the deformation behavior of all systems is captured by the model qualitatively as well as quantitatively. Within a sensitivity analysis, the influence of the initial particle positions on the systems’ response is examined. Furthermore, a comparison of the full three-dimensional model with the often used, simplified two-dimensional approach shows the typical overestimation of resulting interactions in magneto-active elastomers.


Introduction
Field-controllable materials facilitate a variety of engineering applications and are therefore of increasing interest in Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. today's materials sciences. Their particular advantage over other material classes is the potential to tune their macroscopic properties by different external fields. If, for example, micron-sized, magnetizable particles are embedded into a soft elastomer matrix, magnetic fields can be used as an external stimulus that either changes the stiffness of the compound material or induces a macroscopic deformation. Based on these effects, such magneto-active elastomers (MAEs) have recently been proposed for applications as actuators and sensors [1][2][3][4], tunable valves [5] and vibration absorbers [6,7], medical robots [8] as well as prosthetic and orthotic devices with controllable stiffness [9].
The selection of a microscopic continuum approach for the modeling of MAEs allows to account for their microstructure and, with that, identify and gain insight into the underlying mechanisms that influence the materials' macroscopic behavior: local magnetic and mechanical fields are resolved for systems with arbitrary particle distances, shapes and volume fractions. Moreover, this strategy offers the possibility to be easily adaptable for materials with different dissipative and non-dissipative properties of the individual constituents as, e.g. shown in Kalina et al [21] for the modeling of MAEs with magnetic hystereses. While the constitutive behavior of the materials' components is normally known from experiments and used as input data, the model still has to be validated with regard to its capability to reproduce the compound behavior. To the authors' knowledge, none of the models presented in the literature have been checked against real data so far since experiments, see e.g. [33][34][35][36][37][38], are performed with complex MAE samples which involve a lot of variables such as the exact particle shapes and their distribution within the system.
Within this contribution, a different approach is conducted: simulation results obtained by using the microscopically motivated continuum model that is proposed in Metsch et al [22] are compared to experiments for simplified two-, threeand four-particle systems. This procedure facilitates a detailed analysis of the samples' behavior with a manageable number of influencing factors-it makes a validation of the model possible. The applied experimental configuration is similar to the one used in the work of Puljiz et al [39]: the samples are exposed to a constant external magnetic field with varying angle-the motion of the particles is tracked with a CCD camera. However, in this work, the distance of the magnetizable particles is decreased. Thus, in contrast to particle interaction models, the effects of an inhomogeneous, non-linear magnetization can be accounted for. Additionally, the necessity for a finite deformation framework is highlighted.
The organization of the paper is as follows: in section 2, the experimental configuration is described. Details on the materials under consideration and the measuring procedures are given. The applied modeling framework is briefly summarized in section 3 while the model validation is carried out in the subsequent section 4: experimental and simulation results are compared for different systems using the full three-dimensional and widely used two-dimensional modeling approaches. Additionally, a sensitivity analysis regarding the influence of the initial particle positions within the sample is performed. After a discussion of the results, the paper is closed by concluding remarks and an outlook to necessary future work in section 5.

Experimental configuration
In order to perform a validation of the modeling approach presented in section 3, appropriate simplified MAE samples have to be produced and analyzed with regard to their deformation behavior in an external magnetic field. The associated experimental procedures are specified in the subsequent subsection 2.1. Moreover, the magnetization behavior of the embedded particles serves as an important input parameter for the numerical simulations. Its characterization is described in subsection 2.2.

Deformation behavior of simplified MAE samples
For the sample preparation, different requirements have to be met: all particles must be placed in the same layer, to allow a tracking of their displacements with in-plane measurements. Moreover, the particle shape has to be close to spherical and their distances must be sufficiently small, i.e. in the range of one particle diameter or lower, in order to capture effects of an inhomogeneous, non-linear magnetization. Concurrently, the stiffness of the surrounding elastomer material must prevent the particles from colliding since such effects, see the works of Puljiz et al [40] as well as Biller et al [15,41] for details, are outside the scope of this contribution and impede the numerical simulations. Finally, the size of the samples has to be chosen in a way, that the external magnetic field can be assumed as homogeneous and well-defined: this is of particular importance for the application of magnetic boundary conditions within the numerical simulations [42].
In a first step of the sample preparation, a Polydimethylsiloxan (PDMS) layer with a thickness of 2mm and a quadratic cross section of 15mm length is cured at a temperature of 60 • C for 24h. It consists of three components that were purchased from Alfa Aesar and Gelest Inc., respectively: 19wt% of PDMS with molecular weight 770 is added to a mixture of 9,1wt% HMS-151 as crosslinking agent and 90.9wt% DMS-V25 as functionalized polymer. The crosslinking reaction is started with a platinum catalyst. After casting the first layer of PDMS, superparamagnetic nickel particles, also purchased from Alfa Aesar, with nearly identical size and sphericity are selected and placed onto the layer with tweezersusing a micro manipulator device, they are pushed into their final positions which are given in table 1. For all samples, the aim is to achieve comparable distances between the particles, see figure 1. The final preparation step includes the encapsulation of the particles with another PDMS layer of the same size and composition. It is allowed to cure sufficiently long in order to ensure a good connection between the layers. The applied procedure facilitates a positioning of the particles with defined distances in the center layer of the surrounding matrix. Rheology measurements of the bulk material indicate that the  Table 1. Particle centers and equivalent diameters for the two-, three-and four-particle systems shown in figure 1. For each sample, the origin of the coordinate system is located in the center of the upper left particle -due to the symmetry of the experimental configuration with respect to the out-of-plane direction, the z-coordinate is equal for all inclusions. An upper limit for the error of the optical measurements is given by ±4µm. initial shear modulus of the elastomer is in the range of 2 kPa. However, uncertainties in the sample generation and the lack of a full experimental characterization required it to be determined by means of a fitting algorithm during simulation. This procedure is in line with the previous work of Puljiz et al [39]. For the measurements, the 32-magnet Halbach array described by Huang et al in [43] is used to subject the samples to a homogeneous magnetic field with a flux density of 170mT. Bearings allow for a rotation of the array in order to change the direction of the field. To ensure reproducibility of the experiments, a preconditioning step is performed: changes of the sample behavior due to, e.g. the well-known Mullins effect are overcome by applying a full loading cycle prior to the measurements. After an identification of the initial particle positions which are required as an input for the numerical simulations, the samples are placed in the center of the Halbach array. Thus, the magnetic field is applied instantaneously. Using a step size of ∆φ = 5 • , the movement of the particles is then tracked with a Matrix Vision mvBlueCOUGAR-S CCD camera and a Mitutoyo lens with 10x magnification. This setup allows for

Magnetization behavior of the particles
In order to determine the magnetization behavior of the particles, a novel approach is pursued. Instead of putting powder of the material under consideration into a cylindrical sample and using a vibrating sample magnetometer (VSM) [39,45], the behavior of just one particle with high sphericity is analyzed in a superconducting quantum interference device (SQUID). By doing this, the influences of an inhomogeneous distribution of the local magnetic field and enclosed air within the sample are expected to be minimized.
The extraction of appropriate particles is carried out with the assistance of a digital optical microscope, Keyence VHX-100, with integrated zoom lens. Due to a large variance in the particles' sizes, see figure 2, adequate working space for a manual manipulation of the powder is necessary. For the extraction of individual particles, a glass tip is thermally formed into a manipulator and the tip diameter adapted to the size of the selected particle-a suitable, slightly adhesive material at the glass tip enables prehension.
For the measurements within the SQUID magnetometer, the particle is embedded in a resin ball that is cast into a silicone mould. The glass tip with the particle is pressed into a prepared, partially cross-linked resin hemisphere, so that it Additionally, the fit for the same material using data obtained with a vibrating sample magnetometer is shown for a comparison [40].
adheres to the cut surface. Using a microscope with transmitted and simultaneous incident light, the presence and position of the particle on the cut surface are checked. The final spherical shape is achieved by pouring resin into the silicone mould. Since the cross-linked resin does not allow any further optical measurements, the presence of the magnetic particle within the sample is checked again without direct contact using a Matesy CMOS-MagView magnetic field camera.
As the particle cannot be detected by its own magnetism, the weak magnetic field of a large-area permanent magnet is used to generate a homogeneous magnetic background in front of which the particle emerges by a change in the field. In order to obtain the materials magnetization behavior, the diamagnetic background of the resin, which can be identified when the particle is fully saturated, as well as the demagnetization field caused by the particles' spherical shape have to be subtracted from the experimental data. Along with a fit using the hyperbolic tangent function the result is shown in figure 3. The calculation of the saturation magnetization, the scaling parameter and the relative permeability which crucially influences the initial slope and, with that, the magnetization behavior within the linear regime, yields M s = 314.5 kA m −1 , α = 188.6 m/kA and µ r = 5.9, respectively. A comparison with the magnetization behavior that is found for the same material using a vibrating sample magnetometer in the work of Puljiz et al [40] reveals that, although the parameter M s differs only by approximately 5%, the relative permeabilities show a difference of nearly 140%, see figure 3.

Continuum approach
In order to compare the experimental data to numerical results, the modeling framework presented in Metsch et al [22] is applied. It is based on a microscopic continuum approach which fully resolves the fields within the individual samples. Due to negligible magneto-mechanical coupling effects of the constituents, a decomposition of the free energy into a magnetic contribution that only depends on the Lagrangian magnetic field H and a mechanical contribution which is a function of the right Cauchy-Green deformation tensor C is admissible, see [22]. For the mechanical part of the energy, the compressible Mooney-Rivlin model is applied. Herein, the quantities J,Ī C 1 andĪ C 2 denote the Jacobi determinant as well as the first two principle invariants of the isochoric part of C, see the works of Flory [46] and Metsch et al [22] for more details. The split of the energy into a deviatoric term connected to the shear modulus G = E/(2(1 + ν)) and a volumetric term that is linked to the bulk modulus K = E/(3(1 − 2ν)) is an extension of the compressible neo-Hookean material model which is proposed by Wriggers [47]. The dimensionless parameter p controls the energy contributions of the invariantsĪ C 1 andĪ C 2 , respectively. To complete the model, the magnetic part of the energy must still be specified: while ϱ 0 Ψ * mag (H) = 0 holds within the non-magnetizable matrix material, the isotropic relation is used for the embedded particles [19,22]. It is a function of the norm H of the magnetic field H and-in accordance with the experiments shown in the preceding section 2-yields equation (1) for the particles' magnetization behavior. The quantity µ 0 represents the permeability of free space.

Numerical simulations
In order to ensure a good comparability of experimental and numerical results, the full samples are considered within the finite element simulations. Assuming spherical particles with center points located in the same plane x 3 = const., see figure 4, a symmetry with respect to the out-of-plane direction is used to reduce the computational effort. A visualization of the resulting meshes using ParaView [48] is depicted in figure 4: to account for rapid changes of the fields under consideration and prevent volumetric locking effects within the quasi-incompressible matrix [49,50], a high resolution of the particles and their surroundings can be seen. The meshes are coarsened towards the boundaries of the samples. The loading in the simulations is implemented as follows: regarding the mechanical field, the displacement u is fixed on the lower and lateral boundaries of the samples while the symmetry condition u 3 = 0 is imposed on its top layer. Corresponding to the experimentally applied magnetic induc-tionB max = 170mT, the external magnetic field is increased  to a maximum valueH max = 135.3kA m −1 within 10 ramping increments. This is in contrast to the instantaneous application of the magnetic field within the experiments but, since no rate dependent processes are involved, results in the same initial state. During the following 36 increments, the angle φ of the magnetic field is changed in steps of ∆φ = 5 • . This results in a clockwise rotation of 180 • and allows for a pointwise comparison of experimental and numerical results. The magnetic loading is illustrated in figure 5.
In all simulations, a Youngs' modulus E p = 197 GPa, a Poissons' ratio ν p = 0.3 and a parameter p p = 1 are applied for the embedded nickel particles. In conjunction with equation (3), this choice yields the special case of a compressible neo-Hookean material law [47] which allows to capture effects of finite rotations but accounts for the fact that the particles are comparatively stiff. Moreover, the saturation magnetization M s and the scaling parameter α determined with the experiments shown in section 2 are used to describe the particles' magnetization behavior. The deformation behavior of the non-magnetizable elastomer matrix is captured with the full Mooney-Rivlin material law. To account for its quasiincompressibility, the Poissons' ratio ν m = 0.49 is chosen. Assuming equal contributions of the invariantsĪ C 1 andĪ C 2 to the deviatoric part of the energy, see equation (3), the parameter p m = 0.5 is used. Since the Youngs' modulus E m is not known a priori and will be different for the individually prepared samples, the finite element analyses are performed covering a range of different values. Using the normalized measure which identifies the error of the simulation regarding the change in the inter-particle distances ∆d ij , the optimal values E m are determined in a least-squares sense for each sample. In the final comparisons of the experimental and numerical results, these values are used to validate the model. All results presented in this section were obtained with finite element simulations using FEniCS 1 [51,52]; the associated meshes were generated with Gmsh [53].

Two-particle system
The results for the two-particle system are summarized in figure 6: in order to determine the optimal Youngs' modulus of the matrix for a comparison to the experimental data, a range 6600 Pa ≤ E M ≤ 8100 Pa is covered in steps of ∆E M = 150 Pa. In figure 6 (a), the trend of the squared error ε shows a distinct minimum for E M = 7244 Pa -it is found using a fit with a fourth-order polynomial. Since it was the aim of the fabrication process to work with a matrix material which is capable of preventing inclusions with a comparatively low initial distance from colliding, see the work of Puljiz et al [40] for a detailed analysis on such effects, the value seems reasonable.
In figure 6 (b), the change of the inter-particle distance ∆d 12 is shown for varying angles of the external magnetic field. Attraction of the particles can be observed over a wide range. Only in a small band, in this case 65 • ≤φ ≤ 125 • , repulsion occurs. This is in accordance with the findings of Biller et al [15]. The comparison of experimental and simulation data reveals a very good qualitative and quantitative agreement: the results almost coincide for all values of φ. A maximum discrepancy of around 2.5µm can be found between φ = 55 • and φ = 60 • which is close to the transition from attraction to repulsion. The accuracy of the simulation shows that the assumption of magnetizable particles with a spherical shape is Figure 6. Results for the two-particle system: (a) summed squared errors for simulations with different Youngs' moduli E M of the matrix and corresponding fit using a fourth order polynomial, and (b) comparison of experimental and simulation data for the change of the inter-particle distance ∆d 12 in a rotating magnetic field. The comparison is carried out for the optimal Youngs' modulus E M = 7244 Pa. an appropriate first approximation although the experimental configuration in figure 1 (a) illustrates a clear deviation from that.
The comparison of 2D and 3D simulations with the same Youngs' modulus, which is also shown in figure 6 (b), reveals a qualitatively good agreement of both approaches: attraction of the particles is found over a wide range of directions of the external magnetic field and also the transitions from attraction to repulsion occur at nearly identical angles φ. However, it is apparent that the two-dimensional approach overestimates the resulting change in the particle distances for a given magnetic field. This is in good accordance with the findings presented in Metsch et al [22]. The underlying assumption of a cylinder in the simplified two-dimensional setting yields differences in the particle-particle interactions due to changes in the particlevolume fraction and the magnetization. Moreover, the particlematrix interactions are different since, for cylinders, the matrix motion in the out-of-plane direction is hindered. Accordingly, the particle displacements within the plane must be increased.

Three-particle system
The procedure for the three-particle system is similar to the one explained in the preceding section: a range 6150 Pa ≤ E M ≤ 7650 Pa is covered in steps of ∆E M = 150 Pa in order to determine the optimal value for the Youngs' modulus of the matrix material in a least-squares sense. Since the behavior of MAEs is expected to be significantly influenced by the distribution of the magnetizable particles and-due to their nonsphericity, see figure 1 (b)-the specification of their center points can only be accurate to a certain level, this particular system is chosen to additionally analyze the sensitivity of the numerical results with respect to the initial particle positions. To this end, the position of the center particle 2, see figure 1 (b), within the finite element mesh is changed by ±5µm in the horizontal and vertical directions. In view of the system specifications given in table 1, this corresponds to a slight modification of only 2.7% of the particle diameter. The resulting configurations are summarized in table 2.
Regarding the original configuration, the summed squared error ε in figure 7 (a) shows a minimum for E M = 6842 Pa - Table 2. Sensitivity analysis regarding the influence of the initial particle position on the simulation results: updated coordinates of the center particle 2 for the considered configurations. The names of the additional configurations indicate a position change in the horizontal (H) or vertical (V) directions by plus (P) or minus (M) 5µm. this corresponds to a shear modulus of approximately G M = 2300 Pa and is therefore in good agreement with the experimental estimation. With a reached minimum of ε ≈ 60µm 2 , the accuracy of the numerical simulation is also comparable to the one of the two-particle system. A closer look on the configurations with modified initial particle distances reveals that a positive horizontal as well as a negative vertical displacement of the center particle decrease the error of the numerical simulation whereas displacements in the opposite directions tend to increase ε. The best result is achieved with the VM5 configuration: it reduces the minimum error level of the prediction to ε ≈ 32µm 2 while the optimal value for E M remains unchanged. Along with the results of the original setting, it is therefore included in the subsequent comparison to the experimental data. In figure 7(b), the results for different directions of the external magnetic field are shown. The comparison indicates a good qualitative and quantitative agreement for all inter-particle distances. While the numerical and experimental values nearly coincide for ∆d 13 and ∆d 23 , small deviations can be found for the change in the distance ∆d 12 . These could be related to a non-sphericity of the particles or a lack of symmetry with respect to the out-of-plane directiontheir elimination requires a consideration of the full sample with data gained using, e.g. microtomography measurements [54]. Regarding the VM5 configuration which is also shown in figure 7 (b), an improvement of the numerical simulation can especially be seen for ∆d 12 . This implies that the actual initial Results for the three-particle system: (a) summed squared errors for simulations with different Youngs' moduli E M of the matrix and corresponding fits using a fourth order polynomial. The individual curves represent results for different initial positions of the center particle according to table 2. In (b), the comparison of experimental and simulation data for the change of the inter-particle distances ∆d ij in a rotating magnetic field is carried out for the optimal Youngs' modulus E M = 6842 Pa. Additionally, the improvement of the solution for the VM5 configuration is shown-it is indicated with dashed lines. Figure 8. Results for the four-particle system: (a) summed squared errors for simulations with different Youngs' moduli E M of the matrix and corresponding fit using a fourth order polynomial, and (b) comparison of experimental and simulation data for the change of the inter-particle distances ∆d ij in a rotating magnetic field. The comparison is carried out for the optimal Youngs' modulus E M = 7265 Pa. position of the center particle 2 might slightly deviate from the experimentally determined one and underlines the sensitivity of MAE samples with respect to the particle distribution. However, the good agreement between experiment and simulation indicates that the applied modeling strategy allows to capture all relevant effects which govern the samples' behavior.

Four-particle system
The results for the last and most complex sample under consideration are summarized in figure 8. Using a range 6600 Pa ≤ E M ≤ 8100 Pa with steps of ∆E M = 150 Pa and a fourth-order polynomial for the determination of the optimal value E M of the matrix, figure 8 (a) shows a minimum of the summed squared error ε for E M = 7265 Pa. It is apparent that the reached minimum error level is not as low as for the preceding samples. Since every particle increases the probability that one of the underlying assumptions-spherical particles, symmetry with respect to the out-of-plane direction-is not fully met, this effect is unavoidable. However, with a minimum error level ε ≈ 82µm 2 , the accuracy of the simulation is still comparable to the one reached for the two-and three-particle systems.
The comparison of experiment and simulation by means of the change in the inter-particle distances is shown in figure  8 (b). Again, a very good qualitative and quantitative agreement can be found for all values ∆d ij over the whole range of directions of the external magnetic field. The fact that the maximum discrepancies occur in connection with particle 4 could indicate a slightly wrong initial position. Additionally, minor deviations can be seen for ∆d 23 in the range 125 • ≤φ ≤ 160 • : regarding the apparent non-sphericity of particle 3, see figure  1 (c), these errors can again only be eliminated with a detailed analysis of the real sample.

Concluding remarks
Within this work, field-induced interactions in simplified samples are analyzed in order to validate a microscopically motivated, fully coupled continuum model for MAEs with experimental data. To this end, samples comprising two, three and four particles are prepared and the effect of an external magnetic field with varying angle on the resulting mechanical deformations is tracked experimentally. Moreover, a novel approach to determine the particles' magnetization behavior is presented. This reduces uncertainties known from other methods and provides important data for the numerical simulations.
For all samples, a very good agreement between experiments and simulations is found: the change of the interparticle distances within the individual systems is captured in all situations although simplifying assumptions such as a sphericity of all particles and a symmetry of the systems with respect to the out-of-plane direction have been made. This shows that the applied modeling strategy represents an adequate approach to describe the microstructural interactions in MAEs-consequently, it is also applicable to the analysis of more complex, realistic systems. With an additional comparison between the full, three-dimensional simulations and frequently used simplified, two-dimensional approaches, the typical resulting overestimation of the latter is observed. This confirms findings of former studies [22] and points out that two-dimensional simulations are at best only applicable to get a rough impression of the MAEs qualitative behavior. Finally, the performed sensitivity analysis of the samples with respect to the initial particle positions highlights the crucial importance of a profound knowledge of the MAE microstructure in order to be able to provide accurate quantitative predictions of the sample behavior.
Future research must strive for better agreement between experimental and simulation results by analyzing the actual sample geometry. This can be achieved by incorporating data gained from, e.g. microtomography measurements [54] in the meshing procedure of the simulation. However, since realistic MAE specimens comprise multiple particles with complex shapes and distributions, numerical simulations of whole samples with a reasonable computational effort will remain impossible. For this reason, an appropriate computational homogenization framework has to be applied to deduce an accurate macroscopic model for MAEs from data obtained with simulations that are carried out for a representative section of the microstructure. To this end, an extension of the recent contribution by Kalina et al [55] to the general, threedimensional case appears to be promising.