Hydrodynamics of immiscible binary fluids with viscosity contrast: a multiparticle collision dynamics approach

Abstract

We present a multiparticle collision dynamics (MPC) implementation of layered immiscible fluids A and B of different shear viscosities separated by planar interfaces. The simulated flow profile for imposed steady shear motion and the time-dependent shear stress functions are in excellent agreement with our continuum hydrodynamics results for the composite fluid. The wave-vector dependent transverse velocity auto-correlation functions (TVAF) in the bulk-fluid regions of the layers decay exponentially, and agree with those of single-phase isotropic MPC fluids. In addition, we determine the hydrodynamic mobilities of an embedded colloidal sphere moving steadily parallel or transverse to a fluid–fluid interface, as functions of the distance from the interface. The obtained mobilities are in good agreement with hydrodynamic force multipoles calculations, for a no-slip sphere moving under creeping flow conditions near a clean, ideally flat interface. The proposed MPC fluid-layer model can be straightforwardly implemented, and it is computationally very efficient. Yet, owing to the spatial discretization inherent to the MPC method, the model can not reproduce all hydrodynamic features of an ideally flat interface between immiscible fluids.

---

1 Introduction

Multi-phase fluid flows occur ubiquitously in nature and engineering processes. Examples constitute oil–water flows, fluids with air bubbles, emulsions, dairy products, biological fluids, processing of paints, coating, and printing. Owing to their complexity, the theoretical description and efficient modeling of binary fluids pose major challenges, which stimulated a wealth of endeavors to model binary fluids using mesoscale simulations. In the lattice-Boltzmann method (LBM), the implementation of multi-phase flows and phase separation encompasses several variants: the color gradient model,^1–3 the pseudo-potential model,^4,5 the free-energy functional model,^6,7 and the mean-field model,⁸ or combinations thereof.⁹ Dissipative particle dynamics (DPD) simulations, which explicitly account for conservative pair interactions between fluid particles, allow to realize multi-phase fluids via assigning distinct interactions between the particles.¹⁰ Furthermore, the multiparticle collision dynamics (MPC) method, a particle-based hydrodynamic simulation approach which captures hydrodynamic interactions and thermal fluctuations,^11–17 has been proven valuable and efficient for mesoscale simulations, and has been applied in a broad range of studies of biological and active polymers,^18–26 colloids,^27–31 proteins,^32,33 vesicles and blood cells,^34,35 microswimmers,^36–46 and microfluidics.^47,48 To date, various MPC implementations of binary fluid mixtures have been proposed, and their phase behavior has been studied.^49–56 Depending on the applied interaction rule between the different fluid components, the viscosity values of the (two) fluids are equal or individually controlled.^52,55,57,58

Most of the above mentioned simulation methods are aimed to account for both the hydrodynamics and thermodynamics. The large computational costs, which are often necessary to suitably account for the thermodynamics involved in studying phase separation of multi-phase fluids, are dispensable when the dynamics of embedded objects such as proteins, polymers, or living organisms are considered. In fact, a plethora of physical phenomena related to immiscible binary fluids take place under conditions where phase separation is absent or is of no interest, and simulation methods accounting for the hydrodynamics alone suffice here. The MPC approach is very well suited to efficiently simulate hydrodynamic flow properties in the presence of thermal fluctuations (fluctuating hydrodynamics).¹⁴ In particular, MPC allows to tune the viscosity of fluids through the specification of the frequency of MPC collisions, and, hence, to control the viscous properties of immiscible fluids.

In this work, we present a model for planar layers of two immiscible binary fluids A and B using the MPC approach. The fluids, separated by a flat interface, are of distinct shear viscosity, η_A and η_B, whose values are tuned by the corresponding MPC collision frequency. While omitting the thermodynamic and kinetic processes of phase separation, it allows for fluid particle exchange across the interface, associated with a change of the local (collisional) interactions in the arriving fluid layer. No explicit interactions between fluid particles at the interfacial zone are required, although a more sophisticated modeling of the interface properties is possible for future assessment. Shear flow profiles and the shear stress under starting flow conditions are calculated, and the latter is compared with a provided analytical solution of the linearized Navier–Stokes equation of the same composite fluid. Moreover, transverse hydrodynamic velocity correlation functions (TVAFs) are determined for the different layers. In addition, the hydrodynamic mobility/friction properties of a colloidal sphere inside a fluid layer, which moves steadily parallel or perpendicular to an interface, are calculated. The approach recovers the correct flow profiles, fluctuating hydrodynamic properties, and thermal fluctuations of the individual fluid layers. The invoked simplifications in the present MPC treatment lead to a higher computational efficiency compared to other mesoscale techniques and MPC implementations considered so far,^51–56 which is a significant advantage when simulating large-scale systems.

The present simulations constitute a first important step in studying the dynamics, e.g., of monolayers of thermal particles moving near a planar fluid–fluid interface, with full account of the time-resolved (retarded) hydrodynamic interactions of the particles with the interface and among each other. These so-called quasi-two-dimensional systems have been intensely studied recently, since they reveal peculiar dynamic features such as the anomalous hydrodynamic enhancement of lateral collective diffusion,^59–62 and the influence of the interface on the motion of nearby Brownian particle, as reflected in the non-isotropic, hydrodynamic long-time tails of particle velocity correlations.^63,64 Interestingly enough, the motion pattern of microswimmers is also strongly affected by their hydrodynamic interaction with a nearby (fluid) interface.^65,66

Our two-fluids MPC model is also a first step toward mesoscopic simulations of the diffusion and phase behavior of assemblies of interacting proteins attached to or embedded inside a membrane. It should be recognized here that the biophysical properties of the membrane, both in physiological and in vitro conditions, influence the structure and function of many membrane-associated proteins.^67–73 Diffusion properties of single membrane receptor proteins and their orientation-dependent interaction potentials (which can be partially due to local membrane deformations) as obtained from force-field based molecular dynamics (MD) simulations where the lipids and the atomistic structure of the receptor are explicitly accounted for, can be used as input to tune mesoscopic MPC simulations.

The present paper is organized as follows. Section 2 gives the essentials of the single-phase MPC algorithm, outlines its extension to immiscible multi-phase fluids, describes the coupling rules of a colloid with the MPC fluid, and defines the simulation parameters. The two-fluids MPC model is validated in the three subsequent sections. In Section 3, the stationary shear profile of the planar three-layers system and the time-dependent shear stress functions under starting flow conditions are simulated, and compared with our analytic continuum hydrodynamics results. In Section 4, MPC simulated transverse velocity correlation functions (TVCFs) in the bulk regions of the two fluids, and in a region including the interface between them, are contrasted with predictions from the linearized fluctuating Landau–Lifshitz Navier–Stokes equation. In Section 5, the simulated hydrodynamic mobilities of a colloidal sphere moving steadily inside the middle layer of a three-layers fluid system are compared to previous numerical results based on the Stokes equation of low-Reynolds number hydrodynamics. In Section 6 we summarize and conclude our findings, and provide a perspective on future work. The Appendix presents our continuum hydrodynamics results for the time-dependent velocity profiles and for the stress functions of the composite three-layers system under starting shear flow conditions.

2 Model

2.1 Multiparticle collision dynamics (MPC) fluid

A single-phase MPC fluid consists of N point particles each of mass m, typically enclosed in a cubic simulation box of length L with periodic boundary conditions. The dynamics of the fluid particles proceeds through discrete streaming and collision steps.^11–13,74 During a streaming step, the particles move ballistically for a time span h, referred to as collision time. Hence, the position r_i of a fluid particle i, with i ∈ {1,…,N}, is updated according to

r_i(t + h) = r_i(t) + hv_i(t), (1)

where v_i is the particle velocity. In the subsequent collision step, accounting for the interactions between fluid particles, the MPC particles are sorted into cubic cells of size a defining the local interaction environment (collision cells). In the stochastic rotation dynamics (SRD) variant of MPC, MPC-SRD, the relative particle velocities ṽ_i = v_i − v_cm, with respect to the center-of-mass velocity v_cm of a particular collision cell, are rotated around a randomly orientated axis by a fixed angle α.^12,13,75 In three dimensions, the velocity of a particle i after a collision is thus given by

v_i(t + h) = v_cm(t) + ṽ_i,⊥(t)cos(α) (2)

+ [ṽ_i,⊥(t) × ℜ]sin(α) + ṽ_i,‖(t), (3)

where ℜ is a unit vector along the selected rotation axis, and ṽ_i,⊥ and ṽ_i,‖(t) are the parallel and perpendicular components of the velocity ṽ_i with respect to the vector ℜ, respectively. The orientation of ℜ is chosen independently for each collision cell and time step. The MPC-SRD algorithm conserves particle number, energy, and linear momentum. However, angular momentum is not conserved,¹⁷ an aspect irrelevant for the current study.^14,76,77 The employed discretization into collision cells breaks Galilean invariance, which is re-established by a random shift of the collision cell lattice at every collision step.⁷⁸ To maintain a constant temperature, a (simple) cell-level scaling scheme of the relative velocities ṽ_i is employed.⁷⁹ Since mass, momentum and energy are conserved locally, the correct fluctuating hydrodynamic equations for an isothermal compressible fluid are obtained in the continuum limit.

The shear viscosity η = η^k + η^c of a homogeneous MPC fluid can approximately be calculated analytically and comprises contributions from streaming, η^k, and collisions, η^c,^76,79–82 where

Here, 〈N_c〉 is the average number of MPC particles per collision cell, T the temperature, and k_B the Boltzmann constant. Note that the kinetic contribution, η^k, is only approximately valid, because its derivation is based on the molecular chaos assumption.⁷⁹ However, for small collision time steps, η is dominated by η^c.

The viscosity of a MPC fluid can be adjusted in several ways: by changing the mass of a particle, the average particle density in a collision cell, the rotation angle α, and the collision time step h. An implementation of two fluids with different masses has been realized.^57,58 Since the particles are distinctly different, suitable boundary conditions between the two fluids have to be applied. Similarly challenging are simulations of coexisting fluids with a MPC particle density difference. Variations of the collision angle yield a rather limited range of viscosity differences.⁸³ The most suitable strategy to simulate immiscible fluids seems to be a change of the collision time step in the different fluid regimes. In fact, the particles themselves are identical in the various fluids, they only experience more (or less) frequent collisions. This drastically simplifies the numerical implementation and enhances the performance.

2.2 Immiscible binary fluid system

The extension from a single to a layered two-fluid system of phases A and B with distinct shear viscosities η_A and η_B, which are separated by two flat interfaces, is rather straightforward and illustrated in Fig. 1. As discussed above, the viscosity of a single-phase MPC fluid is sensitive to the collision time h,^12,76 and we describe the fluid layers A and B of different viscosities by using accordingly different collision time steps h_A and h_B. Without loss of generality, we take h_A < h_B in the following, implying the viscosity of fluid A to be larger than that of fluid B. Furthermore, A and B particles are assumed to be of equal mass m, and the mean mass densities of both fluid phases are taken to be the same. Hence, the MPC particles in the two regimes are identical. For notational convenience, we will refer to particles A and B and fluids A and B corresponding to MPC particles in the two different domains.

Fig. 1 Schematics of a periodic three-layer MPC system of two immiscible fluids A (red) and B (blue) separated by planar intefaces. MPC collisions are performed independently in cubic collision cells, delineated by the dashed lines, using collision times h_A and h_B according to the fluid type. The width, a, of a cubic MPC cell is indicated.

As for a single-phase fluid, the MPC fluid particles move ballistically, and undergo independent collisions as follows:

• Streaming: A and B fluid particles move according to eqn (1) with collision time step h_A. This ensures an identical “continuous” dynamics of both particle types.

• Sorting: after streaming with time h_A, the particles are sorted in (shifted) collision cells.

• Rotation: particles in the collision lattice of the A domain undergo rotations according to eqn (3), after the time interval h_A. Particles in the collision lattice of the B domain are rotated only after the time interval h_B > h_A, i.e., after h_B/h_A additional streaming steps of time length h_A.

The subdivision of the streaming of particles B does not affect the properties of the bulk part of fluid B, since the particles move ballistically. However, it is important for the particles close to the planar A–B interface, owing to the random shift of the collision lattice normal to the interface, and since all particles in the (shifted) lattice of the A fluid domain undergo rotations after the time interval h_A. The equivalent streaming motion of particles in the two domains ensures an, on average, homogeneous distribution of fluid particles across the system, and a homogeneous density across the interface. However, the random shift of the collision cell lattice normal to an A–B interface broadens the interface at least by the size a of a collision cell.

2.3 Colloid dynamics and fluid coupling

The translational and rotational motions of a neutrally buoyant no-slip hard-sphere colloid of radius R and mass M embedded in the MPC fluid is governed by elastic collisions with the MPC particles, which we account for in a coarse-grained manner.¹⁷ During the streaming step, just as for the solvent particles, a colloidal sphere moves ballistically with center-of-mass velocity V_c(t). Its center-of-mass position vector, R_c(t), changes according to

R_c(t + h) = R_c(t) + hV_c(t). (6)

MPC particles i which (virtually) penetrate the colloid are moved backwards in time by the time interval (h − h_i), where h_i < h follows from the condition |r_i(t) − R_c(t) + h_i(v_i(t) − V_c(t))|² = R². These MPC particles collide then with a virtual colloid at the center position R_c(t) + h_iV_c(t), transfer the momentum p_i elastically to the colloid, and subsequently move with new velocity for the time interval (h − h_i). The linear and angular velocities of the MPC particles and of the colloid before and after collision are related byThe sum extends over all fluid particles colliding with the colloid during the time interval h. Here, Ω_c is the angular velocity of the embedded colloid, n_i = (r_i − R_c)/|r_i − R_c| is the unit vector pointing from the colloid center to the position of fluid particle i, and I = χMR² with χ = 2/5 is the moment of inertia of the spherical colloid.

To realize the hydrodynamic no-slip boundary conditions at the colloid surface, we use the bounce-back rule for the MPC fluid particles, which yields^17,84–87

with the relative velocity, , of a colliding MPC fluid particle i with respect to the according colloid surface point given byHere, μ_red = mM/(m + M) is the reduced mass, and and are the normal and tangential relative velocity parts, respectively, with respect to the colloid surface.

In the MPC collision step, phantom (p) particles are added inside the colloid to enforce the no-slip hydrodynamic boundary condition, which, in addition, act as a thermal bath.⁸⁸ Theses particles are uniformly distributed inside the colloid according to the average MPC–fluid particle density, and their velocities relative to the colloidal translational and rotational velocities are taken from a central Maxwellian distribution function. This yields the updated colloid translational and angular velocities after a collision step

respectively. Here, Δp^p_i denotes the change in the linear momentum of phantom particle i at position r^p_i due to SRD, and V_c(t + h) and Ω(t + h) are the velocities in eqn (8) and (9), respectively.

To further speed up the simulations, we use a common value h_i = h/2 for all MPC particles rather than considering the individual elastic collision events at the exact times t + h_i of each fluid particle-colloid collision. This simplifying step was shown to be as accurate as when the exact h_is are used, especially for small collision time steps.^84,86,89

2.4 Simulation parameters

In what follows, lengths are measured in units of a, mass in units of m, and energy in units of the thermal energy k_BT. We use therefore the unitsfor time, t₀, velocity, v₀, and viscosity, η₀, respectively. Note that t₀ is equal to the ratio of cell size a and thermal velocity . In these units, the sound velocity in both fluids is equal to one. The average number of particles per collision cell is selected as 〈N_c〉 = 10, implying equal mean number and mass densities of the A and B fluids, and the rotation angle is set to α = 130°. The collision time steps are taken as h_A = h_B/5 = 0.02 × t₀. With the mass density ρ = 〈N_c〉m/a³, the corresponding kinematic viscosities are ν_A/ν₀ = 4.12 and ν_B/ν₀ = 0.87, where ν₀ = η₀a³/m. This yields the kinematic viscosity ratio μ² = ν_A/ν_B = 4.74. The related dimensionless Schmidt numbers are Sc_A = ν_A/D_A = 400 and Sc_B = ν_B/D_B = 17, expressing that the viscous diffusion of (transversal) momentum in the fluid is distinctly faster than diffusive mass transport, with the latter characterized by the mass diffusion coefficients D_A = (h_A/h_B)D_B of fluid A and B particles, respectively. Simulations are performed using periodic boundary conditions, applied in Sections 3 and 4 to a cubic simulation box of length L/a = 39 and 80, respectively, and in Section 5 to a rectangular box of lengths 2L_x/a = 2L_y/a = L_z/a = 80. The latter embeds a colloidal sphere of radius R = 2.5a.

3 Shear simulations

As a first example used for scrutinizing the hydrodynamic behavior of our two-fluids MPC approach, we consider a standard stationary shear flow setup as sketched in Fig. 2. The three planar layers of two immiscible fluids A and B are sheared by two walls oriented parallel to the xy-plane, which move oppositely along the x-direction with constant velocities ±u = 0.0975v_th. The lower wall is located at z = 0, and the upper one at z = L_z = L = 39a. No-slip boundary conditions (BCs) at the walls are implemented using the bounce-back rule and phantom particles inside the walls.⁷⁶ The three B–A–B fluid layers are separated by planar fluid interfaces located at z = L_z/4 and 3L_z/4, respectively.

Fig. 2 Schematics of a layered B–A–B fluid system steadily sheared by two parallel no-slip walls moving in opposite direction with velocities ±u = (±u,0,0). A piece-wise linear stationary velocity profile is obtained from hydrodynamics under laminar flow conditions.

The stationary shear velocity v^st(z) = v^st_x(z)e_x, obtained from the stationary Navier–Stokes equation,⁹⁰ is piecewise linear and unidirectional along the x-direction. The flow is uniquely determined by the wall-fluid stick boundary conditions, and the continuity of flow velocity and shear stress across the two clean planar interfaces whose thickness is assumed to be zero. Explicitly,

v^B,st_x(L_z) = −v^B,st_x(0) = u, (15)

v^A,st_x(L_z/4) = v^B,st_x(L_z/4) = u⁻_s, (16)

v^A,st_x(3L_z/4) = v^B,st_x(3L_z/4) = u⁺_s, (17)

where u⁺_s and u⁻_s = −u⁺_s (by symmetry) are the fluid velocities at the upper and lower fluid interfaces, respectively. The interfacial velocities are obtained using the continuity of shear stress across the planar clean interfaces (no Marangoni stress and Laplace pressure), η_A[small gamma, Greek, dot above]_A = η_B[small gamma, Greek, dot above]_B, at z = L_z/4 and 3L_z/4, respectively, which yieldsHere, [small gamma, Greek, dot above]_A = du^A,st_x(z)/dz and [small gamma, Greek, dot above]_B = du^B,st_x(z)/dz are the constant shear rates in the mid-layer of fluid A and in the two layers of fluid B, respectively. The (dynamic) pressure for the unidirectional shear flow is constant throughout the system including the interfaces. The Appendix presents our analytical continuum hydrodynamic result for the transient starting flow v_x(z,t) of the B–A–B system which converges to v^st_x(z) in the course of time (cf.Fig. 6).

The MPC simulation results for v^st_x(z) displayed in Fig. 3(a) reflect the hydrodynamically expected behavior of three linear stationary shear flow regions. Even more, the simulation results agree quantitatively with the hydrodynamic flow profile described in eqn (17), (18), and in eqn (36) of the Appendix. In spite of the non-zero thickness of the interface in the MPC simulations of the order of the collision cell size a, caused by discretization in terms of collision cells and random shift of the collision cell lattice, the MPC results suggest that the interface width is of minor relevance for fluid properties on lengths scales significantly larger than a. The inset of Fig. 3(a) magnifies the stationary velocity profile in the A–B interfacial region. It suggests a continuous change both of v^st_x(z) and its slope across the interface. This indicates also a continuous change of the local viscosity in the interfacial region caused by the discretization.

Fig. 3 (a) Stationary shear velocity profile, v^st_x(z), of a B–A–B fluid layer system, obtained from MPC simulations (open circles) and analytically from continuum hydrodynamics (solid line) according to eqn (17), (18) and (36) in the Appendix. The magnitude of the wall velocity is |u| = 0.0975v_th. Inset: Magnification of v^st_x(z) in the A–B interfacial region (blue rectangle). (b) Moving time average of the external shear stress, 〈σ^e_xz〉_t, at the upper and lower wall, and of the internal shear stress, 〈σⁱ_xz〉_t, in units of thermal stress σ₀ = k_BT/a³ = 13σ^st_xz. Open symbols are MPC data, and solid lines our continuum hydrodynamics predictions (cf. Appendix). Inset: Continuum hydrodynamics results for the instantaneous (solid lines) and moving time average (eqn (22)) (dashed lines) external (red) and internal (green) stresses σ^e,i_xz(t) and 〈σ^e,i_xz〉_t, respectively, as functions of t (as in the main plot).

In order to confirm the hydrodynamically expected relations between stress, shear rate, and viscosity in the MPC simulations of the immiscible fluid, we determine the instantaneous internal (superscript i), σⁱ_xz, and external (superscript e), σ^e_xz, shear stresses. The latter is the stress (modulus) exerted by the fluid particles on the walls, while the former one is the volume averaged stress. Explicitly, the stresses at a given time instant are computed as⁷⁶

where the sums extend over all N fluid particles inside the simulation box, and r_iz is the position of the MPC particle i along the z-axis. Here, the change in the momentum, Δp_i(t), of a particle i in a collision step is given by

Δp_i(t) = m(v_i(t) − v̂_i(t)), (21)

where v̂_i is the particle velocity after streaming and before collision. The superscripts u and l indicate that the considered quantity is calculated at the upper and lower wall, respectively. Note that eqn (19) and (20) account also for momentum exchange due to collisions with phantom particles located inside wall boundary cells (bc). The negative sign in front of the transversal momentum exchange Δp^l_ix accounts for the negative velocity, −u, of the lower wall. In eqn (20), the second term on the right-hand side accounts for the momentum change of B particles “reflected” (bounce-back) at a wall in the streaming step, and Δt_i is the time during which particle i streams in the fluid before colliding with a wall.⁷⁶ The internal stress calculation invokes the momentum exchange of fluid A and B particles, described by the third term on the rhs of eqn (20). Here, time averaging is performed separately for each fluid phase, owing to the different collision times h_A and h_B.

In the simulations, “macroscopic” stress tensors are calculated via averaging over various realizations (ensemble average) as well as averaging over time. For the latter, we determine the moving time averages

of external and internal stresses, which yield the stationary state value, σ^st_xz, in the limit t → ∞. As a particular case, we analyze the starting flow situation, where at time t = 0, the two confining walls suddenly start to move oppositely with constant velocities ±ue_x, respectively. Fluid and walls are at rest for t < 0.

MPC results for the moving time average external and internal shear stresses are presented in Fig. 3(b), as function of the elapsed time t after the two walls started to move. The MPC results for the modulus of the external stress, 〈σ^e_xz〉_t, at the upper and lower wall are equal within the accuracy of the simulations. The external stress decays monotonically towards the plateau value σ^st_xz = 0.077σ₀ where the steady-state regime is reached, characterized by the fully developed, piece-wise linear shear profile in Fig. 3(a) and an uniform shear stress. Its decay reflects the diffusive broadening of the region of changing fluid velocity near the walls with increasing time, as depicted in Fig. 6 of the Appendix. The moving time averaged internal stress increases instead from its minimal value at t = 0, where the bulk fluid is still at rest, towards its steady-state value σ^st_xz. The characteristic transition time for the external shear stress relaxation (and internal stress buildup) towards the uniform steady-state value is estimated as τ_v = (L_z/4)²(1/ν_A + 1/ν_B) ≈ 140 × t₀, which is the viscous diffusion time across half of the simulation box, L_z/2. The MPC moving time averaged external and internal stresses in the sheared B–A–B system are in excellent agreement with the corresponding hydrodynamic stresses (solid curves), the latter obtained from the according analytic expressions presented in the Appendix.

From the limiting steady-state stress value and the steady-state flow profile, the viscosity values η_A/η₀ = 42.9 and η_B/η₀ = 9.1 are deduced, which agree within less than 5% with the viscosity values obtained from analytical theory for respective one-component MPC fluids (cf. Section 2.4). The shear viscosities of the binary fluid model can be easily controlled by a single parameter, namely the collision time step h, but there is a continuous viscosity crossover along the MPC interface of thickness comparable to the collision cell. In general, the continuum hydrodynamic behavior is accurately recovered by the MPC simulations for lengths larger than about 2a.¹⁴

The inset in Fig. 3(b) depicts our continuum mechanics results, summarized in the Appendix, for the instantaneous external and internal shear stresses σ^e,i_xz(t) in comparison with the according moving time averages 〈σ^e,i_xz〉_t. By their definitions, the moving stress averages approach the common steady-state value more slowly than the stresses themselves. As shown in the Appendix, σ^e_xz(t) and 〈σ^e_xz〉_t exhibit the time-dependence near t = 0. In contrast, the internal stress and its time averaged counterpart are finite and minimal at t = 0, with the common value 0.60σ^st_xz according to eqn (43).

4 Hydrodynamic correlations: transverse velocity auto-correlation function

Additional insight into the time-resolved hydrodynamic behavior of the MPC fluid is gained by analyzing the transverse velocity auto-correlation function (TVCF)¹⁴ in the various layers. For a stationary and isotropic Newtonian fluid in a volume V with periodic boundary conditions, the linearized Landau–Lifshitz Navier–Stokes equations yields the single-exponentially decaying TVCF in Fourier space¹⁴Here, u^T(k,t) is the Fourier-transformed velocity part perpendicular to the wave vector k,^14,82i.e., u^T·k = 0. The brackets denote an equilibrium ensemble average, with the fluid system at rest on hydrodynamic time and length scales. The factor 2 on the rhs accounts for the two independent transversal modes. Owing to isotropy, the TVCF depends only on the modulus k = |k| of the wave vector. Simulation results for the TVCF of a single-phase MPC fluid are in excellent agreement with the above hydrodynamic prediction.^14,82,91

To explore thermally induced transverse velocity correlations in our three-layer model of fluids A and B, we perform simulations for a cubic simulation box of size L = 80a, with periodic boundary conditions in all three Cartesian coordinate directions. The higher-viscosity layer A of width L/2 is symmetrically sandwiched between two fluid–B layers, as illustrated in Fig. 2, but now for a system without shear. The period boundary condition along the z-axis implies an alternating pattern of horizontal A and B layers of equal thickness L/2.

We determine the TVCFs of the pure A and B fluids in the three-layer model by considering N_o MPC particles inside an observation cuboid of z-thickness L_o = L/4 = 20a and volume V_o = L × L × L_o, symmetrically located inside the A-fluid and B-fluid layers, respectively. To explore additionally the influence of the fluid interface, the TVCF for another observation cuboid with smaller vertical width L_o = 10a is determined, with the cuboid symmetrically enclosing the A–B interface. The Fourier transform, u(k,t), of the fluid velocity fluctuations in an observation cuboid is calculated according to

where v_i(t) is the velocity of fluid particle i at position r_i(t) inside the considered cuboid. For the cuboid centered around the A–B interface, half of the particles summed over are, on average, of A-type and half of B-type. For the present purpose, we consider only wave vectors k = k_‖ parallel to the xy-plane, with wavelength λ = 2π/k smaller than the cuboid width L_o. This reduces boundary artifacts due to fluid particles leaving or entering the observation cuboid.

Fig. 4(a) displays the normalized TVCFs

of the fluid inside the A and B cuboids, respectively, as well as the TVCF of the mixed-fluid cuboid enclosing the A–B interface. The horizontally oriented wavectors employed here are k = (32π/L)(1,0,0) and k = (32π/L)(0,1,0), of wavelength λ = 5a smaller than the cuboid width. Notice that u^T(k,t) = (1 − k⊗k/k²)u(k,t). Owing to the non-isotropic three-layer structure, C^T_v(k,t) is in principle an anisotropic function of the wave vector, depending also on the vertical location and width of the considered cuboid.

Fig. 4 (a) Normalized TVCFs, C^T_v(k,t), for a wave vector k parallel to xy-plane with wavelength λ = 5a obtained from MPC simulations for the fluid A (red, up triangles), B (blue, down triangles), and the A–B interface cuboid (green, circles). The corresponding solid lines represent the bulk-fluid prediction exp(−ν_A,Bk²t) and the double-exponential expression in eqn (26) with weights ξ_A = 0.541 = 1 − ξ_B (green solid line) and ξ_A = ξ_B = 1/2 (orange dashed line), respectively. Inset: TVCF dependence on ν_αk²t for α ∈ {A,B,int} and ν_int = ν_Aν_B/[ξ_Aν_A + ξ_Bν_B]. (b) Time-integrated normalized TVCFs, T(k,t) (27). Solid lines represent eqn (27) and (28), respectively. Inset: Data collapse for ν_αk²t.

Within the correlation time window t ≤ 3t₀ depicted in Fig. 4(a), the MPC-calculated normalized TVCFs of the pure A and B-fluid cuboids (open symbols) decay exponentially according to exp(−k²νt), with kinematic viscosity values ν_A and ν_B as numerically obtained in Section 3. The reason why the isotropic bulk fluid TVCF form is recovered in the anisotropic three-layer system (within numerical accuracy) is that the viscous diffusion time, τ^A,B_ν = (L/8)²/ν_A,B, over a distance from the cuboid center to the interface is large compared to the resolved correlation time window; the viscous diffusion times are τ^A_ν = 24t₀ and τ^B_ν = 115t₀, respectively. Hence, in the considered time window and the considered wave vector of wavelength λ = 5a, the velocity correlations in the single-fluid cuboids are yet unperturbed by the interfaces.

On the same basis one could expect that the MPC data for the TVCF of the cuboid symmetrically enclosing the A–B interface are for t ≤ 3t₀ decently well reproduced by the superposition of two bulk-fluid exponential TVCFs according to

C^T,int_v(k,t) = ξ_Ae^−ν_Ak2t + ξ_Be^−ν_Bk2t, (26)

for equal weight factors ξ_A = ξ_B = 1/2, and ν_A and ν_B determined from the MPC simulation data for the single-fluid layers. The equal-weight superposition according to eqn (26) somewhat underestimates the decay of the correlation function for t > t₀, reflecting the growing influence of the interfacial region with increasing time. A fit of the simulation data by eqn (26), for unchanged values of ν_A and ν_B, yields the weight factors ξ_A = 0.541 and ξ_B = 1 − ξ_A = 0.469 (dark-green solid line). The asymmetry could be a consequence of the shorter viscous diffusion time across the half-width L_o/2 for the fluid- A part of the two-fluid observation cuboid.

The time integral of the normalized TVCF (25), characterizing a one-component fluid in the hydrodynamic regime, is

where T(k,t), in the limit t → ∞, is related to the Oseen tensor in reciprocal space.^14,90,92 The accordingly time-integrated TVCF in eqn (26) for the cuboid enclosing the interfacial region iswith T^int(k,∞) = 1/k²ν_int and ν_int = ν_Aν_B/[ξ_Aν_A + ξ_Bν_B].

The time dependence of T(k,t) and T^int(k,t) for the three observation cuboids, obtained from the data of Fig. 4(a), are shown in Fig. 4(b). The time-integrated MPC simulation data (open symbols) agree overall well with the analytic expressions in eqn (27) and (28) based on the single-fluid theoretical expressions. As shown in the inset, the time-integrated TVCFs are universal functions of ν_αk²t, as expected by the identical universal behavior of the TVCFs. The factor ν_int = ν_Aν_B/[ξ_Aν_B + ξ_Bν_A] can be considered as a common effective kinematic viscosity of the A and B fluid contributions in the cuboid enclosing the interface. The inset further illustrates that the crossover to the long-time plateau values 1/(k²ν_α) is characterized by the viscous diffusion times τ^α_k = (ν_αk²)⁻¹ = 0.63a²/ν_α. Since τ^A_k < τ^int_k < τ^B_k, the Stokesian regime of inertia-free, quasi-instantaneous hydrodynamics is reached for the considered wavenumber at times distinctly smaller than the viscous diffusion time across the colloid diameter (2R)²/ν_B = 29t₀ ≫ τ^B_k.

5 Mobility of a colloidal sphere near a fluid–fluid interface

The (strong) viscosity difference between two immiscible fluid phases affects the dynamics of a colloidal particle moving near the fluid interface. To explore the hydrodynamic coupling, and to scrutinize the according MPC coupling predictions in our three-layer model, we calculate the mobility coefficients of a sphere embedded in fluid A which moves steadily under low-Reynolds-number conditions parallel or perpendicular to the planar A–B interfaces. The coefficients are determined as functions of the reduced distance d_z = z/(2R) of the sphere center from the A–B interface (see Fig. 5). To reduce finite-size effects due to the periodic boundary conditions in z-direction, different from Section 4, we consider a non-cubic simulation box of lengths 2L_x = 2L_y = L_z = 80a.

Fig. 5 (a and b) Schematics of the three-layer system for determining the lateral and transverse translational mobility coefficients, Γ_‖ and Γ_⊥, of a no-slip sphere of radius R = 2.5a embedded in fluid A, as functions of the reduced distance d_z = z/(2R) from the sphere center to the A–B interface. The employed viscosity ratio is η_B/η_A = 0.21. (c and d) Lateral and transverse mobility coefficients from MPC simulations (symbols) and hydrodynamic force-multipole-expansion calculations (green lines)⁹³ in units of the bulk mobility value 1/(6πη_AR).

As indicated in Fig. 5, the sphere of radius R = 2.5a is subjected to a weak constant force F_‖ (F_⊥) applied to its center, and oriented parallel (perpendicular) to the fluid–fluid interface. Due to the no-slip boundary conditions employed on the sphere surface, the moving sphere drags nearby fluid along, which is compensated by fluid backflow such that the total momentum of the system in any spatial direction is zero (quiescent fluid system assumed).^29,94 Under low-Reynolds-conditions, where in the continuum mechanics picture the fluid flow is described by the quasi-stationary linear Stokes equation, the reduced translational mobilities follow from the relations

by measuring the steady-state mean velocity 〈v_‖,⊥〉 of the sphere for a given constant force F_‖,⊥. In our MPC simulations, the thermal force value F_‖,⊥ = 4k_BT/a is used. After applying the force to the sphere, the steady-state with constant mean drift velocity is reached for times t ≫ R²/ν_A. Note that Γ_‖,⊥ = 1 in the bulk region of the fluid far from any interface or boundary.⁹⁰

Fig. 5(c) and (d) display the MPC results (open blue symbols) for the normalized lateral and transverse mobilities, Γ_‖(d_z) and Γ_⊥(d_z), as functions of the reduced distance d_z = z/(2R). For comparison, according reduced mobilities are shown (green solid lines) as obtained numerically using an elaborate Stokesian dynamics-based hydrodynamic force-multipole expansion method, encoded in the software package HYDROMULTIPOLE.⁹³ The depicted mobility curves by the force multipoles method, valid under creeping-flow conditions, are taken from⁹³ and constitute accurate continuum hydrodynamics results for a no-slip sphere in a half-infinite Newtonian fluid A which moves steadily parallel or perpendicular to an ideally flat and clean interface of zero interfacial viscosity and Marangoni stress. The interface separates the fluid- A half-space from the fluid–B half-space. Note that Γ_‖,⊥ depend on the ratio of the shear viscosities of the two fluids.

Both the MPC simulation and continuum hydrodynamic results predict the lateral sphere mobility to increase with decreasing distance d_z from the A–B interface. They are in good overall agreement, except for small distances where the simulation data are somewhat larger (see Fig. 5(c)). Regarding the transverse mobility depicted in Fig. 5(d), the continuum hydrodynamics curve for Γ_⊥ decreases strongly with decreasing distance, and assumes the value Γ_⊥ = 0 at the sphere-interface contact distance d_z = 0.5 due to lubrication. In contrast, while the MPC simulation data in Fig. 5(d) are in accord with a mild decline of the mobility for deceasing distance d_z ≳ 1.5, they do not reproduce the strong drop in Γ_⊥ at small distances d_z ≲ 1 (i.e., z ≲ 5a). On first sight, this discrepancy is surprising, since friction and lubrication effects for a hard-sphere colloid embedded in a single MPC fluid close to a solid no-slip wall are well reproduced.⁹⁵ However, it can be attributed qualitatively to the mixing of the two fluids in the interfacial region over a thickness larger than a collision cell size a, and to a local perturbation of the hydrodynamic flow field by the no-slip sphere moving normally to the nearby interface. In the HYDROMULTIPOLE calculations, the two fluid half-spaces are taken as ideally incompressible, and the interface as ideally thin and flat, without any sphere-induced perturbation. Notice further that the sphere size is comparable with the MPC interfacial thickness. At any rate, the MPC implementation of immiscible fluids captures the dynamics of the immersed colloidal sphere overall quite well.

6 Summary and conclusions

In this article, we have presented a MPC-based mesoscale hydrodynamic simulation scheme for modeling immiscible (layered) binary fluids with viscosity contrast separated by a flat interface.

Shear flow, external and internal shear stress, fluctuating hydrodynamic velocity correlations, and the hydrodynamic mobilities of an embedded spherical particle moving close to a flat fluid interface have been analyzed for a three-layer MPC fluid, and validated against continuum hydrodynamics predictions. We obtained a piece-wise linear stationary fluid velocity profile in excellent agreement with the continuum hydrodynamics prediction. By computing the shear stress in relation to the shear rate, we confirmed that the analytically obtained viscosity values for single-phase MPC fluids are reproduced by the binary fluid model, in regions distant from the fluid–fluid interface. Considering the build-up of the shear profile from a resting fluid, we obtain excellent agreement between MPC simulation results and hydrodynamic predictions for the moving time averages of the external and internal shear stresses, in the time range assessed in the simulations. For this comparison, we have derived an analytical solution of the linearized Navier–Stokes equation for the shear stress functions of the B–A–B system under starting flow conditions, using the Laplace transformation technique. The analytic expressions have allowed us to assess quantitatively the differences between instantaneous and moving time averaged stress functions.

To examine the predictions by our two-fluids MPC model regarding time-dependent correlations of thermally induced velocity fluctuations, we calculated the transverse velocity auto-correlation function (TVCF) in different observation cuboids. We showed that the TVCFs for the single-fluid cuboids follow closely the expected exponential decay, characterized by the kinematic viscosity of the respective fluid and the considered wave number. In contrast, the calculated TVCF of the cuboid enclosing the A–B interface is overall well fitted by a linear combination of the exponential TVCFs for bulk fluids A and B, using the viscosity values determined in our shear-flow studies. The approximate validity of linear superposition suggests that the TVCF of the cuboid is only mildly affected by the interfacial region. A stronger interfacial influence can be expected for a narrower cuboid of width smaller than the employed value L_o = 10a.

Finally, we have probed the hydrodynamic coupling of a steadily moving no-slip sphere to a nearby flat two-fluids interface by determining its hydrodynamic mobilities. The distance dependence of the lateral mobility coefficient for the three-layers MPC model agrees well with the according mobility result by a hydrodynamic force-multipole expansion method for a sphere moving close to an ideally flat, clean interface separating two incompressible fluids. While decent agreement is observed also regarding the transverse mobility for sphere-interface distances larger than three times the sphere radius, the sharp mobility decline at small distances predicted by the continuum hydrodynamics approach for a non-deformable planar interface of zero thickness, is not obtained by the MPC simulations. We attribute this to the mixing of the two fluids in the MPC interfacial region of thickness larger than the collision cell size a, and possibly a local perturbation of the interface caused by the transverse motion of the sphere. Moreover, and different from what is assumed in the force multipoles calculation, the two fluids in the MPC model are compressible. The non-zero compressibility of the fluids may play a role in particular for transverse (i.e., squeezing) sphere motions. To reduce the influence of the finite interface width on the sphere mobilities, a significantly larger sphere can be considered. Moreover, an alternative method to determine the mobilities may be better suited⁹⁵ in order to reduce fluid perturbations by the translating sphere. Here, further studies will be performed in the future.

A numerical advantage of the two-fluids MPC model is that the desired viscosities of the fluid phases can be easily prescribed using the analytic viscosity expression for a single-phase MPC fluid.¹⁴ Compared to other mesoscale simulation models of immiscible binary fluids, the present model is straightforwardly implemented, since it does not involve the computation of thermodynamic properties and kinetic processes related to phase separation. Hence, the computational cost is comparable to simulating two single-phase MPC fluids with different collision times.

The two-fluids MPC simulation method can be applied to a wide range of biological soft matter systems. For example, the approach can be suitably extended to study interfacial rheological properties including interfacial viscosity^96,97 and interfacial tension. Furthermore, as noted already in the introduction, the model can be applied to investigate the lateral self- and collective diffusion of different in-membrane or membrane-attached proteins. The effects of the viscosity contrast between a membrane and the adjacent cytosol, and hydrodynamic interactions between proteins and membrane, and among the proteins, on protein diffusion can be simulated over several timescales using a simplifying coarse-graining of the system. In this context, we will perform further simulations, specifically of colloids with a diameter comparable to the width of the fluid A domain. In a more refined analysis, lipid molecules and other macromolecules forming the membrane constitute a crowded environment which slows down the diffusion of embedded proteins.⁹⁸ Molecular crowding effects cause so-called sub-diffusion, identified recently to play a vital role in many biological phenomena,^99–102 including neuronal signaling.^103,104 For future assessment, molecular crowding mechanism can be implemented into our three-layer model, with the middle layer playing the role of the membrane, by adding a planar layer of interacting host particles to the middle layer, or alternatively and more realistically, by accounting for visco-elastic effects in the middle layer through semi-atomistic memory function calculations. Work by us in both directions is in progress.

Furthermore, the present MPC model can be extend to interfaces with imposed sinusoidal fluctuations mimicking membrane fluctuations. Moreover, the quantitative control of viscosity values opens the possibility to study systems with designed viscosity gradients. This provides a means to study the dynamics of biological macromolecules or microorganisms responding to viscosity gradients, such as in viscotaxis.¹⁰⁵

Conflicts of interest

There are no conflicts to declare.

Appendix

We present here the analytic hydrodynamic expressions for the transient (starting flow) shear stress and velocity profiles for the layered two-fluid system depicted in Fig. 1, valid under laminar flow conditions where the linearized Navier–Stokes equation of continuum hydrodynamics applies to.

Consider the horizontal planar walls at z = ±L_z/2 to be quasi-instantaneously set into motion at time t = 0, with constant velocities ±ue_x, respectively. The B–A–B layered composite fluid, originally at rest for t < 0, is then gradually set into laminar motion, through the diffusive transport of momentum (viscous stress) away from the moving walls into the fluid interior. The accordingly unidirectional starting flow velocity field, v(z,t) = v_x(z,t)e_x, is incompressible. Since the fluid pressure remains spatially uniform, the linearized Navier–Stokes equation for the starting flow reduces to the one-dimensional transversal momentum diffusion equations

for fluids A and B, respectively. For convenience, z = 0 denotes throughout the Appendix the midplane between the two walls in Fig. 1, so that v_x(−z,t) = −v_x(z,t) by symmetry. The initial and boundary conditions are hereexpressing that fluid A sticks to the walls, and that velocity and shear stress are changing continuously across the A–B fluid interfaces at z = ±L_z/4.

We have solved this linear boundary value problem using the time Laplace transform method and the method of residues, to obtain the composite velocity fields of A and B in the time domain. Skipping the lengthy derivation, we directly present the infinite series solution

v_x(z̄,τ) = v^A_x(z̄,τ)χ_A(z̄) + v^B_x(z̄,τ)χ_B(z̄), (33)

withand denominator functionHere, μ² = ν_A/ν_B = η_A/η_B, since we have assumed equal mass densities of both fluids such as in the MPC simulations. Moreover, z̄ = 2z/L_z is the reduced vertical distance and τ = 4tν_B/L_z² the reduced time. The characteristic functions for the considered slabs [0,L_z/4] and [L_z/4,L_z/2] of A of B fluids, respectively, are χ_A(z̄) = Θ(1/2 − z̄)Θ(z̄) and χ_B(z̄) = Θ(1 − z̄)Θ(z̄ − 1/2), respectively, with Θ denoting the unit step function.

The piecewise linear, long-time stationary velocity field has the fluid A and B contributions

The infinitely growing sequence of relaxation values, 0 < w₁ < w₂ < …, are the purely simple, positive roots of the transcendental equation (for an associated heat conduction problem see ref. 107)associated with the purely imaginary, pairwise conjugate simple poles, {±iw_k}, of the Laplace transformed velocity field ṽ_x(z̄,s) constituting a meromorphic function in the complex s-plane. The exponentially decaying modes in the series solution reflect the overdamped dynamics described by the momentum diffusion eqn (31). In the limiting case μ = 1 of equal kinematic viscosities, w_k = πk with integer values of k. For μ close to or large compared to one, approximate expressions for w_k(μ) are obtained using perturbation theory. However, these are not useful for the intermediate viscosity ratio μ² = 4.74 used in our MPC simulations. For arbitrary μ, we conveniently determined the roots of eqn (36) using Mathematica.

Fig. 6 depicts the starting flow velocity profiles according to eqn (33) and (34), for μ² = 4.74 and values of τ as indicated. For the smallest considered time, 100 terms in the sum over exponentially decaying modes have been accounted for. The stationary profile is reached for τ > 0.2.

Fig. 6 Starting flow velocity profiles of the sheared B–A–B system, obtained from eqn (33)–(35) for the viscosity ratio μ² = 4.74, and for reduced time values τ as indicated.

The shear stress field follows from the velocity field using σ^A,B_xz(z,t) = η_A,Bdu^A,B_x(z,t)/dz. The external stress exerted by the walls on the neighboring fluid layer B is obtained in particular as

whereis the spatially uniform stationary long-time stress in the sheared composite system. The internal hydrodynamic stress is the spatially averaged shear stress, i.e.,which yieldsAs noted in Section 3, in lieu of σ^i,e_xz(τ), we use the moving time averages in the analysis of the MPC data to smooth out statistical errors (see eqn (22)). The respective analytical stress expressions follow from eqn (41) and (22) by the replacementThe moving time-averaged stress converges more slowly toward the stationary stress value σ^st_xz than the stress itself. Both σ^e_xz(τ) and its moving time average decay monotonically in time toward σ^st_xz, which reflects the diffusive broadening of the zone of changing fluid velocity with increasing time (see Fig. 6). At very short times, where the influences of the fluid interfaces and the opposite wall are still negligible, the velocity decays steeply near the wall, from value u at the wall to 0 slightly off the wall. This implies the approximate inverse square-root time dependence of the external stress . The inverse square-root short-time divergence is shared by the moving time average stress, except with a twice as large amplitude for the latter. Different from the external stress, the internal stress is finite at τ = 0 where it attains its minimal valueFor equal viscosities of fluids A and B, where the hydrodynamic effect of the A–B fluid interfaces is absent, the intrinsic stress is constant at all times and equal to the stationary stress.

In the MPC simulation results for the moving time average stress depicted in Fig. 3(b) for μ² = 4.74, time is measured in units of t₀ = a²/ν₀ instead of the reduced time employed here where . Using that L_z = 39a in our MPC study of shear flow, the conversion relation (t/t₀) × 10⁻³ ≈ 0.44τ is obtained. We recall furthermore the relation σ^st_xz = 0.077σ₀ between stationary stress and MPC thermal stress σ₀ = k_BT/a³.

Acknowledgements

ZT thanks Marisol Ripoll (Forschungszentrum Jülich) for valuable discussions. The authors gratefully acknowledge computing time granted through JARA-HPC on the supercomputer JURECA at Forschungszentrum Jülich.¹⁰⁶ This work is embedded in the joint IBI-INM project on mesoscale modeling of intracellular signaling events underlying neuron functions.

Notes and references

1. D. H. Rothman and J. M. Keller, J. Stat. Phys., 1988, 52, 1119.
2. A. K. Gunstensen, D. H. Rothman, S. Zaleski and G. Zanetti, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 43, 4320.
3. Y. Yu, H. Liu, D. Liang and Y. Zhang, Phys. Fluids, 2019, 31, 012108.
4. X. Shan and H. Chen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 1815.
5. X. Shan and H. Chen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1994, 49, 2941.
6. E. Orlandini, M. R. Swift and J. M. Yeomans, EPL, 1995, 32, 463.
7. M. R. Swift, E. Orlandini, W. R. Osborn and J. M. Yeomans, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 5041.
8. X. He, X. Shan and G. D. Doolen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, R13.
9. J. Lombard, I. Pagonabarraga and E. Corvera Poiré, Phys. Rev. Fluids, 2020, 5, 064201.
10. P. V. Coveney and P. Español, J. Phys. A: Math. Theor., 1997, 30, 779.
11. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605.
12. G. Gompper, T. Ihle, D. Kroll and R. Winkler, Adv. Polym. Sci., 2009, 221, 1.
13. R. Kapral, Adv. Polym. Sci., 2008, 140, 89.
14. C.-C. Huang, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 056711.
15. C.-C. Huang, G. Gompper and R. G. Winkler, J. Chem. Phys., 2013, 138, 144902.
16. M. Yang, M. Theers, J. Hu, G. Gompper, R. G. Winkler and M. Ripoll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 013301.
17. M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Phys. Rev. E, 2016, 93, 032604.
18. M. Liebetreu, M. Ripoll and C. N. Likos, ACS Macro Lett., 2018, 7, 447.
19. A. Martin-Gomez, T. Eisenstecken, G. Gompper and R. G. Winkler, Phys. Rev. E, 2020, 101, 052612.
20. N. Kikuchi, A. Gent and J. M. Yeomans, Eur. Phys. J. E: Soft Matter Biol. Phys., 2002, 9, 63.
21. K. Mussawisade, M. Ripoll, R. G. Winkler and G. Gompper, J. Chem. Phys., 2005, 123, 144905.
22. R. Chelakkot, R. G. Winkler and G. Gompper, Phys. Rev. Lett., 2012, 109, 178101.
23. J. F. Ryder and J. M. Yeomans, J. Chem. Phys., 2006, 125, 194906.
24. S. Frank and R. G. Winkler, EPL, 2008, 83, 38004.
25. C.-C. Huang, R. G. Winkler, G. Sutmann and G. Gompper, Macromolecules, 2010, 43, 10107.
26. M. Ripoll, R. G. Winkler and G. Gompper, Phys. Rev. Lett., 2006, 96, 188302.
27. S. P. Singh, D. A. Fedosov, A. Chatterji, R. G. Winkler and G. Gompper, J. Phys.: Condens. Matter, 2012, 24, 464103.
28. S. P. Singh, R. G. Winkler and G. Gompper, Phys. Rev. Lett., 2011, 107, 158301.
29. J. T. Padding and A. A. Louis, Phys. Rev. Lett., 2004, 93, 220601.
30. M. Ripoll, P. Holmqvist, R. G. Winkler, G. Gompper, J. K. G. Dhont and M. P. Lettinga, Phys. Rev. Lett., 2008, 101, 168302.
31. Z. Tan, M. Yang and M. Ripoll, Soft Matter, 2017, 13, 7283.
32. S. Bucciarelli, J. S. Myung, B. Farago, S. Das, G. A. Vliegenthart, O. Holderer, R. G. Winkler, P. Schurtenberger, G. Gompper and A. Stradner, Sci. Adv., 2016, 2, e1601432.
33. S. Das, J. Riest, R. G. Winkler, G. Gompper, J. K. Dhont and G. Nägele, Soft Matter, 2018, 14, 92.
34. H. Noguchi and G. Gompper, Phys. Rev. Lett., 2004, 93, 258102.
35. H. Noguchi and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 14159.
36. J. Hu, M. Yang, G. Gompper and R. G. Winkler, Soft Matter, 2015, 11, 7867.
37. M. Wagner and M. Ripoll, EPL, 2017, 119, 66007.
38. K. Qi, E. Westphal, G. Gompper and R. G. Winkler, Phys. Rev. Lett., 2020, 124, 068001.
39. Y.-G. Tao and R. Kapral, Soft Matter, 2010, 6, 756.
40. A. Zöttl and H. Stark, Phys. Rev. Lett., 2012, 108, 218104.
41. S. B. Babu and H. Stark, New J. Phys., 2012, 14, 085012.
42. J. Elgeti and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4470.
43. S. Y. Reigh, R. G. Winkler and G. Gompper, Soft Matter, 2012, 8, 4363.
44. M. Theers and R. G. Winkler, Soft Matter, 2014, 10, 5894.
45. M. Yang and M. Ripoll, Soft Matter, 2014, 10, 1006.
46. S. M. Mousavi, G. Gompper and R. G. Winkler, Soft Matter, 2020, 16, 4866.
47. M. Yang and M. Ripoll, Soft Matter, 2016, 12, 8564.
48. Z. Tan, M. Yang and M. Ripoll, Phys. Rev. Appl., 2019, 11, 054004.
49. Y. Hashimoto, Y. Chen and H. Ohashi, Comput. Phys. Commun., 2000, 129, 56.
50. Y. Inoue, Y. Chen and H. Ohashi, J. Comput. Phys., 2004, 201, 191.
51. T. Ihle, E. Tüzel and D. M. Kroll, EPL, 2006, 73, 664.
52. E. Tüzel, G. Pan, T. Ihle and D. M. Kroll, EPL, 2007, 80, 40010.
53. E. Tüzel, G. Pan and D. M. Kroll, J. Chem. Phys., 2010, 132, 174701.
54. C. Echeverria, K. Tucci, O. Alvarez-Llamoza, E. Orozco-Guillén, M. Morales and M. Cosenza, Front. Phys., 2017, 12, 1.
55. T. Hiller, M. Sanchez de La Lama and M. Brinkmann, J. Comput. Phys., 2016, 315, 554.
56. T. Eisenstecken, R. Hornung, R. G. Winkler and G. Gompper, EPL, 2018, 121, 24003.
57. S. Meßlinger, B. Schmidt, H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011901.
58. I. O. Götze, H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 046705.
59. R. Pesché, M. Kollmann and G. Nägele, J. Chem. Phys., 2001, 114, 8701.
60. G. Nägele, A. J. Banchio, M. Kollmann and R. Pesché, Mol. Phys., 2002, 100, 2921.
61. S. Panzuela, R. P. Peláez and R. Delgado-Buscalioni, Phys. Rev. E, 2017, 95, 012602.
62. S. Panzuela and R. Delgado-Buscalioni, Phys. Rev. Lett., 2018, 121, 048101.
63. B. U. Felderhof, J. Chem. Phys., 2005, 123, 184903.
64. K. Huang and I. Szlufarska, Nat. Commun., 2015, 6, 1.
65. D. Lopez and E. Lauga, Phys. Fluids, 2014, 26, 071902.
66. J. Hu, A. Wysocki, R. G. Winkler and G. Gompper, Sci. Rep., 2015, 5, 9586.
67. P. Janmey and P. Kinnunen, Trends Cell Biol., 2006, 16, 538.
68. R. Kapoor, T. A. Peyear, R. E. Koeppe and O. S. Andersen, J. Gen. Physiol., 2019, 151, 342.
69. J. Corrie, L. Dey, J. D. Therien and J. E. Baenziger, Nat. Chem. Biol., 2013, 9, 701.
70. L. O. Romero, A. E. Massey, A. D. Mata-Daboin, F. J. Sierra-Valdez, S. C. Chauhan, J. F. Cordero-Morales and V. Vásquez, Nat. Commun., 2019, 10, 1.
71. R. Caires, F. J. Sierra-Valdez, J. R. Millet, J. D. Herwig, E. Roan, V. Vásquez and J. F. Cordero-Morales, Cell Rep., 2017, 21, 246.
72. J. Oates and A. Watts, Curr. Opin. Struct. Biol., 2011, 21, 802.
73. M. Manna, M. Niemelä, J. Tynkkynen, M. Javanainen, W. Kulig, D. J. Müller, T. Rog and I. Vattulainen, eLife, 2016, 5, e18432.
74. A. Malevanets and R. Kapral, J. Chem. Phys., 2000, 112, 7260.
75. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 020201.
76. R. G. Winkler and C.-C. Huang, J. Chem. Phys., 2009, 130, 074907.
77. M. Belushkin, R. G. Winkler and G. Foffi, J. Phys. Chem. B, 2011, 115, 14263.
78. T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 066706.
79. C.-C. Huang, A. Varghese, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 013310.
80. E. Tüzel, T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 056702.
81. H. Noguchi, N. Kikuchi and G. Gompper, EPL, 2007, 78, 10005.
82. M. Theers and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 033309.
83. M. Ripoll, K. Mussawisade, R. G. Winkler and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 016701.
84. J. Padding, A. Wysocki, H. Löwen and A. Louis, J. Phys.: Condens. Matter, 2005, 17, S3393.
85. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101.
86. M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Soft Matter, 2016, 12, 7372.
87. M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590.
88. A. Lamura, G. Gompper, T. Ihle and D. M. Kroll, EPL, 2001, 56, 319.
89. M. Hecht, J. Harting, T. Ihle and H. J. Herrmann, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011408.
90. J. Dhont, An Introduction to Dynamics of Colloids, Elsevier Science, 1996.
91. E. Tüzel, T. Ihle and D. M. Kroll, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 056702.
92. M. Doi and S. Edwards, The Theory of Polymer Dynamics, Clarendon Press, 1988.
93. J. Bławzdziewicz, M. L. Ekiel-Jeżewska and E. Wajnryb, J. Chem. Phys., 2010, 133, 114702.
94. S. P. Singh, G. Gompper and R. G. Winkler, J. Chem. Phys., 2018, 148, 084901.
95. J. T. Padding and W. J. Briels, J. Chem. Phys., 2010, 132, 054511.
96. S. Shkulipa, W. den Otter and W. Briels, Biophys. J., 2005, 89, 823.
97. W. den Otter and S. Shkulipa, Biophys. J., 2007, 93, 423.
98. D. M. Engelman, Nature, 2005, 438, 578.
99. A. P. Minton, Curr. Opin. Struct. Biol., 2000, 10, 34.
100. J. S. Kim and A. Yethiraj, Biophys. J., 2009, 96, 1333.
101. Y. Wang, M. Sarkar, A. E. Smith, A. S. Krois and G. J. Pielak, J. Am. Chem. Soc., 2012, 134, 16614.
102. H. Matsuda, G. Putzel, V. Backman and I. Szleifer, Biophys. J., 2014, 106, 1801.
103. M. Hellmann, D. W. Heermann and M. Weiss, EPL, 2012, 97, 58004.
104. L. E. Sereshki, M. A. Lomholt and R. Metzler, EPL, 2012, 97, 20008.
105. B. Liebchen, P. Monderkamp, B. ten Hagen and H. Löwen, Phys. Rev. Lett., 2018, 120, 208002.
106. Jülich Supercomputing Centre, JLSRF, 2018, 4, A132.
107. H. S. Carslaw and J. C. Jaeger, Conduction of heat in solids, The Clarendon Press, NY, 1988.

---