Pressure-driven flow of oligomeric fluid in nano-channel with complex structure. A dissipative particle dynamics study

We develop a simulational methodology allowing for simulation of the pressure-driven flow in the pore with flat and polymer-modified walls. Our approach is based on dissipative particle dynamics and we combine earlier ideas of fluid-like walls and reverse flow. As a test case we consider the oligomer flow through the pore with flat walls and demonstrate good thermostatting qualities of the proposed method. We found the inhomogeneities in both oligomer shape and alignment across the pore leading to a non-parabolic velocity profiles. The method is subsequently applied to a nano-channel decorated with a polymer brush stripes arranged perpendicularly to the flow direction. At certain threshold value of a flow force we find a pillar-to-lamellar morphological transition, which leads to the brush enveloping the pore wall by a relatively smooth layer. At higher flow rates, the flow of oligomer has similar properties as in the case of flat walls, but for the narrower effective pore size. We observe stretching and aligning of the polymer molecules along the flow near the pore walls.


Introduction
Understanding the behavior of polymers attached to surfaces is of importance in many research areas including biophysics, polymer-induced effective interactions in colloidal suspensions, chromatographic separation, catalysis, and drug delivery [1]. Grafting polymer chains can significantly alter the properties of the surface and make it, for example, biocompatibile or responsive to external stimuli [2]. Due to the large field of potential applications, polymer brushes have been the subject of many theoretical studies. In the seminal works Alexander [3], de Gennes [4,5] have calculated the brush profile and explored analytically the impact of grafting density and molecular weight. Since then, the properties of tethered chains have been investigated by means of self-consistent field theory [6][7][8][9][10], polymer density functional theory [11][12][13][14][15][16], and computer simulation [17][18][19][20][21]. Many theoretical predictions have been confirmed by experiment [22][23][24][25].
Polymer brushes can be used to tailor static properties of surfaces, such as wettability, as well as dynamical, such as hydrodynamic boundary conditions and friction. Fluid flow in polymer grafted nanopores can be described via continuum hydrodynamic equations (e.g., the Brinkman equation [26]) with a priori assumed permeability related to the monomer density profile. The resulting velocity profile is sensitive to the assumed form of the monomer profile [27]. However, the continuum hydrodynamic description of a flow has not been firmly established on the nanoscale [28]. This is important in the context of micro-and nanofluidic devices [29]. Downsizing a channel to the nanoscale, increases the surface-tovolume ratio and introduces new physical phenomena not observed in the macroscale [30]. Covering the surface by a polymer brush may introduce a pronounced reduction of friction, which lowers the pressure difference required to maintain the flow through a nanochannel [31]. Flow in polymer brushes has been the subject of numerous simulational studies in recent years [32][33][34][35][36][37][38].
Recently, the equilibrium properties of binary mixture confined in a slit-like pore decorated with polymer brush stripes were studied by means of dissipative particle dynamics (DPD) [39,40]. It was found that, depending on the geometrical parameters characterizing the system (the size of the pore and the width of the stripes), several different structures (or morphologies) inside the pore can be formed. Such patterned brushes can be fabricated experimentally by means of electron beam litography [41]. In the present paper, we wish to study nonequilibrium properties of such system by considering the pressure driven oligomer flow inside a channel with either flat or brush-modified walls. In particular, we focus on three features such as: (i) the microstructure of a flow depending on the molecular mass of an oligomer and the magnitude of a bulk flow force; (ii) flow-induced morphology changes; and (iii) the effect of the patterned brush decoration of the walls on the properties of the flow. Our paper is arranged as follows: In section 2 we introduce a new simulation method which combines the ideas of fluid-like walls and reverse flow to minimise the near-wall artefacts and maintain constant temperature under flow condition. As a simple test, we apply the method to the case of oligomer flow through the pore with flat walls. In section 3, the analysis is extended to the case when the walls are modified by a polymer brush arranged in a form of stripes. Conclusions are provided in section 4.

Flow of oligomeric fluid through a channel with flat walls
Let us first consider the simulational approach employed in this study. We use the non-equilibrium extension of the DPD technique in a form discussed by Groot and Warren [42]. This is a mesoscopic method that operates at a level of coarse-grained beads, each representing either a fragment of a polymer chain or a collection of solvent particles. The force acting on i th bead due to its pairwise interaction with j th bead can be written as where F C i j , F D i j and F R i j denote the conservative, dissipative and random contribution, respectively. These have the following form [42] and ∆t is the time-step of the integrator. As already discussed in references [28,65], the effective range of friction between beads can be modified by adjusting the shape of the weight functions w D (r i j ) and w R (r i j ). We use the following general form for w R (r i j ): where the exponent β is adjusted, and the weight function w D (r i j ) is set equal to [w R (r i j )] 2 according to Español and Warren [43] arguments. Likewise, it is required that σ 2 = 2γ.
The oligomers and tethered polymer chains (if any) are represented as necklaces of beads bonded together via harmonic springs, the force acting on i th bead from the interaction with its bonded neighbour, j th bead, is

13609-2
where r i j = |r i j |, r i j = r i − r j is the vector connecting the centers of i -th and j -th beads,r i j = r i j /r i j and k is the spring constant. The same bonding force is used to tether the end polymer bead to the surface. The length, mass, time and energy (expressed via T * ≡ k B T ) units are all normally set equal to unity.
Let us now turn to the case where the fluid (or a mixture of fluids) is confined within a slit-like pore. In order to commence a simulation of the pressure-driven flow, it is required to provide a set of rules defining the behaviour of the fluid particles at walls, and a prescription for the construction of the walls. These rules should recover the well known cases of hydrodynamic flow such as the Poiseuille flow (i.e., a flow of a Newtonian fluid with no-slip boundary conditions and a parabolic velocity profile). On the other hand, for the flow of a polymeric fluid (i.e., a non-Newtonian fluid) the set of rules should lead to the slip boundary conditions. The simplest set of rules comprise elastic reflections off the wall [44,45]. Unfortunately, they give rise to a hydrodynamic slip for Newtonian fluids, as well as suffer from nearwall density artifacts at higher density. This can be traced back to the fact that the atoms repelled each other strongly but did not interact with the wall until they attempted to cross [46].
A number of more sophisticated set of rules have been suggested. One option is to form the crystalline walls of a few layers of frozen (or having large mass) particles [47][48][49][50][51][52][53]. The interaction between the bulk fluid particles and those of the wall creates the near-wall drag which leads to the formation of the Poiseuille flow. A drawback of this approach is the propagation of the crystalline order into the nearwall regions of bulk fluid. This effect is perfectly physical for the atomic molecular dynamics simulation, where the solid wall mimics a real crystalline structure. However, for the mesoscopic DPD simulations, each soft bead is assumed to represent a meso-scale portion of the material, on which scale the atomic crystalline structure is smeared-out.
This drawback can be avoided by using the structureless fluid-like walls [46]. The walls in this case are made of the fluid confined in the slabs adjacent to the pore boundary, and the elastic reflections are applied on both sides of the boundary. Therefore, bulk and wall fluid particles are immiscible. Still, the interaction between near-wall beads on the opposite sides of the boundary creates a near-wall drag ensuring no-slip boundary condition for Newtonian fluids.
Another important issue in flow simulation is to avoid the system overheating due to the presence of the body force. This problem was addressed in several studies, cf. for example references [54,55]. In the molecular dynamics simulation, the excessive energy is absorbed by an external thermostat, in either bulk or near-wall form [54]. In DPD simulations, the thermostat is "internal", provided by the balance between interparticle friction and random forces. For the case of a flow, some means for dissipation of additional energy related to the body force should be provided. One of the elegant ways to do this is the concept of a reverse flow [56][57][58]. In this approach, the simulation box contains two sub-flows driven oppositely. The total force applied to the system is equal to zero and a no-slip boundary is formed at the interface between two opposite flows of Newtonian fluids.
In our study, we combine both concepts by employing the fluid-like walls on both boundaries of a pore and initiating contraflows (reverse flows) within them. Separation between the main pore and contraflow-containing walls prevents intermixing between the beads from both regions. This is important both in the case when the flow of a mixture is considered, or in the case of polymer modified walls, where polymer chains are tethered to the boundary between the main pore and fluid-like wall. However, the existence of the reflective boundaries does not prevent a friction between the beads located on the opposite sides of the boundary, enabling the formation of the no-slip boundary condition for Newtonian fluids.
In this section, we consider the pressure-driven flow of oligomeric one-component fluid through the pore with flat walls. The oligomers of length L o = 1, 4, 10 and 20 beads are considered. The aim is twofold.
Firstly, we would like to test to what typical values of bulk force the approach outlined above can be stretched without violation of temperature conservation. Secondly, we aim to study the flow microstructure depending on molecular length of the flowing oligomer and the magnitude of a flow force. The geometry of the system is illustrated in figure 1. Here, X -axis runs from left to right, Z -axis -from bottom to top, Y -axis coincides with the viewing direction. The simulation box is of dimensions L x = 80, L y = 50 and L z = 26.667 with the periodic boundary conditions applied along X and Y axes, the pore size is d = 13.333, the size of the contraflow regions is c = d /2 = 6.667. The chains in contraflow regions are of the same length L o as in the main pore. Therefore, the total number of main and contraflow chains is the same. All beads are assumed to be of the same type, which is reflected in the fact that the parameter but, alternatively, the periodic boundary conditions can be used in Z direction, similarly to the original reverse flow setup [56][57][58].
Each i th bead within a pore is subjected to the flow force of certain amount f i applied along X -axis, where f i > 0, this is indicated by the right-hand side directed arrows in figure 1.
The beads in the contraflow regions are subjected to the force F FL i = − f ix , indicated as reversely directed arrows in the same figure. Several options are available for choosing the amount of f i . The simplest one would be to choose f i ≡ f , the same amount for each bead. However, such an algorithm could lead to less than optimal match of the micro-fluctuations of the applied pressure in real systems, since the polymer molecules tend to form coils with varying distribution of the density. Another, rather extreme option would be to apply the amount f L o to the middle bead only. The other beads feel this force indirectly and are delayed via the elastic spring forces. The latter approach might suffer from large fluctuations of bond lengths and slower relaxation of the intra-chain vibrations, due to the soft nature of the model.
In our view, a reasonable compromise can be achieved by applying a fixed amount of the force f L o to each oligomer, but biasing it towards the middle bead of the chain. Namely, assuming that the beads of an oligomer are numbered sequentially as l = 1, . . . , L o , then the amount of the force applied to the bead number l is found according to the Gaussian distribution: (2.8) Here,l = (L o + 1)/2 is the mid-index of the chain, and the breadth of the distribution is given The shape of the weight function w G (l ) is shown in figure 2 for the cases of L o = 4, 10 and 20. As a result, the total force applied to the oligomer of L o beads is equal to f L o , but it is biased towards the middle beads (illustrated by arrows of different length in figure 1). The acceleration of the fluid beads due to applying the bulk force affects the accuracy of the integrator, as far as the expression for the coordinates at the time instance t +∆t contains the term proportional to v(t )∆t , where v(t ) is the velocity of the particle at the time instance t . The only way to keep the same behaviour. The models describing such non-parabolic profiles exist (see, e.g., reference [58]) and involve an analogue for the viscosity and a number of additional parameters. We found, however, a numerical fitting to these forms impractical. The set of rules defining the behaviour of the particles at walls, as imposed in our simulation, leads here to the slip boundary conditions. The discontinuity of the velocity profile at the wall boundary is clearly visible in figure 4 (c) and is a characteristic feature of the flows of polymeric fluids. In figure 4 (d), we compare two velocity profiles of the flows obtained with applying equation (2.8)-(2.9), i.e., the Gaussian distribution of the bulk force, and a uniform distribution of the bulk force. We note that even for such an extremely large value of the bulk force, the profiles are practically identical. We expect that for very long polymers, the Gaussian distribution of the bulk force would prove  beneficiary and could lead to better stability of the integration of the equations of motion.
An alternative route is to concentrate on the details of the microstructure of the oligomer flow, because these must be responsible for its non-Newtonian behaviour. In particular, comparing to the case of a simple fluid, oligomers have additional conformational degrees of freedom which will affect their flow properties. Therefore, we build the profiles for the average shape anisotropy and the molecular orientation for the oligomers in a flow. The components of the gyration tensor are evaluated for each oligomer of length L o at a given time instance t . Here, α, β denote the Cartesian axes, r i ,α are the coordinates of i th monomer, and R α are the coordinates for the center of mass of the oligomer. In the equivalent ellipsoid representation, the eigenvalues λ 1 > λ 2 > λ 3 of this tensor provide the squared lengths of its semiaxes, whereas the respective eigenvectors u 1 , u 2 and u 3 -the orientation of these axes in space. The shape anisotropy of an individual oligomer can be defined as It is zero for a spherically symmetric body, where λ 1 = λ 2 = λ 3 > 0 and is equal to 1 for an infinitely long thin rod, where λ 1 > 0, λ 2 = λ 3 = 0. The average profile is built for the shape anisotropy in a steady state. It is obtained by first binning the system in Z -axis and averaging κ 2 for individual oligomers found in each bin. Then, time averaging within the steady state is performed. The orientation of each oligomer in space is defined by that for the longest axis of its equivalent ellipsoid. The latter is characterised by the eigenvector u 1 associated with the largest eigenvalue λ 1 . The level of alignment of the oligomer along the flow axis X can be characterised by the order parameter: S x = P 2 (u 1 ·x), (2.12) wherex is defined in equation (2.7) and P 2 (x) is the second Legendre polynomial. The alignment profile is built then in a steady state by averaging S x in each bin and then performing time averaging. It is obvious that both κ 2 and S x can be defined for the case L o > 1 only. f . It is non-zero in the middle of a channel for all f being considered. Therefore, at least for longer oligomers, L o > 4, there is a variation of the oligomer shape and alignment across the channel: the molecules are found to be much more elongated and aligned near the edges as compared to the middle part. The flow-induced deformation of the polymer molecules renders their shape to be more similar to liquid crystals. The effect is detected for longer oligomers L o > 4 and stronger flows, where the effective length-to-breadth ratio of oligomer exceeds a certain threshold. Similar effect is well known for the systems of anisotropic hard bodies, where the orientationally ordered phases are also found above certain threshold length-to-breadth ratio [59][60][61][62]. Using this liquid crystal analogy, we recall the results obtained by Mazza et al. [63,64] reporting the high self-diffusivity of the Gay-Berne-Kihara fluid along the director in the "supernematic" phase. Following these findings, one expects an essential reduction of the friction between the aligned oligomers near the channel edges, as compared to that in the central part. Larger friction between oligomers in the middle of a channel is seen as the reason for the suppression of the velocity profile here and, as a result, its non-parabolic, bell-like shape [cf. figure 4 (b)], and the appearance of the slip boundary condition [cf. figure 4 (c)].

Flow of oligomeric fluid through a channel with polymer modified walls
We turn now to the case when the pore walls are modified by polymer brushes arranged in the form of stripes (see, figure 6). Each chain of a brush is of length L = 20 beads of type A, the pore interior is filled with the oligomer fluid of the length L o beads of type B , the contraflow regions contain an oligomer fluid of the length L o beads of type A. The difference between the bead types is in the value of the repulsion amplitude a in equation (2.2) being set to a A A = a B B = 25 and a AB = 40 for the interaction of similar and dissimilar beads, respectively. Therefore, the oligomer acts as a bad solvent for the brush. The good solvent case, a A A = a B B = a AB = 25, is briefly discussed in the end of this section. The equilibrium properties of the setup depicted in figure 6 for the case of L o = 1 and no contraflow regions are studied in detail in reference [39,40]. Equilibrium morphology was found to depend on the parameters d and w, and is formed as a result of an interplay between the enthalpy and the entropy of the system. In particular, at small w L and any d , the adjacent brush stripes belonging to the same wall merge and form a homogeneous "coat" on the wall resulting in the lamellar morphology. In this case,  Colours follow these in figure 6, contraflow regions not shown. Firstly, at f 0.1, the profiles exhibit two "shoulders" near each internal wall which are characterized by zero values for v x . These are, obviously, the regions occupied by the polymer brush which "envelopes" the internal walls (see, figure 7). The flow is completely suppressed within these layers, rendering the walls thicker and reducing the pore size accessible to the flow. As a consequence, the maxima for v x decrease compared to the case of flat walls. Secondly, the shape of the central part of rrrrrrrrrrrrrrrrrrrr r r r r r r r r r r rr r r r r r r r r r rrrrrrrrrrrrrrrrrrrrrr r r r r r r r r r r r r r r rr r r r rrr This interpretation brings up the possibility to treat a fluid flow within a stationary lamellar morphology similarly to the case of the pore with flat boundaries, discussed in section 2, except for the smaller effective pore size d eff [65]. To evaluate the latter, one can use the expression for an average brush thick- where ρ p (z) is the density profile of the beads that belong to the tethered chains andz is the distance from the nearest pore boundary along the Z -axis. In this case, one obtains d eff = d − 2b h . Alternatively, the effective pore size can be estimated as the distance between the intersection points z 1 and z 2 for ρ p = ρ p /ρ and ρ s = ρ s /ρ, the reduced density profiles for the polymer and the flowing oligomer beads, respectively. This is illustrated in figure 10 figure 10). While d eff can be also calculated at smaller values of f , these results would carry no physical significance due to the pillar morphology. Second, the difference between the values for d eff estimated by means of two alternative methods at the same L o and at the same f , does not exceed 4%. Therefore, either of the estimates for d eff can be used. Third, there is a trend for an increase of d eff with the growth of oligomer length L o , although it is rather modest. For example, for the case of f = 0.4 the value of d eff for oligomer length of L o = 20 is only 10% higher than its counterparts for L o = 1 and 4, and this increase is of the order of the error in the estimates of d eff mentioned above. A slight increase of the brush height with an increasing flow can be attributed to the fact that there is some residual flow of oligomers inside the brush. As the flow increases, the flow-induced elongation of the oligomers leads to an increase of their effective size and this will lead to a slight increase of the brush height.
As was discussed in section 2, flat walls of the setup depicted in figure 1 act as effective "stretchers" and "aligners" for the adjacent oligomer molecules, which results in characteristic profiles for κ 2 and S x shown in figure 5. It is, therefore, of interest to see whether or not the effective walls formed by a flattened polymer brush, as pictured in figure 7 (b) and (c), have a similar impact on the adjacent oligomer molecules. We examine the aligning capabilities of such flattened brushes more in detail, considering both cases of bad and good oligomer solvent. For the former case, the repulsion parameter a in equa-

13609-11
in figure 11 are extremely close to their respective counterparts for the case of flat walls shown in figure 1 in both shape and absolute values, save for being "squeezed" into the smaller pore size d eff . This indicates that the existing roughness of the flattened polymer wall does not reduce its impact on the adjacent oligomer molecules.
For the good solvent case, the brush and oligomer mix well if no flow force is applied, but a flowdriven lamellarization of the system takes place at about f = 0.2, similarly to the bad solvent case. The profiles for κ 2 and S x are shown in figure 12 for the oligomer lengths of L o = 4 and L o = 20. The estimated effective pore size d eff is of the same order but a fraction smaller than that for the bad solvent case. This is indicated in figure 12 by vertical dashed lines. One should remark that despite the strong alignment of the polymer brush and good mixing between the brush and the oligomers, the latter are found less elongated and less aligned along a flow compared to the bad solvent case shown in figure 11. The respective curves for κ 2 are lower by about 0.1 compared to their counterparts for the bad solvent case, whereas these for S x are about 0.2 lower. One should attribute this to the fact that the flattening of the brush in the case of a good solvent requires a higher flow force compared to the case of a bad solvent. In the latter case, lamellarization is also aided by the microphase separation between the brush and the oligomer fluid. Despite these quantitative differences, the qualitative picture emerging for both cases of the bad and good solvent is essentially the same. Namely, at a certain value of the flow force f , the stationary lamellar phase is formed with the flowing oligomer occupying the center part of the pore. The oligomer is found essentially elongated and aligned along the flow near the walls of this channel and much less in the center of the pore. This effect is detected for the oligomer lengths L o > 4 and, due to its impact on the distribution of the local friction across the pore, affects the behaviour of the fluid turning it into a non-Newtonian one.

Conclusions
In this paper, we developed the simulation approach, which allows one to simulate the pressuredriven flow in the pore with flat and polymer-modified walls. It combines the earlier ideas of fluid-like walls and reverse flow. The former enables to avoid highly structured solid walls that usually lead to near-wall artefacts. The latter introduces friction between oppositely flowing streams which makes it possible to conserve the total momenta and keep the temperature constant. Our system geometry contains the central main pore "enveloped" by two fluid-like pores on each side. The flow force is introduced in the main pore and oppositely directed contra-flow of the same magnitude -in both fluid-like walls. Simulation of the oligomer flow through the pore with flat walls is used as a check for the credibility of the method and reproduction of the hydrodynamic boundary conditions. Good thermostatting of the system is achieved when the flow force magnitude does not exceed a certain threshold. For the case of the oligomer length L o > 4, we found the molecules adjacent to the central pore boundaries essentially stretched and aligned along the flow, whereas their shape is more spherical and less aligned in the middle of a pore. This provides the basis for variation of the local friction across the pore and, as a consequence, the non-parabolicity of the velocity profile for the oligomer fluid and the slip boundary condition.
The case of polymer-modified walls is also considered when the polymer brush has the form of stripes arranged perpendicularly to the flow direction. In this case, at a certain threshold value of the flow force, one observes the pillar-to-lamellar transition induced by the flow which leads to the brush enveloping the pore wall with a relatively smooth layer. At higher flow rates, the flow of oligomer is similar to the case of flat walls, although for the narrower effective pore size. The latter is estimated both from the intersection of density profiles for the brush and the flowing oligomer and from the integral equation for the average brush thickness. The effect of local stretching and alignment of oligomers near the walls of the effective pore is detected the same as for the case of flat walls.
The method can be extended to more complex systems, namely: the flow of mixtures and their flowinduced separation; the flow of amphiphilic molecules; the flow of complex macromolecules or their solutions. Combined with the fine-tunable structure of the brush, this opens up a possibility to study various problems of transport of oligo-and macromolecules through a complex structured environment.