Simulating the Structure of Magnetic Fluid Using Dissipative Particle Dynamics Method

1. Introduction

Magnetic fluid (MF) is a type of stable colloidal solution where magnetic nanoparticles are dispersed in a liquid carrier, aided by the presence of surfactants and the Brownian motion[1]. Since their invention, magnetic fluids have been widely applied in various research fields, such as biomedical engineering[2], environmental engineering[3,4], aerospace technology[5], Mechanical Engineering[6,7] and sensing technology[8-11]. Due to the physical properties of magnetic fluids themselves, particularly their relationship with optics and magnetic fluid structure, it is crucial to understand the properties of magnetic fluids at a microscopic level. This understanding can be achieved through modeling and analysis using simulation systems for magnetic fluids. Moreover, numerical simulation studies of the microstructure of magnetic fluids under the influence of an external magnetic field not only offer a more economical and environmentally friendly approach but also allow for the observation of changes that are difficult to capture in experiments, thus compensating for the limitations of experimental.

In recent years, the study of microstructure has attracted significant interest among scientific researchers. To enhance research efficiency, protect the environment, and reduce costs, numerous efforts have been devoted to investigating numerical methods for studying the structure of magnetic fluids. These methods include Molecular Dynamics (MD)[12], Monte Carlo (MC)[13], Lattice Boltzmann (LB)[14], Brownian Dynamics (BD)[15], and Dissipative Particle Dynamics (DPD)[16].In 2013, LV et al. laid the foundation for transmission simulation by using the Monte Carlo method to construct the microstructure of magnetic fluids (MF). Subsequently, they experimentally revealed the movement of magnetic particles in magnetic fluid films under different magnetic fields, demonstrating the feasibility of microstructure simulation and providing effective theoretical support for transmission simulation[17]. In 2016, Li et al. employed the DPD-based algorithm MVVA-MVVA to simulate the 3D microstructure of ferrofluids. They introduced the model of magnetic fluids and systematically studied the influences of magnetic particle interaction strength and volume fraction, yielding various aggregate structures that aligned well with those obtained using other numerical approaches, demonstrating the effectiveness of the DPD-based approach[18]. In 2021, Xu et al. utilized a research method combining mathematical analysis and performance testing to expand the application scope of MR fluids. A micromechanical model considering volume fraction and magnetic induction effects on microstructure evolution is proposed for MR fluids with MWCNTs/GO coated CI particles. The shear yield stress test conducted on self-prepared MR fluids using an MCR302 rheometer is compared with theoretical values from the model to validate its effectiveness[19].In the past 20 years, MD, MC and LB methods were the main approaches in the simulation of structure of magnetic fluids. However, there is a very limited literature dedicated to the study of structure of magnetic fluids using a DPD-based method, but DPD method is a very effective mesoscopic simulation technology for simulating the magnetic fluids structure. Moreover, in the aforementioned studies, there was no in-depth study on the effect of surfactants on the structure of magnetic fluids.

In this work, the Dissipative Particle Dynamics (DPD) method is employed to calculate the microstructure of magnetic fluids. In order to accurately reflect the real state, a two-dimensional simulation model was first established based on Fortran language, capturing various aggregate structures. The simulation results structures qualitatively agree well with actual experimental observations. Subsequently, the system studied the effects of the weight of carrier liquid molecules and the strength of magnetic particle interactions on the microstructure of magnetic fluid particles. This study provides valuable guidance for subsequent experimental re-search.

2. Methods and Frameworks

In the late 20th century, Hoogerbrugge and Koelman first proposed a mesoscopic simulation technique called Dissipative Particle Dynamics (DPD) for simulating complex fluid dynamics phenomena[20]. The popularity of the DPD model stems from its simple algorithm and wide applicability. Since the Dissipative Particle Dynamics (DPD) method is a mesoscale numerical simulation technique, the solvent molecular clusters are considered as dissipative particles. Therefore, in the case of magnetic fluids under consideration, there are two types of particles: magnetic particles and dissipative particles. The advantage of this method is that it gives access to longer time and length scales compared to what is achievable by conventional MD simulations. Figure 1 illustrates a schematic diagram of a typical Magnetic fluid.

2.1. Dissipative Particle Particle Models

We assume a system composed of N magnetic particles with mass m dispersed in a carrier liquid at thermodynamic equilibrium. In this study, the magnetic fluid is set as spherical particles, and we only need to consider the translational motion of the magnetic particles while neglecting their rotational motion.

There are three forces acting on a DPD particle: the conservative force ${\mathit{F}}_{ij}^C$ exerted by other particles, the dissipative force ${\mathit{F}}_{ij}^D$ which produces viscosity effects and drag forces among particles, and the random force ${\mathit{F}}_{ij}^R$ that generates thermal motion of particles. Taking these forces into account, the equation of motion for particle $i$ can be expressed as follows:

$$m_{\mathrm{d}}{\displaystyle \frac{\mathrm{d}{\mathbf{v}}_i}{\mathrm{d}t}}={\sum }_{j(\ne i)} {\mathbf{F}}_{ij}^{\mathrm{C}}+{\sum }_{j(\ne i)} {\mathbf{F}}_{ij}^{\mathrm{D}}+{\sum }_{j(\ne i)} {\mathbf{F}}_{ij}^{\mathrm{R}}$$

(1)

in which

$${\mathbf{F}}_{ij}^{\mathrm{C}}=\mathsf{\alpha }{\mathrm{w}}_{\mathrm{R}}\left(r_{ij}\right){\mathbf{e}}_{ij}$$

(2)

$${\mathbf{F}}_{ij}^{\mathrm{R}}=\sigma w_{\mathrm{R}}\left(r_{ij}\right){\mathbf{e}}_{ij}{\zeta }_{ij}$$

(3)

$${\mathbf{F}}_{ij}^{\mathrm{D}}=-\gamma w_{\mathrm{D}}\left(r_{ij}\right)\left({\mathbf{e}}_{ij}\cdot {\mathbf{v}}_{ij}\right){\mathbf{e}}_{ij}$$

(4)

In the given equations, $m_d$represents the mass of dissipative particle $i$, and $v_i$ represents its velocity. The subscripts indicate that, for example, ${\mathit{F}}_{ij}^C$ represents the force applied to particle $i$ by particle j. Moreover, the constants $\alpha$, $\sigma$, and $\gamma$ denote the strength of the repulsive conservative force, random force, and dissipative force, respectively. The weight functions $w_D\left(r_{ij}\right)$ and $w_R\left(r_{ij}\right)$ are defined in such a way that inter-particle forces decrease as the distance between two particles increases. The expression used to describe $w_R\left(r_{ij}\right)$ is as follows:

$$w_{\mathrm{R}}\left(r_{ij}\right)=\left\{\begin{array}{cc}1-{\displaystyle \frac{r_{ij}}{d_{\mathrm{c}}}}& \mathrm{for}{r}_{ij}\le d_{\mathrm{c}}\\ 0& \mathrm{for}{r}_{ij}>d_c\end{array}\right.$$

(5)

The weight functions $w_D\left(r_{ij}\right)$ and $w_R\left(r_{ij}\right)$, as well as the constants $\gamma$ and $\sigma$, must satisfy the following relationships, respectively:

$$w_{\mathrm{D}}\left(r_{ij}\right)=w_{\mathrm{R}}^2\left(r_{ij}\right)$$

(6)

$${\sigma }^2=2\gamma kT$$

(7)

In the equations provided above, $d_c$ represents the apparent diameter of dissipative particles. ${\mathit{r}}_{ij}$ denotes the relative position($r_{ij}=\left|{\mathit{r}}_{ij}\right|$), ${\mathit{r}}_{ij}={\mathit{r}}_i-{\mathit{r}}_j$; ${\mathit{e}}_{ij}$ is the unit vector indicating the direction from particle $i$ to particle $j$, expressed as ${\mathit{e}}_{ij}={\mathit{r}}_{ij}/r_{ij}$; ${\mathit{v}}_{ij}$ represents the relative, given by ${\mathit{v}}_{ij}={\mathit{v}}_i-{\mathit{v}}_j$; $k$ represents Boltzmann's constant, and $T$ denotes the liquid temperature. Additionally, ${\zeta }_{ij}$ is a random variable that induces the random motion of the particles.

By integrating Eq. (1) with respect to time over a small time interval $\Delta t$ from $t$ to $t+\Delta t$, we can derive the finite difference equations that govern the particle motion in simulations.

$$\mathsf{\Delta }{\mathbf{r}}_i={\mathbf{v}}_i\mathsf{\Delta }t$$

(8)

$$\begin{array}{rr}\mathsf{\Delta }{\mathbf{v}}_i=& {\displaystyle \frac{\alpha }{m_{\mathrm{d}}}}{\sum }_{j(\ne i)} w_{\mathrm{R}}\left(r_{ij}\right){\mathbf{e}}_{ij}\mathsf{\Delta }t-{\displaystyle \frac{\gamma }{m_{\mathrm{d}}}}{\sum }_{j(\ne i)} w_{\mathrm{R}}^2\left(r_{ij}\right)\left({\mathbf{e}}_{ij}\cdot {\mathbf{v}}_{ij}\right){\mathbf{e}}_{ij}\mathsf{\Delta }t\\ & +{\displaystyle \frac{(2\gamma kT{)}^{1/2}}{m_{\mathrm{d}}}}{\sum }_{j(\ne i)} w_{\mathrm{R}}\left(r_{ij}\right){\mathbf{e}}_{ij}{\theta }_{ij}\sqrt{\mathsf{\Delta }t}\end{array}$$

(9)

where ${\theta }_{ij}$ is the stochastic variable that must adhere to the following stochastic properties:

$$⟨{\theta }_{ij}⟩=0$$

(10)

$$⟨{\theta }_{ij}{\theta }_{i^{\prime }{j}^{\prime }}⟩=\left({\delta }_{i{i}^{\prime }}{\delta }_{i{j}^{\prime }}+{\delta }_{i{j}^{\prime }}{\delta }_{j{i}^{\prime }}\right)$$

(11)

in which ${\delta }_{ij}$ is the Kronecker delta. During the simulation, the stochastic variable ${\theta }_{ij}$ is sampled from a uniform or normal distribution with zero average value and unit variance.

2.2. Particle Models

As shown in Figure 2, a magnetic particle is represented as a spherical particle with a central point dipole(positive and negative charges coincide at the same point) and is coated with a uniform steric layer, also known as a surfactant layer. In terms of notation, the diameter of the particle is denoted as $d_s$, the thickness of the steric layer as $\delta$, and the overall diameter including the steric layer as $D(=d_s+2\delta )$.

The magnetic interaction energy between particles $i$ and $j$, $u_{ij}^{\left(m\right)}$, the particle-field interaction energy, $u_i^{\left(H\right)}$, and the interaction energy resulting from the overlap of the steric layers, $u_{ij}^{\left(V\right)}$, are expressed as follows:

$$u_{ij}^{\left(\mathrm{m}\right)}={\displaystyle \frac{{\mu }_0}{4\pi r_{ij}^3}}\left\{{\mathbf{m}}_i\cdot {\mathbf{m}}_j-3\left({\mathbf{m}}_i\cdot {\mathbf{t}}_{ij}\right)\left({\mathbf{m}}_j\cdot {\mathbf{t}}_{ij}\right)\right\}$$

(12)

$$u_i^{\left(\mathrm{H}\right)}=-{\mu }_0{\mathbf{m}}_i\cdot \mathbf{H}$$

(13)

$$u_{ij}^{\left(\mathrm{V}\right)}=kT{\lambda }_V\left\{2-{\displaystyle \frac{2r_{ij}/d_{\mathrm{s}}}{t_{\delta }}}\mathrm{l}\mathrm{n}\left({\displaystyle \frac{D}{r_{ij}}}\right)-2{\displaystyle \frac{r_{ij}/d_{\mathrm{s}}-1}{t_{\delta }}}\right\}$$

(14)

In this context, ${\mu }_0$ represents the permeability of free space, ${\mathit{m}}_i$ denotes the magnetic moment with $m_0=\left|{\mathit{m}}_i\right|$, and ${\mathit{t}}_{ij}$ is a unit vector defined as ${\mathit{r}}_{ij}/r_{ij}$, where ${\mathit{r}}_{ij}$=${\mathit{r}}_i-{\mathit{r}}_j,r_{ij}=\left|{\mathit{r}}_{ij}\right|$. $\mathit{H}$ denotes the applied magnetic field($H=\left|\mathit{H}\right|$), and $t_{\delta }$ is defined as the ratio of the steric layer thickness $\delta$ to the radius of the solid portion of the particle, which can be expressed as $2\delta /d_s$. The nondimensional parameter ${\lambda }_V$, as observed in Equation (14), quantifies the strength of the steric particle-particle interaction relative to the thermal energy, represented by ${\lambda }_V=\pi d_s^2{n}_s/2$, where $n_s$ represents the number of surfactant molecules per unit area on the surface of the particle. The forces acting on particle i can be derived from equations (12) and (14).

$${\mathbf{F}}_{ij}^{\left(\mathbf{m}\right)}=-{\displaystyle \frac{3{\mu }_0}{4\pi r_{ij}^4}}\left[-\left({\mathbf{m}}_i\cdot {\mathbf{m}}_j\right){\mathbf{t}}_{ij}+5\left({\mathbf{m}}_i\cdot {\mathbf{t}}_{ij}\right)\left({\mathbf{m}}_j\cdot {\mathbf{t}}_{ij}\right){\mathbf{t}}_{ij}-\left\{\left({\mathbf{m}}_j\cdot {\mathbf{t}}_{ij}\right){\mathbf{m}}_i+\left({\mathbf{m}}_i\cdot {\mathbf{t}}_{ij}\right){\mathbf{m}}_j\right\}\right]$$

(15)

$${\mathbf{F}}_{ij}^{\left(V\right)}={\displaystyle \frac{kT{\lambda }_V}{\delta }}\times {\displaystyle \frac{{\mathbf{r}}_{ij}}{r_{ij}}}\mathrm{l}\mathrm{n}\left({\displaystyle \frac{D}{r_{ij}}}\right),\left(d_{\mathrm{s}}\le r_{ij}\le D\right)$$

(16)

The movement of magnetic particles is described by Newton's equations and is discretized in time to obtain finite difference equations that govern the particle motion in simulations.

$$\mathsf{\Delta }{\mathbf{r}}_i={\mathbf{v}}_i\mathsf{\Delta }t$$

(17)

$$\mathsf{\Delta }{\mathbf{v}}_i={\sum }_{j(\ne i)} {\mathbf{F}}_{ij}\mathsf{\Delta }t/m_{\mathrm{m}}$$

(18)

In this equation, $m_m$ denotes the mass of the magnetic particles, and it is important to note that

$${\mathbf{F}}_{ij}={\mathbf{F}}_{ij}^{\left(\mathrm{m}\right)}+{\mathbf{F}}_{ij}^{\left(V\right)}$$

(29)

In our approach, each colloidal particle is represented as a collective of dissipative particles. Normally, the interaction between a magnetic particle and the surrounding dissipative particles is considered to be the interaction between the surrounding dissipative particles and the constituent dissipative particles of the magnetic particle. However, in real dispersions, the nature of interactions between colloidal particles and solvent molecules depends on the specific characteristics of the dispersion. These interactions are heavily influenced by the mass-to-diameter ratio of colloidal particles relative to solvent molecules, as well as the properties of the interaction potential. As a result, instead of considering a colloidal particle solely as a group of dissipative particles, it may be feasible to employ a model potential to describe the interaction between magnetic and ambient dissipative particles. Figure 3 illustrates a schematic representation of the collision between a magnetic particle and a dissipative particle.

A suitable model for particle interaction that can be utilized is the Lennard-Jones model. The Lennard-Jones model is characterized by the following equation, which describes the interaction energy $u_{ip}$ between dissipative particle $p$ and magnetic particle $i$:

$$u_{ip}=4\epsilon \left\{{\left({\displaystyle \frac{d_{\mathrm{c}}}{{r_{ip}}^{\prime }}}\right)}^{\mathrm{m}}-{\left({\displaystyle \frac{d_{\mathrm{c}}}{{r_{ip}}^{\prime }}}\right)}^n\right\}$$

(20)

In this equation, $\epsilon$ represents a constant that quantifies the strength of the interaction, ${\mathit{r}}_{ip}^{\text{'}}={\mathit{r}}_i^{\text{'}}-{\mathit{r}}_p$, $r_{ip}^{\text{'}}=\left|{{\mathit{r}}_{ip}}^{\text{'}}\right|$. The position vectors of the magnetic particle $i$ and dissipative particle $p$ are denoted as ${\mathit{r}}_i$ and ${\mathit{r}}_p$, respectively. In addition, ${\mathit{r}}_i^{\text{'}}$,the position vector representing the inscribed sphere is written as:

$${\mathbf{r}}_i^{\prime }={\mathbf{r}}_i-\left(D-d_c/2\right){\stackrel{\mathit{ˆ}}{\mathit{r}}}_{ip}$$

(21)

in which ${\stackrel{\mathit{ˆ}}{\mathit{r}}}_{ip}={\mathit{r}}_{ip}/r_{ip}$, ${\mathit{r}}_{ip}={\mathit{r}}_i-{\mathit{r}}_p$ and $r_{ip}=\left|{\mathit{r}}_{ip}\right|$. In the present simulation, the well-known Lennard-Jones potential is utilized, which is derived from the model potential by setting $m$ to 12 and $n$ to 6 in Eq. (20).The force ${\mathit{F}}_{ip}^{\left(\mathrm{int}\right)}$ exerted on dissipative particle $\mathrm{p}$ by magnetic particle $\mathrm{i}$ is determined by the interaction energy described in Eq. (20), as follows:

$${\mathbf{F}}_{ip}^{\left(\mathrm{int}\right)}=4n\epsilon \left\{{\displaystyle \frac{m}{n}}{\left({\displaystyle \frac{d_{\mathrm{c}}}{{r_{ip}}^{\prime }}}\right)}^m-{\left({\displaystyle \frac{d_{\mathrm{c}}}{{r_{ip}}^{\prime }}}\right)}^n\right\}{\displaystyle \frac{{\stackrel{\mathbf{ˆ}}{\mathbf{r}}}_{ip}}{{r_{ip}}^{\prime }}}$$

(22)

The subsequent quantities are employed for the non-dimensionalization of equations: $d$ for distance, $m_m$ for mass, $kT$ for energy, ${\left(kT/m_m\right)}^{1/2}$ for velocity, $d{\left(m_m/kT\right)}^{1/2}$ for time and $kT/d$ for force. By using these quantities, we have

$$\mathsf{\Delta }{\mathbf{r}}_i^{\mathrm{*}}={\mathbf{v}}_i^{\mathrm{*}}\mathsf{\Delta }{t}^{\mathrm{*}}$$

(23)

$$\begin{array}{rr}\mathsf{\Delta }{\mathbf{v}}_i^{\mathrm{*}}=& {\displaystyle \frac{1}{m_d^{*}{d}_c^{*}}}{\alpha }^{\mathrm{*}}{\sum }_{j(\ne i)} w_{\mathrm{R}}\left(r_{ij}^{\mathrm{*}}\right){\mathbf{e}}_{ij}\mathsf{\Delta }{t}^{\mathrm{*}}-{\displaystyle \frac{1}{{\left(m_d^{*}\right)}^{1/2}{d}_{\mathrm{c}}^{\mathrm{*}}}}{\gamma }^{\mathrm{*}}{\sum }_{j(\ne i)} w_{\mathrm{R}}^2\left(r_{ij}^{\mathrm{*}}\right)\left({\mathbf{e}}_{ij}\cdot {\mathbf{v}}_{ij}^{\mathrm{*}}\right){\mathbf{e}}_{ij}\mathsf{\Delta }{t}^{\mathrm{*}}\\ & -{\displaystyle \frac{1}{{\left(m_d^{*}\right)}^{3/4}{d}_c^{*1/2}}}{\left(2{\gamma }^{\mathrm{*}}\right)}^{1/2}{\sum }_{j(\ne i)} w_{\mathrm{R}}\left(r_{ij}^{\mathrm{*}}\right){\mathbf{e}}_{ij}{\theta }_{ij}\sqrt{\mathsf{\Delta }{t}^{\mathrm{*}}}-{\displaystyle \frac{1}{m_d^{*}}}{\sum }_k {\mathbf{F}}_{ki}^{\left(\mathrm{i}\mathrm{n}\mathrm{t}\right)\mathrm{*}}\mathsf{\Delta }{t}^{\mathrm{*}}\end{array}$$

(24)

where

$$w_{\mathrm{R}}\left(r_{ij}^{\mathrm{*}}\right)=\left\{\begin{array}{ll}1-r_{ij}^{\mathrm{*}}/d_c^{*}& \mathrm{for}{r}_{ij}^{\mathrm{*}}/d_c^{*}\le 1\\ 0& \mathrm{for}{r}_{ij}^{\mathrm{*}}/d_c^{*}>1\end{array}\right.$$

(25)

$${\alpha }^{\mathrm{*}}=\alpha {\displaystyle \frac{d_{\mathrm{c}}}{kT}},{\gamma }^{\mathrm{*}}=\gamma {\displaystyle \frac{d_{\mathrm{c}}}{{\left(m_{d}kT\right)}^{1/2}}}$$

(26)

The number density of dissipative particles is nondimensionalized as

$$n_d^{*}=n_d{d}^2=n_d{d}_c^2{\left(d/d_c\right)}^2={\hat{n}}_d^{\mathrm{*}}/d_c^{*2}$$

(27)

In addition to $n_d^{*}$ , the nondimensional density ${\hat{n}}_d^{*}$ based on the diameter of dissipative particles may be useful for quantifying the packing characteristics of the dissipative particles. The nondimensional number density of magnetic particles is expressed as $n_m^{*}=n_m{d}^2$.

$$\mathsf{\Delta }{\mathbf{r}}_{\mathbf{i}}^{\mathrm{*}}={\mathbf{v}}_i^{\mathrm{*}}\mathsf{\Delta }{t}^{\mathrm{*}}$$

(28)

$$\mathsf{\Delta }{\mathbf{v}}_i^{\mathrm{*}}={\sum }_{j\ne i} {\mathbf{F}}_{ij}^{\mathrm{*}}\mathsf{\Delta }{t}^{\mathrm{*}}+{\sum }_p {\mathbf{F}}_{ip}^{\left(\mathrm{i}\mathrm{n}\mathrm{t}\right)\mathrm{*}}\mathsf{\Delta }{t}^{\mathrm{*}}$$

(29)

$$\begin{array}{ll}{\mathbf{F}}_{ij}^{\left(\mathrm{m}\right)\mathrm{*}}=& -3\lambda {\displaystyle \frac{1}{r_{ij}^{4\mathrm{*}}}}\left[-\left({\mathbf{n}}_i\cdot {\mathbf{n}}_j\right)t_{ij}+5\left({\mathbf{n}}_i\cdot {\mathbf{t}}_{ij}\right)\left({\mathbf{n}}_j\cdot {\mathbf{t}}_{ij}\right){\mathbf{t}}_{ij}-\left\{\left({\mathbf{n}}_j\cdot {\mathbf{t}}_{ij}\right){\mathbf{n}}_i\right.\right.\\ & \left.\left.+\left({\mathbf{n}}_i\cdot {\mathbf{t}}_{ij}\right){\mathbf{n}}_j\right\}\right]\end{array}$$

(30)

$${\mathbf{F}}_{ij}^{\left(\mathrm{V}\right)\mathrm{*}}={\lambda }_V{\displaystyle \frac{1}{t_{\delta }^{\mathrm{*}}}}\times {\mathbf{t}}_{ij}\mathrm{l}\mathrm{n}\left({\displaystyle \frac{1}{r_{ij}^{\mathrm{*}}}}\right)$$

(31)

in which ${\mathit{F}}_{ij}^{*}={{\mathit{F}}_{ij}}^{\left(m\right)*}+{{\mathit{F}}_{ij}}^{\left(V\right)*}$, ${\mathit{n}}_i$ is the unit vector denoting the direction of the magnetic moment ${\mathit{m}}_i,$expressed as ${\mathit{n}}_i={\mathit{m}}_i/m_0\left(m_0=\left|{\mathit{m}}_i\right|\right).$ The nondimensional parameter $\lambda$ in Eq. (30) is the strength of magnetic particle interactions relative to the thermal energy, expressed as $\lambda ={\mu }_0{m}_0^2/4\pi D^{3}kT$.A parameter ${\lambda }_s={\left(D/d_s\right)}^3\lambda \left(={\mu }_0{m_0}^2/4\pi {d_s}^{3}kT\right)$, defined based on the diameter of the solid component. The expression of the force between a dissipative and a magnetic particle is written in nondimensional form as

$${\mathbf{F}}_{ip}^{\mathrm{*}}={\lambda }_{\epsilon }\left\{{\displaystyle \frac{m}{n}}{\left({\displaystyle \frac{d_{\mathrm{c}}^{\mathrm{*}}}{r_{ip}^{\prime \mathrm{*}}}}\right)}^m-{\left({\displaystyle \frac{d_{\mathrm{c}}^{\mathrm{*}}}{r_{ip}^{\prime \mathrm{*}}}}\right)}^n\right\}{\displaystyle \frac{{\stackrel{\mathbf{ˆ}}{\mathbf{r}}}_{ip}}{r_{ip}^{\prime \mathrm{*}}/d_{\mathrm{c}}^{\mathrm{*}}}}$$

(32)

in which ${\lambda }_{\epsilon }$ is a nondimensional parameter representing the strength of the interaction, expressed as ${\lambda }_{\epsilon }=4n\epsilon /\left(kT{d}_c^{*}\right)$.

3. Modeling and Simulation of Magnetic Fluids

When the number of magnetic fluid particles is set to N, the area of this two-dimensional simulation box is represented by ${LL}^2$, $(LL{)}^2={\displaystyle \frac{N\pi }{4\phi }}$when magnetic particle volume fraction $\phi$ is known. Unless explicitly stated otherwise, the designated duration for this task is 1000000 steps. The parameters utilized for the simulations conducted in this study are as follows: ${\lambda }_s=4$, ${\lambda }_{\epsilon }=10$,${\gamma }^{*}=10$,${\alpha }^{*}={\gamma }^{*}/10$,$d_c^{*}=0.4$,$m_d^{*}=0.01$,${\hat{n}}_d^{*}=1$,$\Delta t^{*}=0.0001$, ${\lambda }_V=120$, volume fraction of magnetic particle $\phi$ is 0.2, the number of magnetic fluid particles is 100 and surfactant layer thickness ${\delta }^{*}=0.15$.The total number of simulation steps, ${\mathrm{N}}_{\mathrm{timemx}}$,when the condition of $\Delta t^{*}{N}_{\mathrm{timemx}}=100$ is met, it is anticipated to be adequate.

The time integration method is one of the key factors affecting the accuracy of simulation results in magnetic fluid simulations based on dissipative particle dynamics. Different time integration methods can have an impact on simulation results, so it is necessary to run and compare these methods under the same system parameters. The following are common time integration methods:

The Euler Scheme (ES) is one of the simplest integration methods, but it is still an effective time integration method in some cases. This method is based on obtaining the position and velocity of the next time step from the position and velocity of the previous time step. The expression is as follows:

$${\mathit{r}}_i\left(t+\Delta t\right)={\mathit{r}}_i\left(t\right)+\Delta t{\mathit{v}}_i\left(t\right)$$

(33)

$${\mathit{v}}_i\left(t+\Delta t\right)=\mathit{v}\left(t\right)+\Delta t{\displaystyle \frac{{\mathit{F}}_i\left(t\right)}{m_i}}$$

(34)

$${\mathit{F}}_i\left(t+\Delta t\right)={\mathit{F}}_i\left({\mathit{r}}_i\left(t+\Delta t\right),{\mathit{v}}_i\left(t+\Delta t\right)\right)$$

(35)

In the 1960s, the Verlet algorithm was proposed[21], in which an explicit expression for particle velocity appeared as follows:

$${\mathit{r}}_i\left(t+\Delta t\right)={\mathit{r}}_i\left(t\right)+\Delta t\mathit{v}\left(t\right)+{\displaystyle \frac{{\left(\Delta t\right)}^2{\mathit{F}}_i\left(t\right)}{2m_i}}$$

(36)

$${\mathit{F}}_i\left(t+\Delta t\right)={\mathit{F}}_i\left({\mathit{r}}_i\left(t+\Delta t\right)\right)$$

(37)

$${\mathit{v}}_i\left(t+\Delta t\right)={\mathit{v}}_i\left(t\right)+{\displaystyle \frac{\Delta t\left[{\mathit{F}}_i\left(t\right)+{\mathit{F}}_i\left(t+\Delta t\right)\right]}{2m_i}}$$

(38)

Figure 4 shows the simulation results obtained using different algorithms. The small blue-green circles in the figure represent dissipative particles, while the large orange circles represent magnetic particles. It can be observed from the figure that, under the same conditions, there are a large number of magnetic particles gathering together to form coarse chain structures in the image obtained using the ES algorithm. However, at this time, the interaction between magnetic particles does not play an absolute dominant role relative to thermal energy, and under such conditions, large magnetic particle aggregation structures should not be formed.

The radial distribution function, g(r), defines the probability of two particles appearing at a given distance, r, and serves as an important tool for quantitatively characterizing the microstructure of Magnetic fluid. Specifically, for the two-dimensional radial distribution function, it is defined as follows[22]:

$$g\left(q\right)={\displaystyle \frac{1}{N/A}}{\displaystyle \frac{s}{2\pi q△q}}$$

(33)

In this context, $A$ represents the area of the computational region containing $N$ particles, and $s$ represents the number of particles within a circular shell distance of $q$ from the target particle, in the range of $q$ to $q$ + $△q$.

The radial distribution function corresponding to the ES algorithm has a higher peak near the integer radius r, indicating that compared to the Verlet algorithm, there are more magnetic particles distributed at these distances, corresponding to the chain like structure generated by the magnetic chain in Figure 5(a). The radial distribution function of the Verlet algorithm is relatively smooth, corresponding to the fact that most magnetic particles are still in an independent state.

Figure 6 shows the dynamic evolution of magnetic particle structure over time. The small circles in the graph represent dissipative particles, while the large circles represent magnetic particles. From the figure, it can be seen that at the beginning of the simulation, the magnetic particles were uniformly distributed, and then dissipative particles were filled in the simulation area. With the increase of dimensionless time t*, the initially disordered distribution transforms into movement along the direction of the magnetic field, leading to aggregation and the emergence of long-chain structures. The aggregation structures of these ferromagnetic particles are in excellent agreement with the qualitatively obtained results from experiments in[17].

4. Computational Findings and Analysis

4.1. Effect of the Molecular Mass of Carrier Liquid on the Structure of Magnetic Fluid

When the mass of dissipative particles is small, magnetic particles tend to form chain-like aggregation structures with larger radial sizes, as illustrated in Figure 7a. With an increase in the mass of dissipative particles (Figure 7b), the aggregation morphology of magnetic particles transitions into thinner chain-like structures. Further increasing the mass of dissipative particles leads to the formation of discrete cluster-like distributions, as shown in Figure 7c,d. These structural changes, analyzed through differential mass variations, indicate that the aggregate behavior of magnetic particles is highly sensitive to the mass of dissipative particles. As shown in Fig. 3.8, with the increase of the mass of dissipative particles, the peak of the radial distribution function near r=1 shows a decreasing trend, indicating the gradual dissociation of cluster structures among magnetic particles. This phenomenon suggests a corresponding increase in the proportion of independent magnetic particles in the system. Moreover, the results demonstrate that changes in the weight of the carrier liquid directly influence the motion behavior of magnetic particles in a magnetic field, thereby regulating the morphology of chain-like structures in magnetic fluids.

4.2. Effect of the Magnetic Particle Interaction Strength on the Structure of Magnetic Fluids

In a strong magnetic field environment, the interaction strength between magnetic particles plays a crucial role in regulating their structure. With the increase of interaction strength, the interaction between magnetic particles intensifies, leading to a stronger tendency for magnetic particles to form arranged chain-like structures. This evolutionary process of structure is essential for a profound understanding of the be-havior and properties of magnetic particles in a strong magnetic field.

As shown in Fig. 9(a), when the interaction strength between magnetic particles is low (λs=4), although some particles come into contact with each other to form cluster structures, the structure is relatively loose, indicating a significant influence of Brownian motion. At the same time, there are a large number of particles in a free state without interacting with other magnetic particles. With an increase in λs, as shown in Fig. 9(b) and (c), the cluster structure gradually transitions to a chain-like structure, while the number of free particles begins to decrease. As shown in Fig. 9(d), when the interaction strength between magnetic particles increases to a higher value, its influence on the structure of magnetic fluids relative to thermal energy takes precedence. The majority of magnetic particles form chain-like structures, and are ordered along the direction of the magnetic field. This structural transition process is clearly reflected in the radial distribution function: as shown in Fig. 3.10, with the increase of λs, the peak intensity at r=1 continues to strengthen, reflecting the structural evolution of particles from an independent state to chain-like aggregation.

Figure 9. Effect of the strength of interaction between magnetic particles on the structure of magnetic fluids.

Figure 9. Effect of the strength of interaction between magnetic particles on the structure of magnetic fluids.

Figure 10. Radial distribution function corresponding to the strength of magnetic particle interactions.

Figure 10. Radial distribution function corresponding to the strength of magnetic particle interactions.

5. Conclusions

This study employs dissipative particle dynamics (DPD) method to simulate the structure of magnetic fluids. Initially, a computational model of magnetic fluids consisting of magnetic nanoparticles and dissipative particles is established. A comparison of different time integration methods is conducted, and the numerical results obtained from the simulation are qualitatively consistent with the simulated and experimental results reported in the literature. Additionally, the impact of solvent molecular weight and the strength of magnetic particle interactions on the structure of magnetic fluids is analyzed based on the radial distribution function. The DPD-based method proves to be highly effective in simulating the microscopic structure of magnetic fluids. As for future research, this method will be extended to study the structural characteristics of three-dimensional magnetic fluids.