Exploring the Molecular Distributions in Dilute Polymer Solutions Using a Multi-Scale Numerical Solver

Simulating the rheological behaviors of polymer solutions is intrinsically a multi-scale problem. To study the macroscopic and microscopic characteristics in the fluid flow of dilute polymer solutions, we designed a multi-scale solver, which couples the Brownian Configuration Fields with the macroscopic hydrodynamic governing equations. Numerical simulation results using the multi-scale solver exhibited good accordance with the macroscopic only approach. Through a scalar field D we also quantitatively studied the flow behaviours in 2D planar channels, and analyzed the correlation between the molecular distribution and the macroscopic fluid flow in polymer solutions. Our results verified the correctness of the solver, which could provide valuable guidance for multi-scale simulations of complex fluids based on OpenFOAM.


Introduction
Polymer solutions are frequently used in the production of fibers, films, glues, lacquers, paints, and other items made of polymer materials. It is of great scientific significance and application values to make an in-depth study of the dynamics and the rheological characteristics of the polymer solutions through numerical simulation.
Simulating the rheological behaviors of polymer solutions is intrinsically a multi-scale problem. Multi-scale simulation is a rapidly growing multidisciplinary field, as stated in the review article provided in [1].
Numerical approaches for simulating polymer solutions fall into three categories: the macroscopic, the microscopic and the micro-macro multi-scale methods. In a macroscopic only simulation, the constitutive equation (CE) that relates the viscoelastic stress to the deformation history can be derived from the continuum mechanics. One then solves the constitutive model together with the conservation laws of mass and momentum to predict velocity and stress fields in complex flows. An extensive description of constitutive models is given by Owens and Phillips [2]. The Oldroyd-B model originated from a two-bead dumbbell model with a linear Hookean spring force is employed in our simulations to compute the stress tensor. Dumbbell models ignore the interactions among polymer molecules and they are usually employed to simulate the dilute polymer solutions. Multi-bead chains have been studied as well [3,4]. There are also other constitutive models such as the FENE-P (the closure approximation for finitely extensible nonlinear elastic dumbbell model proposed by Peterlin) model [5], the FENE-L (the closure approximation for finitely extensible nonlinear elastic dumbbell model proposed by Lielens) model and the FENE-LS (the simplification of closure approximation molecular distributions and the macroscopic fluid flow. In this paper, we aim to demonstrate the possibility of exploring the microscopic characteristics using micro-macro simulations.
In the remaining of this paper, we first present the governing equations on the macro-and micro-scale in Section 2. Then, Section 3 discusses the spatial and temporal numerical discretization schemes of the Navier-Stokes, Oldroyd-B and stochastic differential equations. In Section 4, we give the implementation of a multi-scale numerical solver, BCFsolver. The mesh geometry and initial conditions used in the simulations, along with the numerical results, are presented in Section 5. At last, we validate our findings and discuss potential extensions.

Governing Equations
The flow of incompressible and isothermal dilute polymer solutions could be governed by the macroscopic and the microscopic description, respectively. In this section the momentum balance equation, the continuity equation and the stress tensor equation are first presented on the macro-scale, then a Fokker-Planck and a stochastic differential equation are employed respectively to replace the macroscopic stress tensor equation on the micro-scale.

Macroscopic Equations
We consider the isothermal and incompressible fluid of polymer solutions with density ρ. From the macroscopic perspective, the motion of a polymer solution fluid can be governed by the momentum balance and the continuity equation: ρ( ∂u u u ∂t + u u u · ∇ ∇ ∇u u u) = −∇ ∇ ∇p + η s ∆u u u + ∇ ∇ ∇ · τ τ τ p (1) ∇ ∇ ∇ · u u u = 0 where η s is the solvent viscosity and τ τ τ p is the time-dependent viscoelastic stress contributed from the polymer dynamics. The u u u and p are the velocity of the fluid and the pressure field, respectively. The polymer contribution to the stress can be calculated through a constitutive equation. On the macro-scale, the Oldroyd-B model is considered. The Oldroyd-B constitutive equation [31] can be expressed as where the symmetric deformation tensor D D D = 1 2 (∇ ∇ ∇u u u + (∇ ∇ ∇u u u) T ), λ is the relaxation time of the dumbbell system and η p is the polymeric viscosity. The upper convected derivative of the polymer stress tensor ∇ τ τ τ p [31], is expressed as The set of coupled Equations (1)-(3), supplemented with suitable initial and boundary conditions is the so-called macroscopic formulation of viscoelastic flows.
By scaling the equations with the characteristic units L c (characteristic length in macroscopic flow), U U U c (characteristic fluid velocity), ρ c (fluid density, scaling pressure term with 1 (ρU U U 2 c )) and normalizing the polymeric stress tensor with L c (U U U c (η s + η p )) , Equations (1)-(3) can be rewritten in a dimensionless form as: in which the dimensionless parameters De (Deborah number), Re (Reynolds number) and β (viscosity ratio) are defined as For simplicity, the same notation for the knowns are used in the dimensionless formulation as in (1)-(3).

Microscopic Equations
From a microscopic viewpoint, the Oldroyd-B equation originates from a Hookean dumbbell model. In this model, the molecules in polymer solution are represented by massive elastic dumbbells. All of them consist of two Brownian beads connected with a linear spring. The spring denotes inter-molecular forces between the two beads. The length and orientation of the spring describe a dumbbell's configuration indicated by a vector q q q, which has dimension of length. The spring force reads as F F F( q q q) = C q q q (C is the spring constant). In the remainder of this paper, a dumbbell's configuration vector q q q and the spring force F F F( q q q) are scaled as dimensionless quantity q q q and F F F(q q q), respectively. Refer to [28] for more detailed information about how the models are nondimensionalized.

Fokker-Planck Approach
The probability of finding a dumbbell with a configuration q q q at (r r r, t) is indicated by ψ(r r r, q q q, t). ψ is a probability density function (pdf) and fulfills ψ(r r r, q q q, t) ≥ 0, ψ(r r r, q q q, t)dq q q = 1 in the geometrical domain and configuration space at any time of the simulation. By applying Newton's second law to the forces acting on a dumbbell system, we can get the dimensionless Fokker-Planck equation [31].
The Equation (9) describes the evolution of ψ under the dumbbell's spring force F F F(q q q). In Equation (9) we dropped the diffusion term ∆ r ψ in physical space because the diffusion coefficient scales quadratically in the micro-macro length scale ratio L c /l c and is usually only of the order 10 −8 [32].
In this paper, we employed two different spring forces [31] to characterize elastic forces: F F F(q q q) = q q q, (Hooke) (10) F F F(q q q) = q q q 1− q q q 2 /b , q q q 2 ≤ b (FENE) (11) where b = q q q max 2 l c is dimensionless unit for the dumbbell's maximum extension q q q max compared to a characteristic micro length l c . Note that for the FENE spring force in Equation (11), the dumbbell's extension is restricted. The pdf ψ in Equation (9) represents the polymeric configuration of the micro-system. Now Kramers' equation [31] reads as: where · = ·ψ(r r r, q q q, t)dq q q denotes the expectation in configuration space and I I I is a unity matrix [28]. The prefactor α b,d specifies a spring dependent constant and is defined as , for FENE dumbbells (d is the dimension of configuration space).
The conservation of angular momentum leads to a symmetric τ τ τ p .

Stochastic Brownian Configuration Field Approach
The BCF method replaces the Fokker-Planck equation with a corresponding stochastic differential equation [31].
where Q Q Q t is a stochastic process and represents the configuration vector q q q. The 3-dimensional Gaussian Wiener process W W W t used to model Brownian forces is described by its first and second moments W W W t = 0 and W W W t W W W t = min(t, t )I I I. To approximate the actual probability density function, we solve Equation (14) for a number of stochastic realizations, namely Brownian configuration fields, represented by Q Q Q (i) t , i = 1, 2, · · · , N BCF . Using a Monte Carlo integration, the first moment Q Q Q t ⊗ F F F(Q Q Q t ) in Kramers' Equation (12) can be approximately calculated by As a result, the polymer contribution to the extra-stress τ τ τ p can be given by From Equation (16), we can conclude that the accuracy of τ τ τ p and the computing time critically depend on the choice of N BCF .

Numerical Methods
The momentum, continuity and constitutive equations are discretized through the finite volume method, which computes each term of the governing equations via volume integral over a control volume, which assumes a local conservation of physical laws. The 2nd-order Gauss MINMOD and Gauss linear scheme are employed to discretize the spatial terms. Temporal terms are discretised through a simple Euler scheme. The resulting equations would finally reduce to linear systems. As a result, we can get the solutions of these equations at each time step through iterative solvers predefined in OpenFOAM. The conjugate (PCG) and biconjugate gradient (PBiCG) methods are typical linear solvers in the toolbox. More details can be derived from the OpenFOAM Manual [33].
The discrete elastic-viscous split stress (DEVSS) numerical strategy of Guénette and Fortin [34] was employed to solve the governing equations presented in Section 2. And the numerical method modifies the PISO (pressure-implicit split-operator) algorithm by explicitly introducing polymer stress unknowns into the momentum equation as source terms. The momentum equation is rewritten as Discretized through finite volume method, this equation can be written as a linear system: where the polymer stress contribution is included as a source term in the symbol H H H n . The convective and viscous terms are implicitly discretized into the matrix A A A of the linear system. The square brackets [·] indicates numerical approximation of the corresponding variable; n and n + 1 signify the past and present times of the variable, respectively. Multiplying Equation (18) where U U U (n+1) = A A A −1 H H H n . A Poisson equation is derived as Equation (20) for solving the pressure in the pressure-correction step by taking the divergence of Equation (19) and applying the continuity Equations (19) and (20) are the key steps for the PISO algorithm. In the BCFsolver, the consitutive equation part is replaced with a microscopic BCF method to solve the viscoelastic stress tensor. A semi-implicit Eular method is employed to solve the equations of the Brownian configuration fields. The Equation (14) can be rewritten as Equation (21) is used to get the configuration field Q Q Q i (r r r, t) at time t n+1 for the FENE spring force. In case of the Hookean spring force, replace the 2De(1− Q Q Q (i) n+1 (r r r) 2 /b) with 2De. N N N(0, 1) denotes a triple of independent Gaussian random variables with zero mean and variance one.
To compute the polymer stress tensor, Equation (16) can be rewritten using semi-implicit Eular method as Coupling the macroscopic and microscopic parts, the iterative algorithm for solving the multi-scale model is summarised as Algorithm 1.

Algorithm 1:
The iterative algorithm to solve the multi-scale model. Data: mesh data, initial conditions Result: u u u, p, τ τ τ p 1 read the mesh data and the initial conditions; 2 initialization; 3 while t n+1 not reach the end of the simulation time do 4 for i = 1 to N BCF do 5 solve Equation (21) for to get the configuration field Q Q Q i (r r r, t) at time t n+1 ; 6 end 7 solve Equation (22) to get the polymer stress tensor τ τ τ n+1 p (r r r); 8 solve the discretised momentum Equation (17) to obtain the estimated components U U U (n+1) * ; 9 Using the estimated velocity U U U (n+1) * , solve the pressure-correction Equation (20) to obtain the pressure field p (n+1) * ;

BCFsolver Implementation
To solve the micro-macro governing equations, we implement a multi-scale solver named BCFsolver according to the numerical Algorithm 1 mentioned above. In this section, we specify the implementation process of the BCFsolver.

Overall Structure
OpenFOAM is a free, open source platform that not only provides us with a rich set of application solvers but also allows users to write custom solvers for specific applications based on their original solvers. It is often used to study the fluid dynamics of polymer solutions [35,36]. OpenFOAM is written in C++ and provides user interfaces that can describe partial differential equations in a natural language-like way, which makes it easy to extend the theoretical models of existing solvers.
The system organization of the BCFsolver based on the OpenFOAM is described in Figure 1. The solver module designed by us based on the OpenFOAM application framework is mainly the part deepened in gray. According to the solution process, the entire numerical solution process includes three parts: pre-processing, numerical solution and post-processing. The pre-processing part performs mesh generation and mesh decomposition. Mesh generation breaks down a continuous geometric space into tiny grids, and the "blockMesh" tool generates corresponding grids according to the input configuration files. For complex geometries, OpenFOAM can import a mesh generated by third-party software tools. Mesh decomposition will decompose the grid into multiple parts for parallel computing.
The post-processing part mainly includes the data merging and the visualization of the results. After parallel computation, the results are distributed in the corresponding folders of each processor. Data merging consolidates the corresponding calculation results to analyse using the "reconstructPar" tool. OpenFOAM integrates an open source visualization software, "paraView" (5.0.1, Kitware lnc., Clifton Park, NY, USA), which enables basic presentation and analysis of results. At the same time, OpenFOAM can also output a variety of different data formats for third-party software to visualize.
The numerical solution includes four layers: the domain language interface, the discretization method, the linear solvers and the parallel communication interfaces. The parallel communication interfaces are developed based on the MPI library and provides communication facilities in the form of APIs (Application Programming Interface). A variety of parallel solvers are used to solve linear systems, and the iterative solution of a linear system is achieved by calling the underlying parallel communication interface. By using the FVM (Finite Volume Method) method, the differential equations are discretized into the form of a linear system. Domain language interfaces are the key to writing complex fluid solvers and provide the basic interfaces for solving differential equations in theoretical models.

Domain Programming Interface
In this section, we specify the domain programming interfaces that are closely related to the solver implementation. It can be concluded that the core of the numerical solution is to discretize the differential equation according to the approaches described in Section 3. After the partial differential equation is discretized using the finite volume method, it can be transformed into a linear equation.
where A A A is a coefficient matrix and constant vector b b b is usually called source vector. Column vector x x x consists of corresponding values on different spatial grids. The theoretical model consists of a series of partial differential equations, each of which includes multiple terms, such as time derivatives, convection terms, Laplace terms, and so on. These items can be separately discretized by the interface functions provided by the framework. The finite volume method first integrates each term on mesh volumes, which are then linearized in a certain discrete format as part of the linear Equation (23). The discrete format can be specified either by writing a program or by reading the input file at run time. Discrete programming equations in OpenFOAM are defined as static member functions for both classes, finite Volume Method (fvm) and finite Volume Calculus (fvc). The member function of fvm implicitly discretizes the difference equations into the form of a coefficient matrix, which is then added to the left side of the linear system. The member function of fvc returns the corresponding field vector through explicit discretization, eventually joining the right source term.
Major discrete interfaces provided in OpenFOAM are shown in Table 1. The details can be found in the OpenFOAM manual [33]. The function name column contains description of the function parameters, where phi can represent various types of body fields, including scalar, vector, and tensor types; rho is a scalar field; psi is a surface field that can be derived from interpolating the value stored in the volume center into the surface through the finite volume method; chi can represent a surface field or a body field. Note that highest dimension of the data type is a 9-dimensional tensor field, so the gradient term can not operate on the tensor field.
The differential operators corresponding to the governing equations in the solver are included in Table 1, so we can implement the discretization of theoretical models based on the above programming interface.

Implementation
The numerical solver, BCFsolver, consists of two parts. The core part is the solver's main program, as displayed in the left side of Figure 2. Before starting the iteration of the numerical solution algorithm, the main program first needs to read the input parameters through the "argList" module. If the current program is executed in parallel, "parRun" would be called to start the parallel process. "parRun" calls initialization interface of the underlying parallel library to start MPI and complete the initialization of the process and the communication domain. "createTime" initializes the simulation program's time control module, which controls the total simulation time, the time step size and so on. The "createMesh" module constructs the grid data structure based on the grid information entered. All the physical variables in the program are constructed based on the grid data structure. Therefore, the "createMesh" module must be initialized before "createFields". "createFields" initializes the various field variables needed in the simulation based on the grid data, including pressure field p, velocity field u u u, viscoelastic stress field τ τ τ p . The numerical solution algorithm module implements the iterative solution process described in Algorithm 1. The main program of the BCFsolver also needs dynamic running configuration information for execution. These running configurations include three parts, which are respectively stored in the "0", "system" and "constant" file directories. The "0" directory defines the value of each field variable at the initial time. The "system" directory defines the key configuration information for solving the system, which includes the four parts shown in Figure 2. The "constant" folder is configured with mesh data and model parameters.
The controlDict file contains process control variables of the simulation program. The start time of the simulation is generally set to 0. You can also continue from the last simulation saved interrupt. The end time determines the total physical length of the simulation process. The deltaT determines the time step of discrete simulation. To ensure the accuracy of simulation, the time step can not be set too large. It is generally considered that the largest time step should ensure that the fluid advancement distance should not be greater than the width of one grid in a single time step. Therefore, Courant number is defined as Co = v v v m ∆t ∆l to describe the ratio of the distance of fluid movement to the width of the grid in a time step, where v v v m is the maximum fluid velocity, ∆t is the size of the time step, and ∆l is the mesh width. The smaller the Courant number, the higher the time accuracy of the simulation, leading to heavy computing burden. The maximum Courant number is significantly less than one in general. The Courant number is an important reference for setting the time step. Therefore, the controlDict provides the option of dynamically adjusting the time step according to the value of Co. We can set a maximum Courant number "maxCo" and a maximum time step "maxDeltaT", and the size of the time step is dynamically adjusted depending on the fluid velocity. The simulation results need to be written to a file for post-processing. The writeInterval parameter defines the interval at which data is written.
The decomposeParDict file defines a grid-based parallel division. Different partitioning methods can be used to decompose the grid data into individual processors, depending on the complexity of the geometry of the simulation problem.
The fvSolution and fvScheme respectively define the linear equations solver used in the algorithm and numerical discrete format. The most basic linear equations solving algorithm is the conjugate gradient method, and the numerical solver implemented in this paper mainly uses the algorithm. The specific configuration is shown in Table 2.

Field variables Solver Preprocessor Error limit
The essence of a FVM calculation process is to satisfy the conservation law by integrating each differential term over the control volume in the equation. The numerical discretization of equation items determines the process of transforming field data on discrete grids into a system of linear equations, which has an important influence on the accuracy and numerical stability of the solution process. Discrete formats used in the solver are displayed in Table 3. Table 3. Numerical discrete format used in the solver.

Equation terms
Discrete format Accuracy Constant directory contains two parts of information. The first part is the definition of the model constants. All parameters can be set in the transportProperties, such as relaxation time, viscosity, fluid density and so on. The second part is the grid data definition.

Results and Discussion
We ran multi-scale simulations with different configurations on a cluster located in the State Key Laboratory of High Performance Computing of National University of Defense Technology. In this cluster, each computing node contains 12 Intel Xeon E5-2620 2.10 GHz CPU cores and a total main memory of 16 GB. Each calculation presented in the paper costs around 2∼7 days on 8 CPU cores.

Problem Specification
For an unsteady Poiseuille flow in a planar channel, we assume no-slip boundary conditions at the channel wall. We employ a channel of length 1.0 in the x-direction, height 0.1 in the y-direction and of width 0.01 in the z-direction. We prescribe the velocity component u x at the inlet and outlet boundary. Figure 3 shows the two-dimensional geometry and the initial boundary conditions. The parameters of the micro-macro model need to be set in the configuration file, "transportProperties", located in the "constant" folder. Some of the parameters remain constant throughout all of the simulations and are displayed in Table 4: In the simulations, we test the designed micro-macro numerical solver on different meshes. Table 5 lists the different number of grid cells and the corresponding mesh width used in our simulations.

Simulation Results
Since the stochastic approach is much more demanding in terms of memory and computing resources, we first perform the simulations on the mesh level l = 2. The fluid density ρ = 0.1 at first.    From Figures 4-6, we can draw a conclusion that the snapshots of u u u and p fields are almost the same for the three models. However, for the stress tensor field τ τ τ p , the snapshot of τ τ τ p for Oldroyd-B model is different from that for the other two models. For the Hookean dumbbell and FENE model, it is obvious that the τ τ τ p is not so symmetric as that for the Oldroyd-B model. The Brownian configuration field method introduces stochastic noise and the τ τ τ p is derived from the average of local ensemble of sample particles through mathematical statistics. It might lead to additional error to the solution. As a result, the noise-free macroscopic model yields better results than the multi-scale approach.
In Figure 7, we examine the velocity field u u u at position P 1 = (0.498, 0.01) over time. The curve of velocity component u x is plotted with zoomed pictures. And the zoomed region is indicated with red rectangles in the full picture. Compared to the macroscopic result, the curves of multi-scale methods for Hooke and FENE model show minor changes caused by stochastic noise. But they all exhibit overshoot at the beginning of the simulation. We cannot expect to yield better results for the multi-scale model than for the corresponding macroscopic model because of the introduced stochastic noise. However, it is possible that the simulation with higher N BCF performs better with respect to the macroscopic Oldroyd-B solution.

Hooke Model
Since all unknowns solely depend on y when the unknown fields evolve into steady state, we display the unknowns in a steady state of the macroscopic and Hookean dumbbell models along the y-axis. To investigate the effects of different situations, we change only one parameter every time we test the effect of corresponding one. For example, we change the number of stochastic realizations N BCF and the other parameters remain the same with those set in Figure 5. Analogously, we adjust Re (Reybolds number) through changing the density of the fluid. We also change the mesh levels according to Table 5 and carry on the simulation for the Hookean dumbbell model on different mesh levels. Figure 8 shows the evolution of stress component τ xx along y-axis for different parameters of the Hookean dumbbell model.
In Figure 8a, the minima of these curves increase as the number of stochastic realizations N BCF becomes larger. The minima of curve at N BCF = 2000 is close to 0 and the symmetry improved a lot. It can be concluded that the errors of simulation results become smaller with the increase of N BCF . Figure 8b shows that the curves of stress tensor component τ xx have no difference among each other. We can draw a conclusion that the value of Re number imposes little effect on the simulation of Hooke model. In Figure 8c, these curves are close to each other for different mesh levels although there are some errors. The results of stress component τ xx for Hooke model show systematic deviations compared to the equivalent Oldroyd-B constitutive model. Hooke model introduces stochastic noises. However, we employ low-order temporal discretization schemes which may impose less control on the noises than the higher ones. We will explore the reason in further study.  From Figure 9, we can see that the velocity component u x of Hooke model achieved good consistency with the Oldroyd-B model. It can be concluded that the Hooke model can be employed in the simulation for the flow behaviour in 2D planar channel.

FENE Model
For the FENE model, we get the simulation results about stress tensor component τ xx and velocity component u x for different conditions along the y-axis. We first change the Re number of the FENE model by adjusting the density of the polymer solutions. At the meantime, the stochastic realizations are set N BCF = 800 per grid cell and the simulation is implemented on the mesh level l = 2. For the FENE model, the dumbbell's extension parameter b is restricted and we set the extension parameter b = 100. The results are presented in Figure 10. We change the mesh levels according to Table 5. We set ρ = 100, N BCF = 800, b = 100 again for the simulation. The simulation results of stress tensor component τ xx and velocity component u x on different mesh levels are presented and compared in Figure 11. For the stress tensor component τ xx , the minima of these curves increase as the value of N BCF becomes larger. Compared with N BCF = 1000, the symmetry of the curve at N BCF = 2000 has improved. The curves of u x component coincide with each other. From Figure 12, we can draw a conclusion that the accuracy of the solutions to the FENE model improved as N BCF increases.
We also investigate the effect of the dumbbell's extension parameter b. The other parameters are set as follows: ρ = 0.1, N BCF = 800. Again, the simulations run on the mesh level l = 2. Figure 13 displays the results. The results of the FENE model with the dumbbell's extension parameter b = 100 is very close to the Hooke model for the present conditions. So we test and make a comparison. The simulations run on the mesh level l = 2 with ρ = 0.1 and N BCF = 800 for both models. The parameters remaining constant throughout the simulations are displayed in Table 4.
As Figure 14 shows, both the stress component τ xx and the velocity component u x are very close to each other. It can be concluded that the FENE model with the dumbbell's extension b = 100 displays almost the same as the Hooke model in the simulation of the flow behaviour in 2D planar channel.

Molecular Distribution
To investigate the molecular distribution and stretch of 2D viscoelastic Poiseuille flow, We designed a new scalar Field D to describe the length of the extended molecules. We called the scalar Field D as "stretch length". The stretch length D is defined as: where Q Q Q t(x) , Q Q Q t(y) , Q Q Q t(z) represent the component of Q Q Q t in the x, y, z directions, respectively. We implement the simulation with the parameters ρ = 0.1, N BCF = 800, b = 100 on mesh level l = 2 for FENE model. The stress tensor component τ xx and corresponding D component along the y-axis, as well as the distribution of the molecules and corresponding probability density function (pdf), are presented in Figure 15. P 3 is located in the center of the planar channel; P 1 , P 2 and P 4 , P 5 are vertically symmetrically located on both sides of P 3 . The curve of τ xx component shows that the value of τ xx at P 3 is very close to 0. It means that the molecular distribution there was nearly not stretched out. As a result, the distribution of the molecules at P 3 remain in line with the normal distribution. However, the values of τ xx component at P 1 , and P 5 are nearly the same and both are bigger than 0. Therefore, the molecular distribution there are stretched and the molecules are distributed symmetrically. And so do P 2 and P 4 .

Conclusions
The BCF method provides a new perspective to study the rheological characteristics of polymer solutions using a multi-scale numerical solver. In this paper, we gave a full multi-scale model, which includes the macroscopic equations governing the hydrodynamics of the Newtonian solvent, and the microscopic equations governing the motion of the polymer molecules in the solution.
To solve the governing equations numerically, a multi-scale solver was implemented based on the open source CFD (Computional Fluid Dynamics) toolbox, OpenFOAM. The architecture of the solver and the numerical techniques used in the solver were discussed in detail. The simple two-dimensional Poiseuille flow with various configurations was run as the baseline. The key variable in the simulation, such as the velocity component u x , achieved good consistency with macroscopic counterpart. What's more, we designed a new scalar Field D to quantitatively describe the local stretch of molecules, thus the global correlations between the molecular distributions and the macroscopic polymer solutions could be analyzed.
The multi-scale simulation results of the OpenFOAM-based solver showed great potential to study complex engineering problems involving dilute polymer solutions. In the future, we will apply this approach to three-dimensional complex geometries. In order to reduce the computation time, we will transport the multi-scale solver to high performance clusters and focus on parallel optimization techniques.

Conflicts of Interest:
The authors declare no conflict of interest.