Magneto-Permeability Effect in Ferrofluid Flow through Porous Media studied via Multiparticle Collision Dynamics

As more and more promising applications of magnetic nanoparticles in complicated environments are explored, their flow properties in porous media are of increasing interest. We here propose a hybrid approach based on the Multiparticle Collision Dynamics Method extended to porous media via friction forces and coupled with Brownian Dynamics simulations of the rotational motion of magnetic nanoparticles' magnetic moment. We simulate flow in planar channels homogeneously filled with a porous medium and verify our implementation by reproducing the analytical velocity profile of the Darcy-Brinkman model in the non-magnetic case. In the presence of an externally applied magnetic field, the non-equilibrium magnetization and friction forces lead to field-dependent velocity profiles that result in effective, field-dependent permeabilities. We provide a theoretical expression for this magneto-permeability effect in analogy with the magneto-viscous effect. Finally, we study the flow through planar channels, where only the walls are covered with a porous medium. We find a smooth crossover from the Poiseuille profile in the center of the channel to the Brinkman-Darcy flow in the porous layers. We propose a simple estimate of the thickness of the porous layer based on the flow rate and maximum flow velocity.


Introduction
Dynamics and flow of magnetic nanoparticles (MNPs) suspended in nonmagnetic viscous carrier fluids (ferrofluids) have been studied intensively over the last decades [1][2][3][4].These studies have almost exclusively focused on spatially homogeneous solvents.In many applications, however, one is interested in the flow through composite materials.Prominent examples are the transport through sand or other granular matter [5].On a coarse-grained level, transport phenomena through such materials can be described via a porous medium approach [5,6].Similarly, the porous medium approach is also frequently used in biomedical applications to describe transport through biological tissues and cancerous cells (see e.g.[7][8][9][10] and references therein).
Although MNPs find more and more promising biomedical and technical applications [1,2,11], to date, only a handful of studies address the flow properties of MNPs through porous media.One of the experimental studies on ferrofluid flow though sands and sediments observed a strong dependence on external magnetic fields [12].Similarly, experiments on the internal convection of ferrofluids flowing through a capillary tube filled with porous media showed that external fields could significantly enhance the thermal conductivity [13].Also the efficiency of ferrofluids for oil displacements in a sand-filled pipe was investigated experimentally and compared to finite-element simulations [14].
Fluid dynamics simulation of ferrofluids have also been performed to investigate their use as displacing fluid in fractured porous media [15].Other finite-element or finite-volume simulations have addressed some particular flow [16,17] and heat transfer [18] properties of ferrofluid flow through porous media.A porous medium approach was also used in finite-volume simulations of magnetic drug targeting of MNPs, coupling channel flow to adjacent tumor region via the permeable endothelium layer [19].These simulation studies relied on highly simplified constitutive models, typically neglecting internal rotations and corresponding non-equilibrium magnetization components.
In addition to classical fluid dynamics simulations such as finite-volume and finite-element methods, the Lattice Boltzmann scheme has been successfully used to describe flow through porous media [20].In these simulations, explicit scatterers for fluid motion are placed at fixed locations within the simulation cell.Using this method, ferrofluid permeation into a randomly structured porous medium has been simulated and shown to be sensitive to an applied magnetic field [21].As an alternative simulation approach, an extension of the highly versatile multi-particle collision dynamics method (MPC) [22] to transport in porous media has been proposed in Ref. [23] for non-magnetic fluids.In the latter, the effect of porous media on fluid transport is simply modelled as local damping, leading to very efficient simulation methods.Note that the Lattice Boltzmann methods put forward in Refs.[20,21] resolve the detailed fluid dynamics in the vicinity of individual grain boundaries.On the other hand, the MPC modelling proposed in Ref. [23] is suitable for a more coarsegrained level of description where the porous medium can be considered to be locally homogeneous.
Here, we use similar ideas to simulate ferrofluid flow through porous media via a MPC method that is coupled to Brownian Dynamics simulations of the MNP dynamics to model their internal rotations.In the absence of porous media, this approach was proposed and validated in Ref. [24], showing the correct incorporation of a reliable constitutive model for dilute ferrofluids.We here show how this model can be extended via friction forces to model ferrofluid flow through porous media.In the absence of an external magnetic field, we reproduce the Darcy-Brinkman velocity profile and clarify the interpretation of and relationship between the model parameters.Simulations of driven flow through planar channels show that flow properties can be manipulated by external magnetic fields.In particular, we observe an effective permeability that increases with increasing field strength before reaching a limiting value.Using kinetic theories of ferrofluids, we provide a theoretical expression of this magneto-permeability effect in close analogy with the magneto-viscous effect.Finally, we study driven flow through planar channels where only the channel walls are covered with a layer of porous material, with no porous medium present in the center region of the channel.

Continuum Level
Darcy's law predicts the flow velocity v through a porous medium when a pressure gradient ∇p is applied as [5] where η is the dynamic viscosity of the fluid.The empirical proportionality coefficient K is known as permeability of the porous medium.For isotropic porous media K = KI, with I the identity matrix, so that v = −(K/η)∇p.
In a finite domain, the Darcy-Brinkmann model provides a better description than Darcy's law [7].This model can be formulated as the stationary Navier-Stokes equation supplemented with an additional damping term proportional to an empirical parameter α, The density and kinematic viscosity of the fluid are denoted by ρ and ν = η/ρ, respectively.The phenomenological parameter α governs the strength of the damping term and is related to the permeability coefficient by α = ν/K.We here consider Reynolds numbers that are small enough so that the Forchheimer correction [7] is irrelevant.
For spatially homogeneous porosity, i.e. a position-independent α, the exact solution of Eq. ( 2) for one-dimensional channel flow v(r) = v(y)e x reads [20] where L denotes the width of the channel and no-slip boundary conditions on the channel walls have been assumed.The parameters appearing in the flow profile Eq. ( 3) are given by c = −(αρ) −1 dp dx and From Eq. ( 4) and the above relation α = ν/K, we find that the permeability can also be expressed as K = 1/r 2 .These relations will be useful for later analysis.
Note that for weak damping, we recover the usual Poiseuille profile from Eq (3), v(y) = −(2νρ) −1 dp dx y(L − y) + O(r 2 ).Conversely, increasing r leads to stronger and stronger deviations from the parabolic velocity profile.
The amount of fluid transported per unit time through a cross-section of the channel (known as volumetric flow rate in the three-dimensional case) is obtained from V = L 0 v x (y)dy.For the velocity profile (3) we obtain V = cL[1 − tanh(L * )/L * ] with reduced channel width L * = rL/2.As expected, the flow rate is proportional to the applied pressure gradient.For L * ≫ 1 we find V ≈ c[L − 2/r] corresponding to a plug flow, whereas for L * ≪ 1 we recover the equivalent of the Hagen-Poiseuille law for two-dimensional channels, V ≈ −(dp/dx)L 3 /[12ρν](1 + O((L * ) 2 )).

Mesoscopic Level: Particle-based model
Over the past decades, several particle-based methods for simulating fluid flow have been explored (see e.g.[25] and references therein).Contrary to more traditional fluid dynamics simulations, these mesoscopic methods are very flexible, straightforward to implement, and naturally include thermal fluctuations.In the present study, we employ one of these methods known as multi-particle collision dynamics (MPC) [22].One of the advantages of the MPC method is that viscoelastic fluids can be modeled rather straightforwardly [26].In particular, we have already proposed and successfully tested an MPC implementation of ferrofluid flow using a reliable constitutive model [24].To make the paper self-contained, we provide a short description of the standard MPC method, before specifying the extension to porous media.
Within the MPC method, the fluid is represented by a collection of N identical particles, each with mass m.If r i and v i denote the position and velocity of particle i, i = 1, . . ., N , particle dynamics is split into a streaming and a collision step.In the streaming step, particles are advanced for a time ∆t as where f i (t) is the total force acting on particle i at time t.While Eqs. ( 5), (6) are formally identical to those used in Molecular Dynamics simulations, the crucial idea of MPC as a mesoscopic method is to include in f i only external and body forces and to disregard inter-particle interactions in the streaming step.Instead, inter-particle interactions are accounted for via momentum exchange in the collision step.In this collision step, all particles i residing at time t in the same collision cell C i are updated simultaneously as In Eq (7), V Ci (t) denotes the center-of-mass velocity of the collision cell C i and R = R(α) a matrix, describing rotations around a randomly chosen axis by an angle ±α.Equation ( 7) models the effect of collisions among particles as rotations of their relative velocities.A local thermostat is present in Eq. ( 7) and described by the factor β th = T /T Ci , where T Ci is the instantaneous kinetic temperature of the collision cell C i and T a prescribed bath temperature.We use a two-dimensional square grid, so the collision cells are squares of side length a.Using a spatially fixed grid of collision cells violates Galilean invariance.Therefore, we follow common practice [27] and in each step shift the grid of collision cells by a vector with independent random components in [−a/2, a/2].Besides the time step ∆t which determines the mean-free path λ = ∆t k B T /m, the mean number of particles per collision cell Q is the other crucial parameter in the MPC model [22,28].Due to the importance of angular momentum conservation, we here follow our earlier work [24] and implement the angular momentum-conserving algorithm (denoted as MPC-DR in Ref. [28]), where the rotation angle α is not a free parameter but chosen as where The simplified collision rules between particles (7) do not resolve individual collisions, but ensure local conservation laws are obeyed.Therefore, the MPC method is numerically very efficient and hence can simulate hydrodynamic behavior on time and length scales larger than ∆t and a [22,25].
In the form described so far, the MPC method has been successfully applied to various flow problems for viscous and viscoelastic fluids [26], including ferrofluids [24].In order to describe flow through porous media, however, we need to introduce the effect the porous medium exerts on the fluid via inelastic collisions.Here, we follow the ideas put forward in Ref. [20] that on a mesoscopic level, the interaction of the fluid with the porous medium can be described as a local damping of the fluid velocity.Within the MPC method, this approach can be implemented straightforwardly as an additional friction force acting on the particles [23], where ξ(r) is a (possibly position-dependent) friction coefficient.For ξ = 0 we recover the original MPC model, whereas ξ > 0 describes velocity damping due to porous media.For the case of pressure-drive flow that we consider in the following, the force on particle i can be written as , where the external force due to an applied pressure gradient is f ext = −ρ −1 ∇p.
Eqs ( 5) -( 9) describe the MPC model of non-magnetic fluid flow through porous media.This model has essentially been proposed in Ref [23], where instead of adding the friction force ( 9), particle velocities are rescaled by a factor (1 − ξ∆t/m).

Hybrid MPC-BD model for FF flow
For magnetic fluids, the stationary momentum balance equation ( 2) must be supplemented by additional Kelvin-Helmholtz forces [4], where H and M denote the magnetic field and the magnetization, respectively.We assume the fluid to be non-conducting, therefore we must also satisfy the magnetostatic Maxwell equations where B = µ 0 (H+M) denotes the magnetic induction and µ 0 the permeability of free space.For a thorough description of ferrofluid hydrodynamics see e.g.Ref. [4].
To evaluate the force density (10), one needs to employ a model to calculate the field-and flow-dependent magnetization M. Unfortunately, even after 50 years of research, there is no commonly agreed magnetization equation available in the literature (see e.g.[29,30] and references therein).Therefore, we here consider dilute conditions where there is less controversy and adopt the classical kinetic model of Martsenyuk et al. [31] that has been studied frequently since [3,32,33].
To make the paper self-contained, we briefly present the MPC implementation of this model proposed in Ref. [24].Further details can be found in the original reference.In this model, one considers the rotational dynamics of an individual MNP within the rigid-dipole approximation under the influence of external magnetic fields and flow.From the balance of magnetic, flow and random torques, one finds that the orientation of the magnetic moment of particle i evolves to first order in the time step ∆t B by with where Ω Ci and h Ci are one half the local vorticity of the flow and the magnetic field, respectively, both evaluated at the center of collision cell C i .The Brownian relaxation time of a MNP is denoted by τ B and ∆W i are independent, three-dimensional Wiener increments over a time interval ∆t B .In the simulations presented here, we use a weak second order stochastic Heun scheme, where Eqs. ( 12) and ( 13) serve as predictor step.
With the magnitude of the magnetic moment µ of a single MNP, the instantaneous local magnetization in collision cell C i can be calculated as where N Ci (t) is the number of MPC particles in collision cell C i at time t and n denotes the number density of MNPs.We calculate the magnetization M in all collision cells according to Eq. ( 14) and use kernel-smoothing methods to find a discretization of the magnetization field M(r; t).From the discretized magnetization field, we calculate spatial gradients via finite-difference approximations and are thus able to evaluate the Kelvin-Helmoltz force density (10) within each collision cell.Further details on the method are given in Ref. [24].Note that we do not explicitly include the effect of the porous medium on the rotational motion of the MNPs.While it is plausible that inelastic collisions of fluid and MNPs with obstacles lead to an effective damping of the translational motion, their effect on rotations is less obvious, especially since we later consider rotational motion of MNPs on time scales long compare to fluid motion (by choosing τ B ≫ ∆t).One possibility would be to model this effect as an additional rotational friction.In this case, all results presented below still hold when adjusting τ B correspondingly for given ξ.Due to the uncertainties associated with such modelling, we chose to consider τ B constant in this initial study.
3 Non-magnetic fluid flow through porous media All numerical results and parameter values presented below are reported in dimensionless form, with the particle mass m and the linear size of the collision cell a as basic units, together with the reference thermal energy k B T ref .
Consequently, the units for time are t ref = a m/k B T ref and the units for the diffusion coefficient and kinematic viscosity are First, we consider non-magnetic fluids, i.e. the MPC particles experience no other external forces except friction forces (9) and external pressure gradients,

Self-diffusion in unbounded domain, no external forcing
To study self-diffusion under equilibrium conditions, no pressure gradient is applied, f ext = 0. Furthermore, to study self-diffusion in an unbounded domain, we consider in this section a periodic system without any walls present.Under these conditions and for the two-dimensional angular momentumconserving collision model adopted here, the diffusion coefficient was calculated using molecular chaos assumption [28] as with shows the diffusion coefficient as a function of the average number of MPC particles per collision cell Q for two selected temperatures.
We perform MPC simulations for a two-dimensional periodic system of size 30 × 30 with and different values for Q.From the particle mean-square displacement, ⟨[r i (t) − r i (0)] 2 ⟩ = 4Dt, we extract the diffusion coefficient D from a least-square fit and verified that the x-and y-components of the displacement agree with each other within numerical accuracy.Figure 1(a) shows that Eq. ( 15) provides a good representation of the numerical results, even though some quantitative discrepancies can clearly be discerned for small up to moderate values of Q.For this model, similar deviations from the theoretical predictions have been reported in Ref. [28] and attributed to the limitations of the molecular chaos assumption.
Equation (15) has been derived for a standard MPC fluid.We are not aware of a corresponding result for an MPC fluid in a porous material.Heuristically, we can describe the effect of local damping in the MPC scheme for porous materials by an effective time step δt in the free-streaming step [23], δt = (1 − ξ∆t/m)∆t.We assume that the reasoning presented in Ref. [28] leading to Eq. ( 15) remains otherwise unaltered.In particular, we assume that collisions occur sufficiently fast so that they are not affected by local friction effects.Thus, we suggest an approximation to the MPC diffusion coefficient for spatially homogeneous porous media as with the same expression for s K as given after Eq. ( 15). Figure 1(b) shows the variation of the effective diffusion coefficient with increasing value of the friction coefficient ξ.The agreement of the numerical results with Eq. ( 16) is satisfactory.In particular, we find that the effective diffusion coefficient decreases approximately linear with increasing ξ as predicted by Eq. ( 16).

Channel flow
In the following we consider the pressure-driven channel flow of a fluid through a porous medium.Within the MPC model, an applied pressure gradient dp/dx in flow direction is realized by the external force f ext = f ext e x acting on every particle, where f ext = −ρ −1 dp/dx.
To determine suitable values for f ext , we calculate the flow rate V by numerical integration over the flow profile and verify that this quantity varies linearly with f ext in the parameter regime studied here.We consider a planar channel of widths L = 32 and 64 and length 50.From Fig. 2 we find V/f ext approaches a limiting value with decreasing f ext .Within statistical uncertainties, we find the same limiting value for f ext ≲ 10 −3 .Therefore, we will use in the following f ext = 10 −3 unless stated otherwise.
Having chosen a suitable value for f ext , we perform a series of MPC simulations for different values of the friction coefficient ξ and analyze the resulting velocity profiles.For all conditions investigated, we find that the numerical velocity profiles are well described by the analytical profile for the Darcy-Brinkman model, Eq. (3).Fitting the numerical profiles to Eq. ( 3), we extract two fit parameters, c and r, from which the model parameters can be determined as follows.First, as shown in Sec.2.1, c is related to the damping parameter α in the Darcy-Brinkman equation ( 2) by α = f ext /c.Next, the permeability K is directly linked to the parameter r by K = 1/r 2 .Finally, the kinematic viscosity is given by ν = α/r 2 .Figure 3 shows the extracted damping parameter α and permeability K over a wide range of values for the friction coefficient ξ in the MPC model.Within numerical accuracy, we find that the Darcy-Brinkman damping parameter α is equal to the friction coefficient ξ used in the MPC model.Therefore, the newly introduced friction coefficient in the MPC model can be identified with the more familiar damping parameter in the Darcy-Brinkman approach.For a derivation of this result at least for inviscid fluids see Appendix A. From Fig. 3(b) we find that the permeability K decreases with increasing ξ.Except at very small values of ξ, we find that the decrease can be described as K ∼ 1/ξ to a very good approximation.It is interesting to note that the relations α = ξ and K = k 0 ξ −1 also hold for a corresponding Lattice Boltzmann implementation of flow through porous media [20], where the density of scatterers plays the role of ξ.
For small values of ξ, large uncertainties in model parameters extracted from fits to Eq. (3) are found.In this regime, the velocity profiles are close to parabolic, leading to ambiguities in the two-parameter fit to Eq. ( 3).
As a consistency test, we plot in Fig. 4(a) the parameter r governing the velocity profile (3) parametrically versus the damping parameter α that we determined for different values of the friction coefficient.In agreement with Eq. ( 4) we find from our simulations r ∼ α 1/2 .Note that this relation was also confirmed in Lattice Boltzmann simulations [20], while an earlier MPC implementation [23] recovered this relation only over a rather limited interval of ξ.
From the relations α = ξ and K = k 0 ξ −1 extracted from Fig. 3, we conclude that the kinematic viscosity ν = αK (see Sec. 2.1) is given by ν = k 0 , independent of the friction coefficient.In the absence of porous media, we have already established the value of the kinematic viscosity ν = 0.114 ± 0.001 for ∆t = 1.0 in Ref. [24].Since the current model reduces for ξ = 0 to a pure MPC fluid, we have performed some simulations for ξ = 0 and fitted the resulting velocity to a parabolic profile expected for Poiseuille flow.Thereby, we confirm the value of ν obtained earlier for the current choice of parameters and ∆t = 1.0 and find in addition ν = 0.320 ± 0.001 for ∆t = 0.2.In Fig. 4(b) we show the kinematic viscosity ν obtained from fits to the velocity profile (3) via the relation ν = α/r 2 in Sec.2.1.From Fig. 4(b) we indeed find that ν is independent of the friction coefficient within statistical uncertainties, consistent with our conclusions above.Therefore, we can identify ν with the kinematic fluid viscosity in the absence of a porous medium.

Ferrofluid fluid flow through porous media
In this section, we consider ferrofluid flow for the same conditions as studied in Sec.3.2, i.e. the steady-state flow through planar channels filled with porous media.As before, we assume spatially homogeneous porous materials so that the friction coefficient ξ in the MPC method, Eq. ( 9), is constant, independent of the position throughout the channel.Here, we employ the hybrid MPC-BD model of fluctuating ferrohydrodynamics and include the Kelvin-Helmholtz force into the momentum balance as described in Sec.2.3.
Assuming that the rotational relaxation of MNPs is slow compared to fluid motion, we choose τ B = 100 and set ∆t B = ∆t.Figure 5 shows velocity profiles obtained from MPC-BD simulations for selected parameter values of ∆t = 0.2, n * = 0.005, using ξ = 0.01 and ξ = 0.02.As found earlier, with increasing friction ξ, the flow velocity is reduced and deviates more strongly from the parabolic Poiseuille profile found in the absence of porous media.In addition, we observe that the external magnetic field also influences the flow profile.In particular, we find that increasing the magnetic field leads to a reduction of the overall velocity, corresponding to an increase in the effective viscosity of the fluid.This so-called magnetoviscous effect is well-known for magnetic nanoparticles suspended in viscous solvents [3].
Here, we analyze this phenomenon quantitatively for the flow through porous media.To this end, we fit the velocity profiles that we obtain numerically from the MPC-BD simulations to the profile (3) resulting from the Darcy-Brinkman model.For all parameter values investigated, we find that Eq. ( 3) provides an accurate description of our numerical results.Therefore, we can proceed with our analysis of the fitted coefficients c and r.First, we extract the damping parameter α from the amplitude c of the profile (3) via the relation α = f ext /c found in Sec.2.1, where f ext = −ρ −1 dp/dx.From Fig. 6(a) we find that the effective damping parameter α for ferrofluid flow is independent of the magnetic field strength and still given by the MPC friction coefficient ξ as in the non-magnetic/field-free case.Thus, within statistical uncertainty, the friction coefficient ξ in the MPC model is equal to the damping parameter α in the Darcy-Brinkman model (2) also in the magnetic case.In Fig. 6(b), we show the effective permeability K eff obtained from the parameter r in the fitted velocity profile (3) versus the magnetic field strength h for different values of the friction coefficient ξ.We observe that K eff increases with increasing field strength, reaching a limiting value for large h.In order to explain this finding, we use Eq. ( 4) to define the effective viscosity ν eff and evaluate this quantity with the values for the damping parameter α and the fit parameter r that we have determined.In Fig. 7(a), we show the viscosity change ∆ν eff (h) = ν eff (h) − ν(0), scaled with the viscosity ν(0) at zero field.We verified that ν(0) agrees with the viscosity ν obtained in Sec.3.2 for non-magnetic fluids with otherwise identical model parameters.From Fig. 7(a), we find that the relative viscosity change does not depend on the value of the friction coefficient ξ.Furthermore, the increase of ∆ν eff (h)/ν(0) is well described by the classical result for ferrofluids [31], where L(h) = coth(h)−h −1 denotes the Langevin function and ϕ the magnetic volume fraction.Therefore, the magnetoviscous effect can be seen in porous media in very much the same manner as in viscous solvent.
It should be mentioned that the model of Martsenyuk et al. [31] employed here can be justified only for dilute ferrofluids.Therefore, in future applications, smaller values for the number density n should be chosen with corresponding smaller viscosity changes.Here, we have chosen larger values of n to better illustrate the field-dependent effects captured by the simulation method.
Having rationalized the effective viscosity, we replot in Fig. 7(b) the effective permeability K eff from Fig. 6(b) not versus the magnetic field h but versus ν eff (h).Multiplying K eff with the friction coefficient ξ, all data fall nicely on the diagonal, showing that the effective permeability of ferrofluid flow through porous media is given by K eff (ξ, h) = ν eff (h)/ξ, which is well approximated by with ν the kinematic viscosity in the absence of magnetic fields.In the absence of a magnetic field, we recover K = K eff (α, h = 0), where we made use of the equality ξ = α.For increasing magnetic field strength h, we find that the effective permeability K eff increases monotonically to approach an asymptotic value for h ≫ 1.In analogy to the well-known magnetoviscous effect, one might speak of a corresponding "magneto-permeability effect" in ferrofluid flow through porous media.From Fig. 6(b) we observe that the magnetopermeability effect for the present model is well-described by Eq. ( 18).We close this section by considering the most important dimensionless groups determining the nature of fluid flow.First, the Reynolds number is defined as Re = U L/ν, where U is a characteristic flow velocity and L the channel width.For the present choice of parameters U ≈ 0.1 for L = 64 (see Fig. 5), so that we find typical Reynolds numbers Re ≈ 20, well in the laminar regime.Next, the Schmidt number Sc = ν/D is a measure for the importance of viscous versus molecular diffusion.For Q = 100 ≫ 1, we find D ≈ T ∆t/2 ≈ 0.01 so that Sc ≈ 30, indicating that we are indeed operating in the relevant regime for fluids, where collisions dominate over kinetic transport.Finally, the Mach number is defined as Ma = U/c s , where c s = 5k B T /3m denotes the speed of sound.In our case, we find typically Ma ≈ 0.25.It is well-known that particlebased methods like MPC do not strictly obey the incompressibility condition.The smaller the Mach number, the better incompressibility is restored.Finally, for viscoelastic fluids like ferrofluids, the Weissenberg number Wi = τ B U/L gives the ratio of viscous to elastic forces.Here, we typically find Wi ≈ 0.2 and therefore the simulations are performed in the Newtonian regime.
5 Flow through channels with walls covered by porous media In this section, we consider again driven flow through a parallel channel.This time, however, the porous medium is only present within a layer of width ℓ p on both walls, whereas the fluid is unperturbed in center.Within the MPC scheme described in Sec.2.2, this situation can be implemented conveniently by a position-dependent friction coefficient ξ(r) in Eq. ( 9).The porous layers are described by setting ξ(r) = ξ for y ≤ ℓ p or y ≥ L−ℓ p and ξ = 0 in the center, ℓ p < y < L − ℓ p .Fig. 8 Flow profile in a parallel channel where walls are covered with a porous layer, indicated by the grey shaded regions.Black and blue symbols correspond to h = 0 and h = 5, respectively.In panels (a) and (b), the width of the porous layer is ℓp = 0.2L and 0.3L, respectively.The model parameters are chosen as T * = 0.1, ∆t = 0.2, Q = 100, ξ = 0.02, f ext = 0.0002.The dot-dashed lines indicate the estimated value for ℓp based on Eq. ( 20) and dashed lines show the corresponding simplified profile (19).
Figure 8 shows typical simulation results of the velocity profile obtained in such a situation.From Fig. 8 we observe a smooth transition between a parabolic Poiseuille profile in the center and a more flat Darcy-Brinkman profile in the porous layers.These observations are in qualitative agreement with earlier simulations [23].In the center region, the magnetic field increases the effective viscosity, therefore slowing down the flow.In the porous region, however, the magnetic field is found to have a much smaller effect on the flow.
In a simplified description, we can approximate the velocity as constant, v = v 0 , within the porous layers of width ℓ p , and a parabolic profile in the center of the channel.Insisting on a continuous velocity profile at the interface between the two at y = ℓ p and y = L − ℓ p , we arrive at with L the channel width.The quantity u can be expressed in terms of the maximum velocity v max in the channel center as u The flow rate V = L 0 v x (y)dy for the flow profile ( 19) can be calculated as For constant v 0 < v max , increasing the width ℓ p of the porous layer leads to a proportional decrease in the flow rate.
In view of possible applications, it might be useful to be able to estimate the thickness of the porous layer ℓ p without the need to determine the detailed flow profile.We assume one can measure the total flow rate V and the centerline velocity v max .In our simulations, we obtain v max from fitting a parabola to the flow profile in the center region.Furthermore, we approximate the velocity within the porous layer by Darcy's law (1), v 0 = −(K/η)dp/dx = f ext /ξ.Knowing the channel width L, we can then determine the thickness of the porous layer within the simple model ( 19) by solving Eq. ( 20) for ℓ p .
Vertical dot-dashed lines in Fig. 8 indicate the value of ℓ p obtained in this way, while dashed lines show the corresponding approximate profile (19).Note that changing the magnetic field strength alters the flow profile, but the estimate for ℓ p remains unchanged within numerical accuracy.From Fig. 8, it is apparent that Darcy's law captures the mean velocity, but approximating the flow as constant within the porous layer is a rather crude approximation.Consequently, the layer thickness is underestimated by around 15-30% within the parameter range investigated.In spite of this inaccuracy, such a simple estimate might be a useful first step in determining the thickness of porous layers.

Conclusions
In this communication, we propose an extension of the hybrid MPC-BD scheme developed in Ref. [24] to describe fluctuating ferrohydrodynamics in porous media by additional friction forces.We performed numerous computer simulations and validated the model and its implementation in several ways and over a considerable range of parameters.In particular, we verified that the newly introduced friction coefficient is identical to the damping parameter in the Darcy-Brinkman model.Using an Irving-Kirkwood approach, we argue that this identity is expected, at least on large enough scales where hydrodynamics emerges from the MPC model.We also verified the expected dependence of the permeability on the friction coefficient.Therefore, the current MPC model can serve as an alternative to the Lattice-Boltzmann implementation proposed in Ref. [20] for the non-magnetic case.We mention that our implementation of the MPC method for non-magnetic fluids improves on an earlier approach [23], which suffered from a number of misconceptions concerning model parameters.
In the present work, we benefit from the flexibility of the MPC approach and couple the stochastic rotational motion of MNPs to the MPC scheme for flow through porous media.From simulations of channel flows, we find a "magneto-permeability effect" in porous media, i.e. the field-dependent change of the effective permeability.This effect is analogue to the traditional magnetoviscous effect in ferrofluids [3] and can be described theoretically in the same manner.
As an application of the method, we consider the flow through a planar channel with walls covered by a layer of porous material.In this situation, a parabolic velocity profile develops in the center of the channel, as is wellknown for Poiseuille flow of viscous fluids.On approaching the porous layers, the parabolic profile gives way to the more plug-like flow, which is typically observed in porous materials.We propose a simple method to estimate the thickness of the porous layer, based only on the total flow rate and the centerline velocity.Such rough estimates might be useful in a number of practical applications.
The present study can be extended in various ways.First, the extension to fully three-dimensional flows and more complicated geometries is straightforward, thanks to the flexibility of the MPC method.These extensions, together with a generalized model to describe non-dilute ferrofluids have already been explored for flow in non-porous media [34].From a more conceptual point of view, future studies are needed to investigate the role of porous media on the rotational motion of nanoparticles.Once these effects are better understood, they can then be included within the present hybrid MPC-BD scheme.A further extension of the present work could be the investigation of strongly heterogeneous, e.g.fractured porous media and to revisit earlier studies [16] with a more reliable constitutive model.The current approach is also well-suited to be incorporated within multi-scale simulation schemes [9,35], which are needed e.g. to study the highly interesting phenomenon of colloidal deposition in porous media [36,37].

Fig. 2
Fig.2The flow rate divided by the strength of the driving force, V/f ext , is shown versus f ext on a semi-logarithmic scale.Symbols denote simulation results obtained for parameters T = 0.1, ∆t = 0.2, Q = 100 and ξ = 0.005.Dashed line indicates the limiting value for weak driving forces.

Fig. 3 (
Fig. 3 (a) The damping parameter α in the Darcy-Brinkman equation (2) extracted from fits to the velocity profile is shown versus the friction coefficient ξ in the MPC model for non-magnetic fluids.Note that a double-logarithmic scale is used.Open black squares and filled blue circles correspond to ∆t = 1.0 and 0.2, respectively.The dashed line indicates the relation α = ξ.(b) Permeability parameter K versus MPC friction coefficient ξ on a doublelogarithmic scale.The same model parameters have been chosen and the same symbols are used as in (a).Dashed lines shows the power-law relation K = k 0 ξ −κ with exponent κ ≈ 0.99 and prefactor k 0 ≈ 0.13 for ∆t = 1.0 and κ ≈ 1.02 and k 0 ≈ 0.30 for ∆t = 0.2.

Fig. 4 (
Fig. 4 (a) The parameter r is shown parametrically versus α for different friction coefficients and non-magnetic fluids.The same color coding is used as in Fig. 3. Dashed lines are powerlaw fits r ∼ α b with exponents b ≈ 0.49 and 0.50 for ∆t = 1.0 and 0.2, respectively.(b) The kinematic viscosity ν extracted from fits to the velocity profile (3) versus the friction coefficient ξ in the MPC model.The same color coding is used as in Fig. 3.The dashed lines indicate the kinematic viscosity determined in the absence of a porous medium (ξ = 0).

Fig. 5
Fig. 5 The velocity profiles vx(y) across the channel are shown for ferrofluid flow through a porous medium with h = 0 (open squares) and h = 5 (filled circles) for ∆t = 0.2, n * = 0.005.Top curves correspond to ξ = 0.01, bottom ones to ξ = 0.02.Dashed lines show fits to the profile (3) in the Darcy-Brinkman model.

Fig. 6
Fig.6(a) The effective damping parameter α in the Darcy-Brinkman equation(2) extracted from fits to the velocity profile is shown versus the applied magnetic field strength h.The same model parameters are chosen as in Fig.5with ξ = 0.01 (lower) and ξ = 0.02 (upper).(b) Effective permeability parameter K eff versus magnetic field strength h.The same model parameters have been chosen and the same symbols are used as in (a).Dashed lines show the theoretical result (18) explained below.

Fig. 7 (
Fig. 7 (a) The relative change of the effective viscosity, ∆ν eff (h)/ν(0), versus the applied magnetic field strength h.The model parameters are chosen as in Fig. 6.Dashed line shows a fit to Eq. (17).(b) The permeability K scaled with the friction coefficient ξ is shown versus the effective viscosity ν eff for the same data as in panel (a) and in Fig. 6(b).The dashed line shows the relation ξK = ν eff .