High fidelity, discrete element method simulation of magnetorheological fluids using accurate particle size distributions in LIGGGHTS extended with mutual dipole method

We simulate magnetorheological fluids (MRF) using open source LIGGGHTS soft sphere discrete element method code, extended by us to include a mutual dipole magnetic model. Our simulations take advantage of the many pair forces available in the LIGGGHTS framework, including SJKR cohesion, friction, and rolling resistance. In addition, we have included an uncoupled, Couette flow background carrier fluid. The simulated particles in this work are polydisperse, with distributions made to match the distributions used to produce magnetorheological fluids in literature, increasing the fidelity of the simulations. Using the accurate particle size distributions, high heritage contact models, and an uncoupled fluid model, we are able to match experimental MRF yield stress results more closely than with monodisperse simulations.


Introduction
Magnetorheological fluids (MRF) are mixtures of paramagnetic particles, with sizes on the order of microns to 10s of microns, in suspension with a carrier fluid and small amounts of various additives [1][2][3][4]. An MRF flows similarly to a nearly Newtonian fluid when there is no applied magnetic field. When a magnetic field is applied, the paramagnetic particles form structures that impede the flow, resulting in a non-Newtonian, shear thinning, Bingham plastic (figure 1) behavior. The MRF acts like a solid below a critical yield stress, then flows with a linear shear rate-stress curve offset from the origin by the yield stress. When the magnetic field is removed, the MRF reverts to its nearly Newtonian behavior [5]. The yield stress of MRF increases with magnetic field strength until saturation, at which point increasing the magnetic field does not significantly increase the yield stress. The maximum yield stress for MRF is ∼100kPa, at saturation of ∼700−1000mT [6][7][8].
The properties of MRF depend on the particle size and packing fraction. As packing fraction increases (f), the off-state viscosity increases and the on-state yield stress increases [10,11]. Most MRF formulations tend to fall between f = -0.3 0.45 representing a good compromise between off-state viscosity and yield stress. Similarly, as particle size increases, the off-state viscosity and on-state yield stress increase. However, as the particle size gets too large, the particles quickly settle out of suspension. Typical MRF formulations tend to have particle sizes between 1-10μm. Interestingly, a bidisperse mixture of particles enhances the yield stress beyond that of a purely larger particle mixture, while also having an off-state viscosity lower than that of the smaller particles [12][13][14]. The lower off-state viscosity may be due to the reduced ratio of the packing fraction to the maximum packing fraction with bidisperse particle distributions but the mechanism for the yield stress enhancement is not well understood.
The unique viscosity characteristics of MRF make them interesting in a wide range of applications, including solid state valves, flow control, power transmission, polishing, and robotics [1][2][3][4]. For many applications, it is desirable to maximize the yield stress for an applied magnetic field, while minimizing the off-state viscosity.
Many of these properties are empirically determined. While empirical, bulk properties of MRF are sufficient to achieve good results in applications, better understanding of the MRF microstructure in use and how different formulations affect the characteristics can allow for better prediction of performance and optimization of designs, as well as development of novel applications.
In order to better understand MRF, it is useful to see what is happening at the microscopic level. Attempts have been made to capture the microstructures experimentally in both electrorheological and MR fluids by using resin as the carrier fluid and allowing it to set inside a field or making quasi-2D cells [5]. However, resin cannot capture dynamic or transient effects and quasi-2D cells affect the physics of the interactions. To overcome these shortcomings, many attempts have been made to simulate MRF using the Discrete Element Method (DEM) [15][16][17][18][19][20]. DEM takes a brute force strategy; summing forces on individual particles and integrating their equations of motion [21]. Because the particle dynamics are fully determined from integrating Newton's 2nd Law, additional physics are easily implemented during the force summation step.
Early simulations were limited to quasi-2D models due to computational limitations. Such models are subject to similar limitations as the experimental work with monolayer cells [16], but were able to recreate the  . Two log-normal distributions used to create bidisperse MRF mixtures from Trendler and Bose with mean particle sizes of approximately 1.8μm and 6.7μm for the small and large particles, respectively [12].
structures seen in the experimental cases. Mohebi extended simulations into full three-dimensional space showing the formation of three-dimensional particle chains, however the magnetic induction model assumed fixed dipoles unaffected by the fields of adjacent particles and there was no attempt to characterize the rheology of the fluid [15]. Keaveny developed a mutual dipole method where the magnetic fields from nearby particles are included to calculate the magnetic moment of each particle, iteratively converging to a more accurate value [22]. This model has been used in several works [17,20]. The effect of the hydrodynamics model on simulating MR fluids has also been examined. Often, uncoupled Stokes drag is assumed, where the carrier fluid is not affected by the particles and the drag on the particles is independent of its neighbors. In the case of a single grain at low Reynolds number, Stokes drag is very accurate; however, because of the dense packing of grains in contact with each other in an MRF, the assumptions break down. To get around the breakdown of the assumptions, attempts have been made at using coupled fluid dynamics where the effect of the particles on the fluid and the fluid on the particles is calculated using computational fluid dynamics (CFD) [19,20]. This increases the fidelity of the simulation; however, it also drastically increases computation time. Thus, for simple geometries with low shear rates, where relative velocities between the grains and the carrier fluid are low, uncoupled models are preferable for their lower computing costs.
One area that has not been fully explored is the effect of particle size distributions on DEM simulations of MRF. Most previous work has assumed a monotonic particle size distribution [15-17, 19, 20]. Kittipoomwong investigated bidisperse distributions, but again, used two monotonic particle distributions. Despite the simple particle distribution model, they were able to show a yield strength enhancement for the bidisperse particle distribution, though the absolute yield stress was considerably lower than experimental values [18]. Actual particle distributions tend to follow log-normal distributions (figure 2) [12,13]. Due to the cubic relationship between volume and diameter, this can result in mixtures that have many smaller particles. A theoretical lognormal distribution was examined by Sherman using a simple fixed dipole model [23]. In this work we examine the effect of using hi-fidelity contact models and reproductions of the particle size distributions found in experimental investigations of MRFs on the simulation of magnetorheological fluids using DEM software.

Model
For our simulations we use LIGGGHTS an open-source Soft Sphere Discrete Element Method software [21]. LIGGGHTS, which is built on very mature LAMMPS code, is highly parallelizable and uses a modular code structure that makes it suitable for extension with custom models. We use 3D shear cells to examine the full structure and behavior of the MRFs with high fidelity to experimental and applied geometries.

Contact model
The particle-particle contact forces govern how particles resist interpenetrating and rebound off each other in collisions. There are two common contact models used in DEM simulations with magnetic particles. The highest fidelity model is the Hertzian model, which treats the particles as an elastic material with a defined Young's modulus and Poisson ratio. Given two spherical particles with radii R R , i j , Young's modulus Y Y , i j , Poisson ratio n n , i j , and coefficient of restitution e, co-penetrating each other a distance Dn ij , with a normal component of their relative velocities v ⊥ , the normal force between particles is given by (equation (1)). The tangential force is similarly formulated with an additional sliding friction coefficient m x which limits the maximum Hertzian tangential force m  F F t n x [21].
where Y * , m * , and R * are the reduced Young's modulus, reduced mass, and reduced radius of the particle pairs respectively. This contact force has the advantage of creating realistic separations between particles in contact, which is important due to the d 4 dependence of the magnetic force, however the spring force ramps up very quickly, creating a very stiff system. The stiffness of the system requires very short timesteps in order to maintain numerical stability in the system. To reduce the stiffness, some simulations have implemented an exponential contact force model (equation (2)).
where Q is chosen to cancel the maximum attractive magnetic force between the particles and k is chosen to balance the tradeoff between the stiffness of the system with the additional, unphysical repulsion force that extends past the particle contact distance [17].
To maximize the fidelity of our simulations, we use a Hertzian contact model in all our simulations. We also include several additional standard LIGGGHTS contact forces. A sliding friction coefficient m x limits the maximum Hertzian tangential force m  F F t n x . Rolling resistance is included to mimic some of the mechanical   effects of some asphericity in the iron grains using the Elastic-Plastic Spring-Dashpot (EPSD) model. This creates a static restoring spring torque back to original contact point, which ramps up from zero to a maximum value as displacement increases, and then falls back along the same slope to zero as if the particles start to roll backwards (figure 3) [24]. Finally, we include a cohesion term using the SJKR cohesion model F c =k c A, where A is the contact area and k c is a volume energy density. Values used in our simulations can be found in table 1.

Magnetic induction
Paramagnetic grains with susceptibility χ, such as those in an MRF, become magnetized when in the presence of an external magnetic field H  . In the case of a sphere in a uniform magnetic field, the grain is uniformly magnetized, resulting in a dipole moment (equation (4)) [25].
The force on the resultant dipoles is found in equation (5).
Notably, there is no force on a magnetic dipole in the presence of a uniform magnetic field. There is a torque term in the magnetic dipole interaction; however, because we use soft magnetic grains, the dipole of the particles will always be aligned with the local magnetic field. The resulting torque will always be zero. Therefore, we omit the magnetic dipole torque equation in our simulation.
In disperse mixtures, it is often acceptable to use the dipole calculated from the background field in simulations (fixed dipole method). However, in the dense mixtures of MRFs, the magnetic fields from neighboring grains begin to have a significant effect. To capture the effect of neighboring grains, we implement a mutual dipole method [22]. In the mutual dipole method, all of the grains are magnetized using the background magnetic field. Then the dipoles are recalculated using the background field plus the magnetic fields created by all near neighbors. This calculation is then iterated until the dipole moments of the individual particles converge (equation (6)).
Increasing the distance neighbors are considered will increase the accuracy of the model while incurring computational costs increasing with the square of the neighbor length. An infinite neighbor-length recreates a 2nd order approximation (terms higher than dipole ignored) of the entire system. A neighbor length of three particle diameters was used by Sherman simulating non-disperse particle [26]. Han, Feng, and Owen demonstrate that using the mutual dipole method with a two-particle chain, parallel to the background field, increases attractive forces by 77.78% while two particles adjacent to each other, aligned perpendicular to the background field, have their repulsive force reduced by 20.99% [20]. This model converges rapidly, with forces on an L shaped chain of particles converging to within one part in 10 5 within five iterations. Due to the rapid convergence and extremely short timestep from the stiffness in the contact forces, we make the assumption: and only iterate the mutual dipole method once per timestep. This, by the assumption that the motion of a particle in one timestep results in a very small change in the local magnetic field, does not introduce appreciable error, and greatly reduces the computational cost of the mutual dipole method.

MRF simulation
In order to test the effects of particle size distribution on MRF simulations, we attempted to recreate the properties of the MRFs in Trendler and Bose (figure 5) [12]. Trendler examined the effect of size and bidisperse mixtures on MRF properties. They reported precise curves for their particle size distributions (figure 2) found through laser diffraction, allowing us to recreate their MRF mixtures in LIGGGHTS. Trendler and Bose do not specify whether their distributions are volume or number densities, so simulations are performed with both interpretations. The two distributions they used to make their mixtures have mean particle size of 1.8μm and 6.7μm that were mixed in varying proportions for their experiments. Due to computation time limitations, we only attempt to reproduce two mixtures from their experimental data; a 00:100 ratio and a 33:67 ratio of small to large particles with a packing fraction/vol% of f = 0.3. The parameters for the particle contact forces are given in table 1. Most of the properties are chosen as accepted values for lubricated iron. The Young's modulus notably varies from the accepted value for iron at approximately 1% of the accepted value. Reducing the modulus reduces stiffness and allows for an increased timestep. Chen et al found that for moduli down to 0.1% of the actual material modulus, the behavior of a granular mixing simulation is statistically unchanged [27]. Moreover, in our simulations, particle chains will be largely under tension, further reducing the dependence on the Young's modulus.
The carrier fluid dynamic viscosity is 10mPa·s. This matches the silicone oil used by Trendler and Bose as their carrier fluid.
Unless otherwise specified, all simulations are carried out in a shear cell 250μm×250μm×250μm, with periodic boundary conditions in the X and Y directions, and fixed boundary conditions in the Z direction.
The shear cell is initialized by filling the simulation volume with particles randomly packed with a particle distribution pulled from a digitized version of Trendler's data ( figure 2, table 2). A constant magnetic field is then applied to the volume in the Z direction, and the particles are allowed to reach an equilibrium state, comprised of particle chains aligned along the magnetic field. Once the particles have settled, the top and bottom 15μm are frozen to create particle rafts to act as the shearing surfaces. Prior work has used either densely packed, crystalline particle rafts, randomly distributed particle rafts, or monolayers as their shearing surfaces. This tends to create a weaker attachment of particle chains, which slip across the shearing surfaces. By using the magnetic structures, the surface of the raft has strong magnetic field concentrations, which create strong attachment points for the active volume particle chains and are more representative of physical shear cells.
Simulations were run with 600mT background fields for the two different mixtures with the two different particle distributions. With just the 00:100 mixture MRF, the following additional simulations were run: with a monodisperse particle distribution; with each dimension of the shear cell doubled; without cohesion; with a stationary carrier fluid; and with no rolling resistance. We also ran simulations with the large particle mixture  Figure 6(a) shows monodisperse 6.7μm particles forming vertical chains when the magnetic field is applied. As the cell is sheared, the disordered chains begin to form sheets with crystallographic domains. Figure 6(b) shows a lognormal, number density distribution. During the initial magnetization, the structures formed are thick chains tending towards the diameter of the larger particles, sometimes coming together to form large clumps. As the cell is sheared, the chains also shear over with the cell but remain morphologically similar to the unsheared chains. Figure 6(c) shows a lognormal, volume density distribution. The initial magnetization still produces structures with a thickness comparable to the larger particles despite there being relatively fewer large particles. The chains again shear along with the cell though some additional clumping occurs. and varying magnetic fields at shear rates down to 100 s −1 corresponding to the shear stress/magnetic field experiment performed by Trendler and Bose (data recreated in figure 5(b)).

Particle structures and morphology
Our simulations formed particle chains along the applied magnetic field after approximately 1e5 timesteps. In the case of the lognormal distributions, rope-like chains and thick walls are formed at approximately the diameter and width of the largest particles (figures 6(b) and (c)). These thick chains are stable under shear until they break and are reformed into new morphologically similar chains. In contrast, the monodisperse particles initially form one particle thick, disordered chains, which quickly coalesce into sheets composed of crystallographic domains ( figure 6(a)). The domains rotate under the shear, but do not break and reform as obviously as the chains formed by the lognormal distributions. Videos of shearing behavior for a Y cross section of the 00:100 mixtures can be found in Supplementary Videos 1, 2, and 3 available online at stacks.iop.org/ MRX/8/085701/mmedia.  Figure 7(a) shows monodisperse 6.7μm particles forming extremely convoluted, single particle wide structures when the magnetic field is applied. They then straighten into single particle wide lamellar structures, aligned with the direction of shear as the cells are sheared. Figure 7(b) shows a lognormal, number density distribution. During the initial magnetization, the structures formed are thicker tending towards the diameter of the larger particles. As the cell is sheared, there is some alignment along the direction of shear, though the formation of lamellar sheets is muted. Figure 7(c) shows a lognormal, volume density distribution. The initial magnetization still produces structures with a thickness comparable to the larger particles. There is little alignment into lamellar structures during shearing.
Observing the Z-plane cross-sections, it can be seen more clearly that the monodisperse particles form structures that are predominantly one particle wide, particularly after they have been sheared ( figure 7(a)). The structures also align themselves prominently along the direction of shear, with the crystallographic sheets seen in the Y-Plane cross-sections forming roughly parallel bands across the Y direction. The lognormal distributions show considerably less banding and continue to exhibit thickness similar to the largest particle diameter (figures 7(b) and (c)). The gaps between structures also appear as approximately the same width as the largest particles in the lognormal distributions and single particle wide in the monodisperse case. Videos of shearing behavior for a Z cross section of the 00:100 mixtures can be found in Supplementary Videos 4, 5, and 6.  8(a)) The number density distribution produces strands and clumps similar morphologically to those from the volume distribution from the 00:100 mixture. (figure 8(b)) breaks from the trend and has structures somewhat smaller than the largest particles with more cross linking from the smallest particles and more numerous and smaller gaps between strands. In the 33:67 distributions many of the same trends continue ( figure 8). With the addition of smaller particles, the particle chains become tighter and more numerous, with smaller and more numerous gaps particularly pronounced in the volume density distribution. The 33:67 mixture volume distribution is the first simulation where the particle chains are smaller than the largest particle diameters and instead tend towards the diameter of the more numerous, medium-large diameter particles. The volume density mixture also creates smaller and more convoluted structures than any of the other simulations with considerably smaller gaps between structures ( figure 9). There is virtually no plane formation along the shear direction in the case of the volume distribution. Video of the shearing behavior for the 33:67 number density distribution can be found in Supplementary Videos 7 and 8.
One interesting feature is the way that the chains behave at large strains. Past works have characterized the chains as reaching a critical strain where the tensile stress in the chains exceeds the magnetic forces, at which point the chains break and then reform with other broken chains with lower shear and strain [26]. In our simulations using particle distributions, the deformation and destruction of chains are considerably less distinct. During shear, the gaps between chains shrink. As the shear angle of the chains increases, the component of the background field across the chains increases. Neighboring chains then link together and the attachment points split into new, lower shear chains ( figure 10). Figure 10. The progression of a chain of particles (red, figure 10(a)) from a disjoint, low shear pillar to a higher shear chain, which then attaches itself to adjacent chains (green, figure 10(b)). The attachment points then go on to form sections of a new set of chains ( figure 10(c)). Note that the new chains have a pattern of red sandwiched between two layers of green, indicating that a left side and right side attachment point go on to form the top and bottom of a chain segment. (Y cross-section 00:100 volume density). Figure 11. Plots of shear stress obtained from simulations compared to Trendler and Bose experimental data at 600mT showing: (figure 11(a)) Shear stress for several different initial particle packing random seeds being sheared at rates ranging from g =  1100s −1 to g  =9 s −1 . (figure 11(b)) Mean shear stress from the multiple seeds in figure 11(a). Linear regressions for the shear rates above 200 s −1 are used to calculate the yield stress.

Yield stress
To compare the simulation results to the data from Trendler and Bose, the force on the upper raft particles was stored every 5K time steps, summed and divided by the area to calculate the stress (equation (8)).
The last 100 calculated shear stresses for each shear rate were averaged to obtain the shear stress at each shear rate. Figure 11 shows the results of four different initial seeds for the particle packing using the number density distribution, compared against the Trendler and Bose experimental results. Simulations with number particle distributions for 00:100 mixtures were computationally inexpensive and agreed with Trendler and Bose's experimental data very well. The four simulations were run using the standard shear cell size and the 00:100 number distribution mixture. The different seeds were fairly consistent, though Seed3 had somewhat larger deviation from the experimental data, and showed some morphological differences ( figure 11(a)). Still, the maximum standard deviation in the simulated shear stress for the four seeds was 7.16% of the mean shear stress, at a shear rate of 36.4 s −1 . The experimental data was within one standard deviation of the mean of the four seeds at all data points but the lowest shear rate, 9 s −1 , which was too high in the simulations ( figure 11(b)). In addition  to the regular size shear cells, two additional simulations were run with a 500μmX500μmX500μm shear cell, to check the sensitivity to the shear cell size in simulations ( figure 12). The double sized shear cell simulations produced similar results to the 250μm cells. The volume distribution, in addition to being considerably more computationally expensive, does not match the experimental results as well as the number density distributions ( figure 13). Yield stress was calculated using a linear fit to data points with shear rates g >  200 s −1 ( figure 11(b)). The yield stress and errors are shown in table 3. The mean yield stress for all 00:100 number density simulations was 43.9kPa compared to 44.9kPa in Trendler and Bose's data, and the standard deviation in yield stress for these data sets was 1.28kPa or 2.9% of the mean yield stress. The volume distribution produced a 49.6kPa yield stress.
Our simulation with the 33:67 mixture, using number density distributions, was able to capture the shear stress enhancement found experimentally in bi-disperse MRF mixtures ( figure 14). The enhancement to the shear stress was somewhat over-pronounced however, overshooting the experimental value by approximately 10%. There was also some reduction in the knee where the shear stress drops off at low shear rates. Simulations using volume distributions for the 33:67 mixture were too computationally expensive for us to produce a full shear profile in simulations and are only examined qualitatively in this work. In addition to the shear stress versus rate simulations, we also simulated various field strengths and measured the resultant shear stress at 100 s −1 . Two seeds were used, and they produced very consistent shear stresses with a similarly shaped trend to the experimental data, however the deviation from experiment is considerably larger than the yield stress data, growing as the magnetic field strength approaches 300mT, with a maximum error of 31%. At 32mT the error also spikes, though the magnitude of the shear stress is very low at that point, therefore, small deviations in the magnetic forces on the particle raft can create large errors in the simulated shear stress ( figure 15).
To examine the sensitivity of the simulated shear stresses to the contact force model we ran simulations without some of the contact forces enabled. These simulations were run using a 00:100 mixture with a number density particle distribution and were all run using the same initial particle packing seed. Eliminating the rolling resistance caused the largest reduction in shear stress of any of the forces and resulted in a less smooth roll-off at Figure 15. Trendler and simulation shear stress versus magnetic field for a 00:100 number density distribution. The shape of the trend is comparable, though the simulated shear stress droops in the middle compared to the experimental data from Trendler and Bose. The deviation from experimental is also large near zero, though the magnitude of the shear stress is very low at low fields. Figure 16. 00:100 Number density distribution simulations run with various model forces neglected. Neglecting rolling resistance or or using a monodisperse particle distribution resulted in reduced yield stress and some morphological differences in the shape of the shear rate/stress curve. Removing cohesion or using a stationary background fluid (though still including Stoke's drag) produced slightly increased shear stress but remained very similar morphologically to the full force simulations. low shear rates. Using a monodisperse mixture resulted in a high initial shear stress that dropped off quickly. Simulations run using a stationary background fluid with Stoke's drag had slightly increased shear stress, as did simulations run without cohesion ( figure 16). The closest match to the experimental data came from the simulations with all the forces considered, demonstrating the importance of maintaining a high fidelity contact model in MRF simulations.

Computational costs
Simulations for the 00:100 mixture with 250μm cube shear cells were run on the Deep Thought 2 High Performance computing cluster at University of Maryland consisting of dual Intel E5-2680v2 processors and 128 GB of 1866 MHz DDR3 memory per node. The number density distributions simulations were run with 18 cores on a single node. The simulations used approximately 15 MB of RAM per core and took 4.64ms per timestep with about 12k particles. The average number of neighbors used to calculate magnetic interactions, per particle, using a neighbor length of 8μm plus the radii of the particles being evaluated as neighbors, was approximately 100. The monodisperse 6.7μm particle simulations were also run on 18 cores, using approximately 29 MB of RAM per core, and took 3.90ms per timestep with approximately 30,000 particles. Despite the larger particle counts, the absence of smaller particles resulted in only 25 neighbors per particle with which to calculate magnetic interactions. Volume density distribution simulations are considerably more computationally intensive because of the larger number of small particles in the simulation. The volume density simulations ran on 60 cores on three nodes and used 17 MB of memory per core taking 31ms per timestep to run the simulations with ≈60,000 particles. On average, approximately 500 neighbors per particle were used to calculate magnetic interactions. Double sized shear cell simulations with number density distributions were run on a standalone workstation with dual Intel E5-2699V4 processors and 128 GB of 2400 MHz DDR4 RAM. The simulations used 27 cores and 17 MB of RAM per core, taking 30.1ms per timestep with ≈93,000 particles. An average of 100 neighbors per particle were used for magnetic interaction calculations.
Simulations for the 33:67 mixture were run on the standalone workstation. The number distribution used 40 cores and 33 MB per core for a runtime of 72ms per timestep with ≈200,000 particles. On average, approximately 220 neighbors per particle were used for magnetic interactions. For the volume distribution, the neighbor length was required to be shortened to 2μm plus the radii of the particles to prevent neighbor list overflows. On 40 cores, 80 MB per core were used, taking 1152ms per timestep to run the simulation with ≈910,000 particles and an average of approximately 300 neighbors per particle. Due to the computational resources required for the volume distribution, the simulation was only run for a single shear rate to observe the morphology of the structures and to analyze how that shear stress compared to the values for the same shear with a number density distribution.

Conclusion
Including realistic contact forces and particle size distributions has allowed us to accurately model magnetorheological fluids with higher fidelity than previous efforts. We have identified qualitative differences in the behavior of MRF composed of a distribution of particle sizes compared to the behavior of more uniformly sized mixtures, such as the resistance to the formation of sheets, lack of crystalline structures, and the reorganization of particle chains under shear. We have also captured the yield stress enhancement effect of MRFs composed of bidisperse mixtures of particles. The tradeoff to the high-fidelity model is the potential for much higher computation times in the case of particle volume size distributions which contain many, much smaller particles, increasing the total particle number and the number of close neighbors included in magnetic interaction computations. For the 33:67 mixture, considerably more computing resources would need to be devoted to get data for volume distributions. Despite the computational limitations, the increased fidelity improves our understanding of the microstructures and behaviors of magnetorheological fluids, some of which are contrary to results obtained from uniform particle sizes. Furthermore, the LIGGGHTS base code allows for the implementation of arbitrary MRF mixtures, allowing for rapid evaluation of new formulations. We've also matched experimental data for bulk properties of MR fluids, providing high confidence validation of the magnetic model and its integration with the LIGGGHTS framework. Additionally, the model has been used to simulate non-standard MRF operation modes. The behavior of an electropermanent magnet-based jamming, magnetorheological valve was replicated using LIGGGHTS simulations with the same magnetic model as was used in this work. [28].
Due to the use of the LIGGGHTS open source framework, it is possible to easily integrate new force models into simulations, creating a magnetic, granular simulation testbed, including a coupled fluid model using the supported CFDEM framework which integrates LIGGGHTS with OpenFoam. In keeping with the open source