# Deterministic particle approach of Multi Bead-Spring polymer models

Amine Ammar, Francisco Chinesta, David Ryckelynck

<figure><img src="image_1.png" /></figure>

## To cite this version:

Amine Ammar, Francisco Chinesta, David Ryckelynck. Deterministic particle approach of Multi Bead-Spring polymer models. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2006, 15 (5), pp.481-494. 10.3166/remn.15.481-494. hal-01006737

HAL Id: hal-01006737

https://hal.science/hal-01006737

Submitted on 6 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Deterministic particle approach
of Multi Bead-Spring polymer models

Amine Ammar* — Francisco Chinesta**

David Ryckelynck**

*Laboratoire de Rhéologie, UMR 5520 CNRS – INPG – UJF
1301 rue de la piscine, BP 53 Domaine universitaire, F-38041 Grenoble cedex 9
Amine.Ammar@ujf-grenoble.fr

**LMSP UMR 8106 CNRS-ENSAM-ESEM
151 boulevard de l'Hôpital, F-75013 Paris
{francisco.chinesta, david.ryckelynck}@paris.ensam.fr

**ABSTRACT.** Kinetic theory models of complex fluids involve the resolution of the advection-diffusion 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.

**RÉSUMÉ.** Les fluides complexes décrits dans le cadre de la théorie cinétique ont généralement été traités à l'aide des techniques stochastiques. Ils ont permis, entre autres, la simulation du comportement des haltères de type FENE (Keunings, 1997; 2004). La difficulté majeure introduite par les modèles Multi-Haltères des polymères est reliée au caractère multidimensionnel de l'espace des conformations résultant, qui nécessite un grand nombre de trajectoires stochastiques (Somasi et al., 2002). Dans un travail antérieur (Ammar and Chinesta, 2005), une approche déterministe a été proposée pour traiter le modèle cinétique des suspensions de fibres courtes où, après introduction du terme de diffusion dans celui d'advection, l'équation de convection pure résultante est discrétisée à l'aide d'une technique de particules. Dans ce travail, cette approche déterministe est appliquée à la simulation du modèle Multi-Haltère pour des écoulements homogènes afin d'évaluer sa pertinence ainsi que son impact sur le nombre de particules pour de grandes dimensions de l'espace des conformations.

**KEYWORDS:** Fokker-Planck equation, Multi Bead Spring models, method of particles.

**MOTS-CLÉS :** équation de Fokker-Planck, modèles Multi-Haltères, méthodes particulaires.
# 1. 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.

# 2. 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><img src="image_3.png" /><figcaption>**Figure 1.** MBS model for polymer chain</figcaption></figure>

The dynamics of the chain is governed by hydrostatic, Brownian and connector forces. If we denote by $\dot{r}_j$ the velocity of the bead 'j' and by $\dot{q}_j$ the velocity of the spring connector, then we obtain:

$$ \dot{q}_j = \dot{r}_{j+1} - \dot{r}_j, \quad \forall j = 1, \dots, N \qquad [1] $$

The dynamics of each bead can be written as:

$$ \underbrace{-\zeta(\dot{r}_j - v_0 - \text{grad } v \cdot r_j)}_{\text{Hydrostatic effects}} - k_b T \underbrace{\frac{\partial}{\partial r_j} \ln(\Psi)}_{\text{Brownian effects}} + \underbrace{F_j^c - F_{j-1}^c}_{\text{Interactions Forces}} = 0 \qquad [2] $$

where $\zeta$ 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 $\Psi$ is the probability distribution function. From Equations [1] and [2] it results:

$$ \dot{q}_j = \text{grad } v \cdot q_j - \frac{1}{\zeta} \sum_{k=1}^{N} A_{jk} \left( k_b T \frac{\partial}{\partial q_k} \ln(\Psi) + F_k^c \right) \qquad [3] $$

where $A_{jk}$ is the Rouse Matrix defined by:
$$A_{jk} = \begin{cases} 2 & \text{if } j=k \\ -1 & \text{if } j=k \pm 1 \\ 0 & \text{otherwise} \end{cases}$$

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

$$\frac{\partial \Psi(q_1, \dots, q_N, t)}{\partial t} = - \sum_{j=1}^{N} \left( \frac{\partial}{\partial q_j} (\dot{q}_j \Psi(q_1, \dots, q_N, t)) \right) \quad [4]$$

and introducing Equation [3] it results:

$$\begin{aligned} \frac{\partial \Psi}{\partial t} = & - \sum_{j=1}^{N} \left( \frac{\partial}{\partial q_j} \left( \left( \text{grad} \cdot q_j - \frac{1}{\zeta} \sum_{k=1}^{N} A_{jk} F_k^c \right) \Psi \right) \right) + \\ & + \frac{k_b T}{\zeta} \sum_{j=1}^{N} \sum_{k=1}^{N} A_{jk} \frac{\partial^2 \Psi}{\partial q_j \partial q_k} \end{aligned} \quad [5]$$

Moreover, Equation [3] can be rewritten excluding Brownian effects which are retained in Equation [5]. Then, Equations [3] and [5] can be expressed, in a dimensionless form (see Bird et al., 1987, for more details), as:

$$\underline{\dot{q}} = \text{grad} \cdot \underline{q} - \frac{1}{4} \underline{\underline{A}} \underline{\underline{F}}^c(\underline{q}) \quad [6]$$

$$\frac{\partial \Psi}{\partial t} = - \frac{\partial}{\partial \underline{q}} (\dot{\underline{q}} \Psi) + \frac{1}{4} \frac{\partial}{\partial \underline{q}} \left( \underline{\underline{A}} \frac{\partial \Psi}{\partial \underline{q}} \right) \quad [7]$$

Where $\underline{q}$ denotes the vector containing the $N$ spring connector vectors $q_j$, with $j = 1, \dots, 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:

$$F^c(q_j) = \frac{1}{1 - q_j^2 / b} q_j \quad [8]$$

where $\sqrt{b}$ is the maximum stretching of each spring connector of the chain.

# 3. Deterministic discretization approach using the method of particles

## 3.1. 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 pseudo-advection field:

$$\tilde{\underline{q}} = \dot{\underline{q}} - \frac{1}{4} A \frac{\partial \Psi / \partial \underline{q}}{\Psi} \quad [9]$$

being the associated Fokker-Planck equation:

$$\frac{\partial \Psi}{\partial t} = - \frac{\partial}{\partial \underline{q}} (\tilde{\underline{q}} \Psi) \quad [10]$$

Problem [9]-[10] is equivalent to the initial problem [6]-[7].

## 3.2. 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 $\Psi = \sum_{i}^{Cn} \alpha_{i} \delta(\underline{q} - \underline{q}^{i})$, where 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:

$$ \zeta_{\epsilon}(d) = \frac{\zeta\left(\frac{d}{\epsilon}\right)}{\epsilon} \quad \text{with} \quad \zeta\left(\frac{d}{\epsilon}\right) = \frac{e^{-\left(\frac{d}{\epsilon}\right)^2}}{\sqrt{\pi}} \qquad [11] $$

where $\epsilon$ denotes a smoothing parameter and ‘d’ denotes a distance in the $M \times N$-dimensional space of molecule configurations.

Thus we can write:

$$ \zeta_{\epsilon}\left(\|\underline{q}-\underline{q}^{i}\|\right)=\zeta_{\epsilon}(d^{i})=\frac{\zeta\left(\frac{d^{i}}{\epsilon}\right)}{\epsilon} \qquad [12] $$

where $d^i = \|\underline{q} - \underline{q}^i\|$ is the Euclidean norm in the $M \times N$-dimensional space.

The smoothed distribution can then be written as:

$$ \Psi(\underline{q}) = \sum_{i=1}^{C_n} \alpha_i \zeta_{\epsilon}(d^i) \qquad [13] $$

Where $\alpha_i$ is the weight related to each chain, which, in this work, is considered the same overall the population ($\alpha = 1/C_n$ for satisfying the normality condition of the distribution function).

### 3.3. Equilibrium conditions

Before starting simulations, an equilibrium state $\Psi_{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 $\epsilon$ by comparing $\Psi_{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 $2 \times N$ dimensional case considered in this work, the expression of $\Psi_{eq}$ (for the connector force defined in Equation [8] is given by:

$$ \Psi_{eq}(\underline{q}) = \frac{\prod_{j=1}^{N} \left( 1 - \frac{q_j^2}{b} \right)^{b/2}}{\left( 2\pi b / (2+b) \right)^N} \qquad [14] $$
### 3.4. Simulation algorithm

The simulation algorithm proceeds as follows:

- Consider $C_n$ particles (molecules) consisting of $N$ connectors whose orientation must describe the initial known distribution:

$$q_j^{i(0)}, \quad i = 1, \dots, C_n; \quad j = 1, \dots, N$$

that is, $q_i^{i(0)}$ is assumed known.

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 [9]:

$$\dot{\underline{q}}^{i(n)} = \dot{\underline{q}}^{i(n)} - \frac{1}{4} A \frac{\partial \Psi^{(n)}}{\partial \underline{q}|_{\underline{q}^{i(n)}}}$$

where $\Psi^{(n)}(\underline{q}) = \sum_{i=1}^{C_n} \alpha_i \zeta_{\varepsilon}(d^{i(n)}),$ being $d^{i(n)} = \|\underline{q} - \underline{q}^{i(n)}\|$.

- Update the molecular conformation: $\underline{q}^{i(n+1)} = \underline{q}^{i(n)} + \dot{\underline{q}}^{i(n)} \Delta t$

- Update the distribution function: $\Psi^{(n+1)}(\underline{q}) = \sum_{i=1}^{C_n} \alpha_i \zeta_{\varepsilon}(d^{i(n+1)})$

### 4. 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:

$\Delta t = 0.001, \varepsilon = 0.02, b = 50, We = du/dx = 10$. At the initial time the particles are distributed for verifying the equilibrium initial condition given by:
$$ \Psi(q) = \left( \frac{1}{1 - q^2 / b} \right)^{-b/2} \qquad [15] $$

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><img src="image_4.png" /><figcaption>Figure 2. Evolution of the distribution function related to a single 1D FENE model</figcaption></figure>

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$). Similar results were obtained by using only around 600 particles.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Zoom of the distribution functions in the region close to the maximum
moldeule stretching</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Histogram related to the particle position at t=1</figcaption></figure>

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 following parameters: $C_n=500$ (number of particles), $b=10$ ($\sqrt{b}$ being the dimensionless parameter giving the maximum stretching), $N = 2$ (number of branches), $M=1$ (dimension of the physical space) and $We = du/dx = 0.35$.

Figure 5 depicts the reference steady solution as well as the solution computed by using the method of particles with a smoothing parameter $\epsilon = 0.4$. The steady distribution function depends on the stretching of both branches $q_1$ and $q_2$.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Distribution function related to a 2 branches FENE MBS molecule model conformation: (left) Reference solution; (right) Solution computed using the method of particles</figcaption></figure>

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 $\epsilon = 0.1, 0.2, 0.3, 0.5$ ($\epsilon = 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.
<figure><img src="image_8.png" /><figcaption>Figure 6. Particles distribution related to the solution depicted in figure 5(right)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Conformation distribution of a 2 branches FENE MBS model related to $\epsilon = 0.1$ (a), $\epsilon = 0.2$ (b), $\epsilon = 0.3$ (c) and $\epsilon = 0.5$ (d)</figcaption></figure>
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$, $du/dy = 1$, $\Delta t = 0.01$, $N = 5$, $M=2$, $C_n = 52$, $\epsilon = 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.

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:

$$ \Psi_j = \sum_{i=1}^{C_n} \alpha_i \zeta_{\epsilon} (\|q - q_j^i\|) = \frac{1}{C_n} \sum_{i=1}^{C_n} \alpha_i \zeta_{\epsilon} (\|q - q_j^i\|), \quad j = 1, \dots, N \quad [16] $$

where $q_j^i$ refers the connector $j$ of chain $i$.

<figure><img src="image_10.png" /><figcaption>Figure 8. MBS conformation in a simple shear flow</figcaption></figure>
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 $\Psi_1$ and $\Psi_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.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** *Example of chain representation at time t2*</figcaption></figure>

**5. 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.

**6. References**

Ammar A., Chinesta F., A Particle Strategy for Solving the Fokker-Planck Equation Governing the Fiber Orientation Distribution in Steady Recirculating Flows Involving Short Fiber Suspensions, Lectures Notes on Computational Science and Engineering, Springer, 43, 2005, p. 1-16.

Bird R.B., Crutiss C.F., Armstrong R.C., Hassager O., *Dynamic of polymeric liquid*, vol. 2, *Kinetic Theory*, John Wiley & Sons, 1987.
Chaubal C.V., Srinivasan A., Egecioglu O., Leal L.G., "Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems", Part 1, *Method Journal of Non-Newtonian Fluid Mechanics*, 70, 1997, p. 125-154.

Chauviere C., Lozinski A., Simulation of dilute polymer solutions using a Fokker-Planck equation. *Computers and Fluids*, 33, 2004, p. 687-696.

Chinesta F., Chaidron G., Poitou A., "On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions", *Journal of Non-Newtonian Fluid Mechanics*, 113/2-3, 2003, p. 97-125.

Keunings R., "On the peterlin approximation for finitely extensible dumbbells", *Journal of Non-Newtonian Fluid Mechanics*, 68, 1997, p. 85-100.

Keunings R., "Micro-Macro Methods for the Multiscale Simulation of Viscoelastic Flow using Molecular Models of Kinetic Theory", *Rheology Reviews* 2004, D.M. Binding and K. Walters (Eds.), *British Society of Rheology*, 2004, p. 67-98.

Lozinski A., Chauviere C., "A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model", *J. Computational Physics*, 189, 2003, p. 607-625.

Somasi M., Khomami B., Woo N. J., Hur J. S., "E.S.G. Shaqfeh. Brownian dynamics simulation of bead-rod and bead-spring chains: numerical algorithms and coarse graining issues", *Journal of Non-Newtonian Fluid Mechanics*, 108, 2002, p. 227-255.