Deterministic particle approach of Multi Bead-Spring polymer models

Kinetic theory models of complex fluids involve the resolution of the advectiondiffusion Fokker-Planck equation that is usually performed using stochastic approaches. Stochastic simulation for Finitely Extensible Non Elastic (FENE) dumbbells has been successfully applied (Keunings 1997; Keunings 2004). The main difficulty found in the simulation of Multi Bead-Spring (MBS) polymer models using that approach is related to the high number of realizations required because of the highly multidimensional conformation space (Somasi et al., 2002). In a former work (Ammar and Chinesta 2005), a deterministic approach was proposed for treating the kinetic theory description of short fibers suspensions which operates by introducing the diffusion term into the advection one, then applying the method of particles. In this work, that deterministic approach will be applied to the Multi Bead-Spring polymer models in simple flows in order to conclude about its performance as well as about the impact of the conformation space dimension on the number of particles that must be considered in the simulation.


Introduction
As a consequence of the increasing use of composite materials, micro-macro simulation comes back into fashion nowadays with the development of computer performance. However, it remains crucial that simplicity and performance of numerical techniques affect the usability of microscopic models.
In the case of dilute polymer, simple dumbbell models, such as the FENE model, are not representative of the behaviour of a polymeric chain. In fact, both, the flexibility of the chain and the description of its complex configurations state are better accounted using a Multi Bead-Spring (MBS) representation.
Simulation of MBS has been widely performed by using stochastic techniques, whose accuracy is related to the number of realizations considered. In the context of deterministic techniques, the use of standard discretization techniques such as finite differences or finite elements are forbidden, because of the extremely large number of degrees of freedom involved due to the multidimensional character of the associated Fokker-Planck equation. Thus, the direct resolution of the Fokker-Planck equation in the finite element framework only concerned simple dumbbell models (Lozinski and Chauviere 2003;Chauviere and Lozinski 2004). An alternative to the standard discretization techniques of the Fokker-Planck equation consists of the use of the method of particles. However, this technique has never been used until now for solving models defined in a multidimensional space, as in the case of the MBS polymer models. The main aim of this work is the analysis of the application of this technique to those models.
Results of such simulation highlight the flexibility effects of the chain (especially in shearing flows), which cannot be observed in the case of simple dumbbell models. However, a number of difficulties persist, some of them restricting seriously its application in some simulation conditions.

Mechanical model
As shown in Figure 1, the MBS chain consists of N+1 beads connected by N springs. The bead serves as an interaction point with the solvent and the spring contains the local stiffness information depending on local stretching (see Bird et al., 1987, for more details).

Figure 1. MBS model for polymer chain
The dynamics of the chain is governed by hydrostatic, Brownian and connector forces. If we denote by j r the velocity of the bead 'j' and by j q the velocity of the spring connector, then we obtain: The dynamics of each bead can be written as: Hydrostatic effects Interactions Forces Brownian effects where ζ is the drag coefficient, v is the velocity field, v 0 is an average velocity, k b is the Boltzman constant, T is the absolute temperature and Ψ is the probability distribution function. From Equations [1] and [2] it results: where A jk is the Rouse Matrix defined by: In the Rouse model connector force F c is a linear function of the connector orientation vector. But in our case a Finitely Extensible Non Linear spring will be used so that each spring branch behaves as in FENE model whose connector force will be detailed later.
Using now the evolution of the distribution function given by Where q denotes the vector containing the N spring connector vectors j q , with 1 j N = , , .
The problem [6]-[7] is defined in a MxN-dimensional space, being M=1,2 or 3 (the dimension of the physical space) and N the number of springs that compose the polymeric chain. In the present work the connector force is given by: where b is the maximum stretching of each spring connector of the chain.

Advection form of the Fokker-Planck equation
For applying the method of particles to discretize the Fokker-Planck equation, we firstly introduce its diffusion term into the advection one (Ammar and Chinesta 2005; Chinesta et al., 2003;Chaubal et al., 1997). Thus, it results the pseudoadvection field: being the associated Fokker-Planck equation: Problem [9]-[10] is equivalent to the initial problem [6]-[7].

Smoothed particle approximation
In order to discretize Equations [9]-[10] one must define the approximation of the unknown field which is, in our case, the distribution function. The method of particles uses a Dirac representation, consisting of ( ) the superscript refers to each chain of the population involved in the simulation that consists of Cn chains (particles). However, from a computational point of view the derivative of Dirac masses cannot be defined in a proper way. So, usually, a smoothed approximation of the Dirac mass is introduced according to: where ε denotes a smoothing parameter and 'd' denotes a distance in the MxNdimensional space of molecule configurations.
Thus we can write: The smoothed distribution can then be written as: Where α i is the weight related to each chain, which, in this work, is considered the same overall the population (α = 1/Cn for satisfying the normality condition of the distribution function).

Equilibrium conditions
Before starting simulations, an equilibrium state Ψ eq can be defined, and then used for defining the initial condition for start-up problems, as well as for adjusting the value of the smoothing parameter ε by comparing Ψ eq with the solution computed from the steady state related to the transient solution of the Fokker-Planck equation without the advection term.
According to Bird et al., (1987) and after the application of the normality condition in the 2xN dimensional case considered in this work, the expression of Ψ eq (for the connector force defined in Equation [8] is given by: SPH in polymer flow simulations 487

Simulation algorithm
The simulation algorithm proceeds as follows: -Consider Cn particles (molecules) consisting of N connectors whose orientation must describe the initial known distribution: Index i refers to the chain (particle) involved in the simulation, j to the connector in the chain, and (n) indicates the time step.
-For each time step: -Compute the advection field according to Equation

Numerical test
In order to validate the approach just proposed we consider firstly a problem whose solution can be obtained analytically. This is for example the case of a single FENE dumbbell in a one-dimensional physical space. This problem was deeply analysed in (Keunings, 1997). The number of particles involved in the simulation is in the order of 2400, being the other simulation parameters the following: that is equivalent to Equation [14] after normalization. We prefer using in this numerical example Equation [15] as initial condition instead of Equation [14] because in this manner the computed solution can be directly compared with the results reported in (Keunings, 1997).
The evolution of the distribution function is depicted in Figure 2 which represents the solutions at times: t=0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,1. Figure 3 zooms the region located around the maximum molecule stretching. These distributions are in perfect agreement with the exact solution as well as with the numerical solutions computed using the stochastic technique that can be found in (Keunings, 1997). Figure 4 depicts the histogram related to the particle position at t=1 (we recall that the portion of the particle indicates the position of one of its ends, the other one remains at the origin x=0). Simular results were obtained by using only around 600 particles.  Now, after proving the accuracy of the proposed strategy, we consider higher dimensional problems. For this purpose we consider a MBS -FENE model defined in a one-dimensional physical space composed of two dumbbells. In this case a reference solution can be easily obtained by using a more experienced numerical technique as is the finite element method. We consider the problem defined by the  The distribution depicted in Figure 5(right) is related to the particle distribution depicted in Figure 6. In order to conclude about the incidence of the smoothing parameter on the computed solution we show in Figure 7 the distributions computed by using different values of that parameter, in particular 0.1,0.2,0.3,0.5 ε = ( 0.4 ε = being the optimum value whose associated distribution is depicted in Figure 5(right)). We can notice the significant impact of the smoothing parameter on the particles distribution and then on the conformation distribution. Moreover this parameter must be identified in each particular application because it depends on the number of particles considered in the simulation; on the solution itself (more localized is the solution smaller must be the smoothing parameter, but at the same time it must be large enough for coupling particles in the regions with low particles density); and on the dimension of the physical and conformation spaces. The definition of an appropriate adaptive procedure for evaluating the local value of the smoothing parameter is not an easy task. Finally, to analyse the capabilities of the approach just described, for treating more realistic models, a MBS -FENE model composed of many dumbbells in a 2D physical space, is considered. This test is carried out in a simple shear flow with the following simulation conditions: b = 9, 1 du dy = , 0.01 t ∆ = , N = 5, M=2, Cn = 52, ε = 0.8 and being the initial condition given by Equation [14]. In all the simulations the flow kinematics is assumed known and uncoupled with the polymer presence.

Figure 2. Evolution of the distribution function related to a single 1D FENE model
The solutions at times t1 = 3 and t2 = 10 (at t2 the steady state is reached) are depicted in Figure 8. To simplify the solution post-processing, a smoothed distribution involving the connectors belonging to the same level (having the same index in the chain) is defined according to: where i j q refers the connector j of chain i. From the computed results, we notice that symmetric levels (in our case levels 1 and 5, or 2 and 4) exhibit the same behaviour. Connectors are more stretched in the centre of the chain than on its ends. Figures 8(a) and 8(b) depicts Ψ 1 and Ψ 3 at time t1 and Figures 8(c) and 8(d) the same distributions at time t2. We can also notice the perfect symmetry of the computed distributions since the initial distributions respect symmetry. Figure 9 depicts the chain configurations at time t2 of some of the chains involved in the simulation. In this figure the origin of the first spring connector q 1 of all chains is located at the centre of the figure.

Conclusions
Despite of its simplicity, the deterministic particle technique just described seems to be robust and accurate when the chain model consists of a reduced number of connectors. However a number of difficulties persist, being the most important: -the number of particles increases significantly with the dimension of the conformation space as was also found in the context stochastic simulations; -the smoothing parameter has a strong impact on the computed results, even if it can be identified from the equilibrium conditions, when the distribution becomes highly localized an adaptation of this value is required to guaranty the solution accuracy. From some numerical experiments we have conclude that each particle must have a minimum number of neighbours that seems to be around 3; -it is well known that in highly multidimensional spaces the use of a Euclidean distance results delicate, because in a high dimensional space the measure of the surface of a hyper-sphere becomes higher than the measure of its volume.