Simulation of Magnetically-Actuated Functional Gradient Nanocomposites

Magnetically-actuated functional gradient nanocomposites can be locally modulated to generate unprecedented mechanical gradients that can be applied to various interfaces and surfaces through following the design principles of natural biological materials. However, a key question is how to modulate the concentration of magnetic particles using an external magnetic field. Here, we propose a model to obtain the gradient concentration distribution of magnetic particles and mechanical gradients. The results show that three states exist when the magnetic force changes in the z direction, including the unchanging state, the stable gradient state, and the over-accumulation state, which are consistent with experiment results. If both radial and axial magnetic forces are present, the inhomogeneity of magnetic–particle distribution in two dimensions was found to break the functional gradient. Furthermore, the size effects of a functional gradient sample were studied, which indicated that adjusting the magnetic force and diffusion constant would enable larger nanocomposites samples to generate functional gradients.


Introduction
Biological materials with functional gradient designs within interfacial and surface regions are finding increasing use in bioapplications [1][2][3][4][5][6][7][8][9].These can be inspired by many biomaterials in nature that are reported to be mechanically gradient, such as tree frog toe pads, squid beaks, sandworm jaws [10][11][12][13], etc.Meanwhile, interface/surface failures for synthetic materials are ubiquitous in various fields, including the coating industry [14,15], biomedical implants and restorations [16,17], flexible and wearable electronics [13,18], additive manufacturing, bioinspired adhesion [19,20], etc.Most of these failures come from mechanical/thermal mismatches in integrated systems containing dissimilar materials or structures.However, it remains a challenging scientific and technological task to transplant these concepts of functional gradients into the design of synthetic interfaces/surfaces for extreme performance.It was found that a variety of biomedical applications can be achieved in systems with magnetic particles [5].Recently, Wang et al. developed a new manufacturing technique for controllable functional gradient nanocomposites (FGNCs) via use of combined magnetic and mechanical actuations; the results show that a near one order of magnitude smooth change of the local elastic moduli (~2.2 GPa to ~9 GPa) can be enabled within regions as narrow as 10 microns [21].These results indicate that precise control over the spatial distribution of magnetic particles in a monomeric fluid is important, and further understanding is required.
The evolution of magnetic particles in an applied field can be affected by several factors, including magnetic field, viscous drag, gravity, Brownian dynamics, particle/fluid interactions, and particle-particle interactions (dipole-dipole interactions) [22][23][24].Two different approaches (Eulerian approach and Lagrangian approach) are used to model magnetic-particle transport in fluid, based on whether the Brownian motion has a significant impact on particle dynamics.In the Eulerian approach, particle evolution is characterized by particle concentration, and it is governed by a diffusion equation that accounts for both magnetic force-induced drift and Brownian diffusion [25][26][27].In the Lagrangian approach, particles are treated as discrete entities, as in molecular dynamics, and the Newtonian equation of each particle is solved independently at every time step [28].Generally, the effects of Brownian motion are negligible in the Lagrangian approach.Based on the Eulerian approach, the transport and accumulation of magnetic nanoparticles in a passive magnetophoretic system are studied [29], which is commonly used to discover the accumulation based on the assumption of ultra-low concentrations.
However, precise control of time-dependent spatially varying particle concentrations is needed in an experiment to acquire good, functional gradient materials, especially when the concentration is not low enough.To acquire tunable mechanical gradients in FGNCs, a permanent magnet was used to drive the magnetic-responsive nanofillers in a monomeric fluid until a steady-state concentration gradient was formed.In order to understand and design experiments, a more precise theoretical prediction is needed.Here, the drift-diffusion model, with a high concentration, was proposed to simulate functional gradient nanocomposites in polymers, in order to control the accumulation of magnetic particles in space at different times.Different magnetic forces will influence not only the accumulation of particles, as proposed in previous research, but also particle distribution in different times or spaces.At the same time, particle motion will be affected by a demagnetization field, which is induced by magnetic particle accumulation.In this paper, we completed a series of simulations in order to obtain the gradient concentration of nanoparticles at different scales and with different force distributions.

Theoretical Methods
Magnetic particle transport can be treated as a drift diffusion process, supposing that a magnetic particle is extremely small with a diameter less than 100 nm in simulations and experiments.According to Fick's second law of diffusion, magnetic particle volume concentration n can be described by: where J is the total flux vector of the particles; t is time; and S is the source term, which includes exchanges that occur at material interfaces or other chemical reactions.The total flux is the sum of chemical term J D and drift term J F .J D = −D∇n is the result of chemical concentration diffusion, and J F = vn is the result of the action of magnetic forces on the particles, where D = µk B T is the diffusion coefficient and v is particle drift velocity.Here, k B is the Boltzmann's constant; η is the fluid viscosity (Stokes' approximation); and µ is the mobility of a particle, which can be calculated by µ = 1/(6πηR p ), where R P is the effective hydrodynamic radius.Drift velocity v can be written as where v f is the fluid velocity.F m and F g are the external magnetic and fluidic forces on particles, respectively.It should be noted that the gravitational force is usually negligible compared with the magnetic force.Therefore, the conservation equation (Equation (1)) can be rewritten as [28,29]: When the particle concentration is high, it needs to be taken into account in the drift-diffusion process that particle motion actually takes place on spatial sites occupied by the atoms of materials.Hence, the governing equation should be described as: where f (n) is the factor that can be assumed to be (1 − n).
The magnetic force is calculated using the "effective" dipole moment method, in which the particle is modeled as an "equivalent" point dipole with an effective moment m p,eff .The magnetic force F m on the dipole is given by: where µ f is the permeability of the fluid and H a is the externally applied magnetic field at the center of the particle.Moment m p,eff can be determined using a magnetization model that takes into account self-demagnetization and the magnetic saturation of the particles.The dipole moment is m p,eff = V p M p , where V p is the volume of the particle, and: where M p is the magnetization of the particle; χ p = µ p /µ 0 − 1 is the susceptibility of the particle; and µ p is permeability.H in is the internal magnetic field, which can be calculated by applied magnetic field H a minus demagnetization field H demag .H demag is the self-demagnetization field that accounts for the magnetization of the particle.Here, H demag = M/3 is used for a uniformly magnetized spherical particle with magnetization M in free space.Thus: where f (H a ) is a function of magnetic susceptibilities, and can be expressed by [24]: where χ p and χ f are the magnetic susceptibilities of the particle and fluid, and M sp is the saturation magnetization of the particle.Thus, the magnetic force on the magnetic particle can be written as: For example, the magnetic field distribution (H az ) and magnetic force (F mz ) along the axial of a cylindrical rare-earth magnet can be calculated using [30]: where z is the distance from the sample to the magnet, and L m and R m are the length and radius of the magnet, respectively.If we need to estimate the radial force, the method in Appendix A can be used to calculate both radial and axial force.
Assuming that the scale of the sample is relatively small (~20 µm), the strength of the magnetic force can be seen as a uniform distribution due to the fact that distance Z 0 is constant.
Substituting Equation (11) to Equation ( 3), the drift-diffusion equation can be simplified as: Since the particle concentration is high, the drift-diffusion transport equation should be modified, compared with the dilute case.First, the mobility coefficient should be modified due to particle accumulation.According to the literature [31], in a polymer liquid, the mobility of nanoparticles can be affected by several factors, including particle size, solvent viscosity, and solute concentration.The mobility of particles decreases with solution volume fraction with the power relation D(n) ~n−α , where α is the parameter to be determined under different experimental conditions.
During the accumulation of magnetic particles in space, there will be a long-range interaction coming from the magnetic charge interaction.Here, we consider the magnetic charge interaction as a demagnetization field (H md ) in the whole matrix.Then, the magnetic force on the dipole is: The dipole moment of the magnetic particle m p,e f f = V p M p should be modified, and: The magnetic force on a magnetic particle can be written as: where H is the total magnetic field and H = H a + H md .
In order to calculate demagnetization field H md in the matrix, Gauss's law was introduced to solve the magnetic scalar potential and demagnetization field: where M(r,H) is the spatially, position-dependent, local magnetization; H a is the applied magnetic field; µ(r) = δ + χ m (r) is relative permeability; and δ is the Dirac delta function.
Here, the magnetic scalar potential ϕ m (r) is introduced to solve for demagnetization field H mdj (r) = −∇ j ϕ m (r).
In this paper, the diffusion equation is solved numerically using the finite volume method (FVM).Details of the FVM scheme have been described previously [29].Briefly, for one-dimensional analyses, the computational domain is discretized at an interior node, as shown in the following equation: where dx is the length of a computational cell; N z dx is the length of the sample; and N z is the number of computational nodes.F i+1/2 is a discretized version of the flux at the edges of the computational cell.Due to the convection term in the diffusion equation, the upwind numerical scheme is used, which can be described as:

Results and Discussion
In our previous experimental work [21], prepared core-shell nanoparticles (Fe 3 O 4 @SiO 2 , ~20 nm in diameter) could be well dispersed in polymer matrices at a volume fraction as high as 36 vol %.An initial particle concentration of 18 vol % and a specimen thickness of ~10 µm was used to produce function gradient nanocomposites under a magnetic field from a magnet.Two states of concentration distribution were obtained for a low and medium magnetic field (50 mT and 150 mT), respectively, corresponding to an almost-unchanged and smooth gradient filler distribution.However, when the field was further increased to, for example, 300 mT, the computed concentration gradient became notably steeper; almost all of the particles would accumulate on one side of the sample.
In this work, the FVM method was employed to complete the simulation of one-dimensional and two-dimensional models.The grid sizes are 100 dx and 100 dx × 100 dx, respectively.The parameters are shown in Table 1, which were obtained from the experiments [21].Here, the magnets were made from neodymium-iron-boron (NdFeB), and the nanoparticles (NPs) were Fe 3 O 4 with an initial concentration of 0.1.We should note that the fluid viscosity of the monomeric fluid used here was much higher than that of the liquid.The concentration-dependent diffusion coefficient can be assumed as , where D 0 is the diffusion coefficient in pure solvent, n ξ d is the crossover concentration, and α, β are coefficients to be determined from experiments.Here in our simulation, α = 0.53, n ξ d = 0.01, and β = 1.52 were assumed to take into account the concentration-dependent mobility [31].To solve Equation ( 16), we assumed that the magnetic charge could be evaluated as ρ m = −∇M = −∇nm p,e f f .We firstly completed a one-dimensional simulation that considered a 20 µm width sample; then, a two-dimensional model with 20 × 20 µm was used.The size effects were studied, with dx that varied from 2 × 10 −7 m to 2 × 10 −5 m, and, consequently, the total sample sizes were from 20 µm to 2 mm.As shown in Figure 1a, we can see that with the distance changing from 1 mm to 5 mm causes the magnetic field strength under the magnet to drop markedly, and the magnetic force on the nanoparticle (NP) decreases sharply, from 1.2 fN to 0.01 fN. Figure 1b,c shows the axial force and radial force distributions at different positions under the magnet.It can be seen that the radial force is very small compared with the axial force in the center region.The force will decrease in the axial direction, but will increase in the radial direction.Therefore, the magnetic force is determined by both the external magnetic field and the distance from the sample to the magnet.We completed three simulations of both one dimension and two dimensions using three different magnetic forces: 0.02 fN, 0.4 fN, and 1.2 fN.According to the magnetic force calculation of Equation ( 9), a magnetic force of 1.2 fN is the point at 1 mm under magnet.As shown in Figure 1d, the magnetic force is so large that the particle can almost be assembled in the right half of the sample.If the magnetic force is smaller (0.4 fN), due to the diffusion effect, the stable particle concentration gradually grows, from left to right.If the magnetic force is smaller than 0.02 fN, there is almost no change in the particle concentration in the sample.To compare the results with the experiment, we changed the particle concentration to the numbers of particles existing in the domain (Figure 2).The distribution of particles in Figure 2a-c comes from the assumption that the magnetic force only exists in the z direction; thus, the distribution of the cross-section of the sample (the plane perpendicular to the x direction) can be simplified to a one-dimensional model, which corresponds with the distribution in Figure 1d.It is easy to see the difference among the three states, and the results of the three different concentration distributions were consistent with experiments.
It is obvious that if the sample under the magnet is not at the center point, there will be a radial force, as seen in Figure 1c.Here, we consider the sample close to the magnet as having both the radial and axial force components.For simplicity, we assume that the magnetic force is constant in the y and z directions to study magnetic particle accumulation evolution.In Figure 3, the particle concentration shows high inhomogeneity when the force exists in both the z and y directions.If the force component in the y direction is zero, the accumulation will only exist in the z direction.If the force in the y direction is small (Figure 3b), accumulation in the y direction will exist.Due to the drift in the z and y directions, most magnetic particles will accumulate in the upper right corner.When the force in the y direction increases, the accumulation in the upper right corner will be larger, until it is sufficient to destroy the concentration gradient and mechanics gradient.This is the reason why in experiments, a larger cylindrical magnet is needed to make the magnetic force homogeneous in the sample in the y direction (radial direction) [21].However, it is interesting that by adjusting the magnetic force in the z and y directions, one can precisely control the highly inhomogeneous or particular magnetic particle concentration distributions, and thus the high mechanics inhomogeneity in a more complex design.As the sample scale increases, the diffusion process will not be the dominant process, and the magnetic force will change at different positions of the sample.Therefore, the force will gradually decrease when the distance from the magnet increases.Here, we consider the distance between the sample and the bottom of magnet as 2 mm.The scales of the samples are changed from 200 µm to 2 mm; other parameters used are the same as those listed in Table 1.
As shown in Figure 4, when the scale is small, a gradient distribution of concentration will be formed due to the diffusion process.However, when the size scale increases (e.g., 2 mm), the magnetic force will be the dominant driving force.A particle will be driven to the left side, where the magnetic force is larger than that on the other side of the sample.Another important problem is that particle accumulation is much slower than in a small-scale case.As seen from Figure 4a (2 mm case), the gradient distribution of the whole scale is difficult to obtain.
We derived the gradient concentration distribution in large-scale sample, and several points should be noted.First, the diffusion process should be accelerated, which is the reason why, in an experiment, the mechanical vibrations that can accelerate the diffusion process are needed.Second, the force should be increased to accelerate the accumulation process, and it can be described as: Increasing the radius or magnetic field gradient can increase the magnetic force.For verification, we increased the diffusion process speed 2000 times, and increased the radius of the magnetic particle to 100 nm.For the 2-mm size scale, the result can be seen in Figure 4b; the gradient distribution of nanoparticles will be obtained at 3 × 10 4 s.For comparison, we fixed the particle radius and accelerated the diffusion process using different times (in Figure 4c).The drift process dominates the evolution, and most particles will be driven to the right side of the sample.As the diffusion constant is gradually increased, the diffusion process competed with the drift process.A gradient concentration will be formed when the acceleration is sufficiently large.Similarly, the radius of the nanoparticles was increased, from 20 nm to 100 nm, and the diffusion constant was fixed.As shown in Figure 4d, the force was small due to the smaller particle size, and the diffusion process was the dominant process.As a nanoparticle increases in size, the driving force also increases, and the stable gradient state forms in some specific cases.Theoretical predictions indicated that, for the material compositions adopted in this study, the FGNCs should be able to be scaled up to ~200 µm using a realistic magnetic field and processing duration.For larger length scales, such as millimeters to even centimeters, the concept of FGNCs still works; however, the material parameters (e.g., particle size, viscosity of resin matrix, etc.) and the processing conditions (e.g., magnetic field intensity, processing duration, etc.) need to be adjusted in order to achieve an optimal concentration gradient.Additionally, the method of layer-by-layer deposition, as presented by Libanori et al. [32], may be combined with FGNCs to create large objects with continuous and significant mechanical gradients.

Conclusions
A drift-diffusion model with a high particle concentration was proposed to study the concentration evolution of the magnetic nanoparticles under magnetic force, and the three stable states obtained from the simulations were consistent with experiments.In addition, the size effect of the sample was also studied, and it was found that large-scale gradient distribution of magnetic particles can be obtained by precisely adjusting the diffusion process and the drift process.Increasing the magnetic force and diffusion constant will enable gradient distribution in large-scale samples.The two-dimensional model also demonstrated that the nanoparticle concentrations can be highly inhomogeneous.The results showed that drift-diffusion can be an effective model to predict and design functional gradient concentrations of magnetic-responsive nanoparticles.

Figure 1 .
Figure 1.(a) Magnetic field and force distribution in the z direction of magnet; (b,c) magnetic force in the axial and radial directions in different locations; (d) concentration distribution after 1.2 × 10 3 s with different applied forces.

Figure 2 .
Figure 2. (a) Cross section of the nanoparticles (NPs) distribution in the sample after 1.2 × 10 3 s evolution when the magnetic force is 0.02 fN; (b) the gradient concentration distribution under the force of 0.4 fN after 1.2 × 10 3 s evolution; (c) the over-accumulation state when the magnetic force is 1.2 fN after 1.2 × 10 3 s evolution.

Figure 3 .
Figure 3. Two-dimensional nanoparticle concentration distributions after 1.2 × 10 3 s when the force in the z direction is fixed with 0.4 fN and the force in the y direction is selected as (a) 0, (b) 0.008 fN, (c) 0.2 fN, and (d) 0.4 fN.

Figure 4 .
Figure 4. (a) Nanoparticle concentration along the normalized distance (e.g., distance from sample bottom divided by sample thickness) for sample thicknesses ranging between 0.1 and 2.0 mm at a sufficiently long processing duration of 1.2 × 10 4 s; (b) real-time evolution of the particle concentrations for a 2.0-mm (2000 µm) thick sample when the diffusion constant is 2000 times larger, and the magnetic particle radius is 100 nm; (c) nanoparticle concentration evolution after 1.2 × 10 4 s when particle size is fixed, and the diffusion constant is changed; (d) nanoparticle concentration evolution after 1.2 × 10 4 s when the diffusion constant is fixed, and the nanoparticle radius is changed.