*2.1. Magneto-Hydrodynamics*

Magneto-Hydro-Dynamics (MHD) equations consist of equations of electromagnetics and hydrodynamics. Electromagnetics equation relates the current density *J* to the magnetic field strength *H* [48–50].

$$
\vec{\text{curl}} \, \vec{H} = \vec{\text{J}} \tag{1}
$$

$$\begin{array}{c} \stackrel{\rightarrow}{div} \stackrel{\rightarrow}{J} = 0 \end{array} \tag{2}$$

In the case that the substance such as fluid has a velocity, *J* is calculated as follows:

$$
\overrightarrow{J} = \sigma(\overrightarrow{E} + \overrightarrow{v} \times \mu\overrightarrow{H})\tag{3}
$$

where σ is the electrical conductivity, *v* is the velocity of flow, μ is the magnetic permeability, and *E* is the electric field intensity. The body force or external force → *Fe* which will be used in DPD equations, Section 2.2 (as hydrodynamic equations) is calculated as follows:

$$
\overrightarrow{F}\_{\mathfrak{c}} = \overrightarrow{\mathfrak{J}} \times \mu \overrightarrow{H} \tag{4}
$$

#### 2.1.1. Analytical Solution for MHD in Simple Channel

By using the source term *Fex* calculated in Equation (4) in the steady equation of fluid motion, applying the assumption of uniform conductivity of liquid between two parallel walls (simple channel), implementing *H*0 perpendicular to the walls in *z*-direction while implementing *E* in y-direction, and ignoring the gravity effect, a fluid motion should occur in *x*-direction satisfying the following momentum equation [47,51].

$$0 = \frac{-\partial p}{\partial \mathbf{x}} + \mu f\_y H\_0 + \eta \frac{\partial^2 v}{\partial z^2} \tag{5}$$

where *p* is pressure term, η is the dynamic viscosity, and *Jy* = σ(*<sup>E</sup>* − *<sup>v</sup>*μ*H*0).

Substituting the *Jy* in Equation (5) and assuming that velocity is zero at the walls and shear stress is zero at the middle of channel, would result in the following solution of MHD flow in the simple channel.

$$v = V\_v (1 - \frac{\cos \mathbf{h}(\frac{Ha \ z}{L})}{\cos \mathbf{h}(Ha)}) \tag{6}$$

where, *Vv* = ( *E*μ*H*0 + −∂*p* ∂*x* σμ2*H*0<sup>2</sup> ) and *Ha* = μ*H*0*<sup>L</sup>*( σμ )0.5.

If the parameter *H*0 varies while the other parameters are kept constant, the value of *Ha* (or Hartmann number) will change linearly. Hartmann number or *Ha* is the ratio of magnetic force to viscous force.

$$\text{If we assume that } V\_v = 1 \text{ then } E = (1 - \frac{\frac{-\partial p}{\partial x}}{\sigma \mu^2 H\_0^2})\mu H\_0 \text{ and } J\_y = \sigma ((1 - \frac{\frac{-\partial p}{\partial x}}{\sigma \mu^2 H\_0^2})\mu H\_0 - \mu v H\_0).$$

#### *2.2. Dissipative Particle Dynamics Method*

The DPD method treats the simulated system as a cloud of beads each having the mass, *mi*, position vector → *ri*, and velocity vector →*v i*. The evolution of beads' velocity vectors follows the basic kinematic and dynamic laws of motion [47].

$$
\overrightarrow{v}\_i = \frac{d\overrightarrow{r}\_i}{dt} \tag{7}
$$

$$m\_i \frac{d\overrightarrow{\boldsymbol{\upsilon}}\_i}{dt} = \overrightarrow{\boldsymbol{F}}\_i \tag{8}$$

The net force exerted on bead *i* and *Fi*, can be decomposed in two parts, namely the force exerted by an external force-field → *Fe* like gravitational, electrical, or magnetic fields, and the intermolecular forces exerted by polymer → *<sup>F</sup>*p,*<sup>i</sup>* and fluid particles *<sup>j</sup>i* → *f ij*:

→

$$
\overrightarrow{F}\_i = \overrightarrow{F}\_{\text{P},i} + \sum\_{j \neq i} \overrightarrow{f}\_{i\hat{j}} + \overrightarrow{F}\_e \tag{9}
$$

In DPD method, it is assumed that the intermolecular force exerted by bead *j* on bead *i* and → *f ij*, consists of three parts, namely the conservative force, the dissipative force, and the random force [52].

$$
\stackrel{\rightarrow}{f}\_{ij} = \stackrel{\rightarrow}{f}\_{ij}^{\mathbb{C}} + \stackrel{\rightarrow}{f}\_{ij}^{\mathbb{D}} + \stackrel{\rightarrow}{f}\_{ij}^{\mathbb{R}} \tag{10}
$$

Of course, → *f ij* would be equal to , as the third law of motion indicates. → *Fext* is the external force such a electro-osmotic force or magnetic force [32].

The conservative force can be described using the following force filed:

$$\stackrel{\rightarrow}{f}\_{ij}^{\mathbb{C}} = \begin{cases} a\_{ij}(1 - r\_{ij}/r\_{\mathbb{C}})\mathfrak{f}\_{ij} & r\_{ij} < r\_{\mathbb{C}} \\ 0 & r\_{ij} \ge r\_{\mathbb{C}} \end{cases} \tag{11}$$

where *aij* is the maximum repulsion force between beads *i* and *j*, *rij* is the distance between those beads, i.e., *rij* = → *rij* , and *r*ˆ*ij* = → *rij*/ → *rij* is the unit vector pointing from bead *j* to bead *i*, and *r*c is the cut-o ff radius beyond which, the intermolecular forces are assumed to diminish e ffectively.

The dissipative intermolecular force is calculated as follows:

$$\stackrel{\rightarrow}{f}\_{ij}^{\mathcal{D}} = -\gamma \omega^{\mathcal{D}} \ r\_{ij} (\mathfrak{f}\_{ij}, \overrightarrow{\boldsymbol{v}}\_{ij}) \mathfrak{f}\_{ij} \tag{12}$$

where γ is a constant that determines the strength of dissipative force. The weighting function ω<sup>D</sup> is calculated as follows:

$$
\omega^{\rm D}(r\_{ij}) = \begin{cases}
\left(1 - r\_{ij}/r\_c\right)^2 & r < r\_{ij} \\
0 & r \ge r\_{ij}
\end{cases}
\tag{13}
$$

On the other hand, the random force is calculated as follows:

$$\stackrel{\rightarrow}{f}\_{\,ij}^{\mathbb{R}} = -\sigma^{\mathbb{R}} \omega^{\mathbb{R}} \, r\_{i\bar{j}} \theta\_{i\bar{j}} \mathfrak{f}\_{i\bar{j}} \tag{14}$$

in which, the constant σ*<sup>R</sup>* determines the strength of random force. → *v ij* is the relative velocity vector between beads *i* and *j*. θ*ij* is a random number chosen from a symmetric Gaussian distribution having the zero mean and unit variance. σ*<sup>R</sup>* and ω<sup>R</sup> relate to γ and ω<sup>D</sup> as follows:

$$
\omega^{\mathbb{R}} = \sqrt{\omega^{\mathbb{D}}} \tag{15}
$$

$$
\sigma^R = \sqrt{2\gamma k\_\mathcal{B} T} \tag{16}
$$

where *k*B is the Boltzmann constant and *T* is the temperature.

In this study, the modified velocity-Verlet algorithm is used for time integration of position and velocity of each particle is calculated explicitly. The Ref. [53] is suggested for more information. Implementation of wall boundary condition has been follow based on the ref [48,54–56].

## *2.3. Polymer Chain*

Polymer chain model consists of a number of masses and springs which are connected together. The mechanism of polymer chain motion influenced from interaction of fluid particles and polymer beads is depicted in Figure 1. According to number of di fferent types of particles, we have di fferent interactions. In this case, we have three types of interactions, including fluid to fluid particles, polymer and fluid, and polymer and polymer beads interactions. Each mass has its own repulsion force to fluid particles or other masses. Thus, the conservative force *Fc* should be modified for this type of simulation. In this paper, the harmonic spring force, → *f* S *i*,*i*±1, or → *<sup>F</sup>*p,*<sup>i</sup>* is used between beads which is added to *F*c in DPD formulation. For the beads consisting a polymer chain, the following conservative harmonicforcepresentstheintermolecularbonds[32,34,53].

$$\stackrel{\rightarrow}{\hat{f}}\_{i,i\pm1}^{\mathcal{S}} = -\mathcal{K}(r\_{i,i\pm1} - r\_{i,i\pm1}^{\mathcal{eq}}) \\ \hat{r}\_{i,i\pm1} = \stackrel{\rightarrow}{\hat{F}}\_{\mathcal{P}^{\hat{J}}} \tag{17}$$

where *req i*,*i*±1 is the equilibrium (zero-force) distance between two beads in the polymer chain and *K* is the harmonic bond constant.

**Figure 1.** Mechanism of polymer and fluid particles interaction and motion of polymer chain.

Due to difficulty of studying the properties of polymer chain in motion, some physical properties such as mean square radius of gyration and velocity of mass center are employed for finding effect of MHD on polymer chain motion.

In equilibrium, the radius of gyration is defined as [30,51]:

$$
\langle R\_{\rm G}^2 \rangle = \frac{1}{2N^2} \langle r\_{ij}^2 \rangle \tag{18}
$$

where, *rij* = →*ri* − →*r j* and *N* is the number of beads of polymer chain. The end to end distance (Ee) is defined as *R*2*G* = *r*2*N*1. Also, the center of mass velocity of system is calculated as follows:

$$w\_{\rm cm} = \frac{\sum\_{1}^{N} m\_{i} v\_{i}}{\sum\_{1}^{N} m\_{i}} \tag{19}$$

where *mi* is the mass of bead *i*. Also, the average kinetic temperature is calculated as follows [51].

$$
\langle k\_B T \rangle = \frac{m}{3n - 3} \langle \sum\_{1}^{n} v\_i^2 \rangle \tag{20}
$$

where *n* is the number of particles.
