A multi-scale method for complex flows of non-Newtonian fluids

1.   Introduction

Modelling and computational simulation of non-Newtonian fluids is a challenging problem, since these fluids exhibit complex effects, such as shear thinning or thickening, viscoelasticity or flow-induced phase separation. A detailed analysis of the rheology of complex fluids can be obtained by particle-based simulations. Clearly, such micro-scale description is accurate but computationally very expensive and cannot be applied to engineering-scale problems. Consequently, new mathematical algorithms and hybrid multi-scale approaches have been proposed in recent years.

In the literature, we can find several hybrid models combining particle dynamics with a macroscopic continuum model, such as the heterogeneous multi-scale methods ^{[7,8,9,11,15,34]}, the triple-decker atomistic-mesoscopic-continuum method ^[13,14], the seamless multi-scale methods ^[10,12,30], the equation-free multi-scale methods ^[19,20,21] or the internal-flow multi-scale method ^[1,2]. In ^[27] software requirements and design principles are presented and illustrated for a prototype coupling between molecular dynamics and Lattice Boltzmann methods. Note that, in general, hybrid multi-scale approaches are successful when processes occurring on a small scale are only loosely coupled with the large-scale behavior, that is, in the presence of an effective scale separation. In keeping with this, our method is applicable when the local flow conditions experienced by a fluid element change on a much larger time scale than the microscopic relaxation time necessary to achieve a statistically steady state of the molecular conformation and interactions.

The aim of this paper is to illustrate a new heterogeneous multi-scale method for the simulation of flows of complex fluids in generic geometries and different conditions of motion. We extend and generalize the idea presented in ^[32], where a heterogeneous multi-scale method is developed for polymeric solutions subjected to simple shear flows. Simple shear is very well investigated for it is more easily realized in experiments than other flows. It provides fundamental information that is often sufficient to characterize simple fluids. However, complex fluids show a molecular structure that is able to change with respect to different conditions of motion, geometries, but also time scales and deformation rates. Therefore, to retrieve important rheological information about the stress response of such fluids, their behavior is studied under conditions of flow different from steady simple shear, for example extensional motions, startup flows and oscillating shear ^[26]. Moreover, it is also fundamental, for these types of materials, to consider flows in complex geometries that contain holes, barriers, contractions or other irregular features, because in close proximity of these structures the fluid is subject to local motions that are not equivalent to simple shear and can manifest some unexpected behavior.

We focus our efforts on the development of a method to link the micro-scale rheological information, available from simulations in different local flow types, to macro-scale flows in complex geometries. The main aim of this work is to show that it is necessary to correctly take into account the local flow type, which plays an important role in modifying the stress response in non-viscometric flows of non-Newtonian fluids featuring a flow-type dependent rheology, such as polymeric fluids. While simulations based on constitutive assumptions do take into account, by construction, the possible flow-type dependence incorporated at the level of constitutive laws, multi-scale methods for fluid mechanics application have so far considered simple shear flows as the only source of information for the micro-to-macro coupling. With our approach, we overcome this severe limitation. The method is applied to planar flows, but it can be extended to more general three-dimensional flows. The most challenging task for such an extension will be to implement MD simulations under arbitrary local flow conditions.

The paper is organized in the following way. Section 2 discusses the essential question that arises in building a heterogeneous multi-scale method. Namely, how micro- and macroscopic models are linked together. Section 3 is devoted to the description of micro-scale simulations with an emphasis on the decomposition of the average stress tensor and the reconstruction of data-driven material functions. For the macroscopic simulation of incompressible flows we use, as described in Section 4, a numerical method based on a semi-implicit time evolution scheme and mixed $\mathsf{P}_2$-$\mathsf{P}_0$ finite elements for the approximation of the velocity-pressure pair under the incompressibility constraint. Results of numerical simulations obtained by our heterogeneous multi-scale method in different geometries are presented in Section 5. The paper is closed by Section 6 with an outlook on future research.

2.   Link between micro-scale data and macro-scale simulations

At the heart of our proposal there is the understanding that the local kinematics of a generic planar flow must be identified by at least two parameters, one measuring the speed of the deformation and the other measuring the local flow type (see ^[24],Chapter 1.4] for a discussion of this point). We introduce the symmetrized gradient $\boldsymbol{\mathsf{D}}\equiv\tfrac{1}{2}(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^ \mathsf{T})$ of the velocity field $\boldsymbol{u}$ and, for the case of planar flows in which $\boldsymbol{e}_3$ is normal to the flow plane, we can define the kinematic parameters

While there is ample agreement about measuring the deformation rate in generic situations by means of $\dot{\varepsilon}$, different choices can be made for the flow-type parameter, depending on the type of material response that we envision (see ^[17], Sec. Ⅲ). Here we choose $\beta_3$ for its simplicity and we stress that, to make it a frame-indifferent parameter, we consider it as measured relative to the vanishing spin of a fixed inertial frame of reference. In other contexts, different choices may be appropriate (such as the relative rotation rate proposed by Schunk and Scriven ^[31]), but the structure of the scheme that we are going to present remains the same.

The local kinematic parameters $\dot{\varepsilon}$ and $\beta_3$ are the macro-scale input of micro-scale simulations that will in turn give an ensemble-averaged stress tensor $\boldsymbol{\sigma}$. To be able to use the micro-scale information encoded in $\boldsymbol{\sigma}$ we must express it in a form that is independent of the coordinates, because coordinates that are chosen for computational convenience cannot always be coherent when we pass from the micro- to the macro-scale simulation. We can achieve this independence by decomposing the stress on an orthogonal tensorial basis built upon $\boldsymbol{\mathsf{D}}$ (see ^[17] for details and significance of this decomposition). For the case of planar incompressible flows, we have

where $\boldsymbol{\mathsf{I}}$ is the identity tensor and

The coefficients of the linear decomposition (2.1) represent a pressure, a generalized viscosity and a reorientation factor related to normal stress differences and can be retrieved from the computational data in any coordinates by computing

Here the double scalar product and the tensorial norm are defined by

For the sake of clarity, we indicate with $\tilde{\eta}$ and $\tilde{\lambda}_3$ the two material functions of the kinematic parameters $(\dot{\varepsilon}, \beta_3)$, reconstructed by sampling the two-dimensional parameter space. In this way, we obtain from computational rheological measurements the constitutive expression

for the stress tensor, that can be used in performing the macro-scale continuum simulations.

There is no need to define a pressure function because it is determined, at the macro-scale, as the Lagrange multiplier of the incompressibility constraint.

3.   Micro-scale simulation

The NEMD simulation at the micro-scale is performed using LAMMPS [29,lammps.sandia.gov]. We consider both monomeric and polymeric aggregates in a three-dimensional computational domain with periodic boundary conditions and undergoing two different planar flows: simple shear and extensional flow. The average velocity field is imposed by suitably deforming the computational box and the canonical ensemble statistics of the velocity fluctuations is achieved via the classical Nosé–Hoover thermostat ^[16]. For the case of simple shear we use the LAMMPS command $\verb|fix deform|$ in conjunction with the SLLOD algorithm ^[5], while the extensional flow is treated by means of the UEF package ^[28]. The latter implements boundary conditions developed by Dobson ^[6] and Hunt ^[18], that are an extension of the boundary conditions of Kraynik and Reinelt ^[22]. We will discuss results for monomers, $N_\mathrm{pol} = 1$, and polymers, $N_\mathrm{pol} = 20$, where $N_\mathrm{pol}$ is the number of monomers per molecule.

3.1.   Interaction potentials

We model polymers as bead-spring systems with two types of interaction potential. The first one is active between each pair of monomers and is the Weeks-Chandler-Andersen (WCA) potential, the repulsive part of a Lennard-Jones potential with minimum energy $-\epsilon$ when the centers of the beads are at distance $r = \sqrt[6]{2}\sigma$:

The second one, only active between successive monomers in a polymeric chain, is the Finitely Extensible Nonlinear Elastic (FENE) potential

with $K$ the elastic constant and $R_0$ the maximum extent of the bonds.

3.2.   Average stress coefficients

We considered values of the deformation rate $\dot{\varepsilon}$ ranging from 0.001 s$^{-1}$ to 1 s$^{-1}$. For each value of $\dot{\varepsilon}$ we performed 21 simulations with different initial configurations of the particles. At each time, the average stress given as output by LAMMPS is projected onto the basis tensors through the formulae (2.3) and these projections are averaged over time. Finally, a mean over all initial configurations is taken. Figure 1 presents the data obtained for the generalized viscosity $\eta$ and the reorientation factor $\lambda_3$. By fitting those data with Carreau and power functions we obtain, for the monomeric case with $N_\mathrm{pol} = 1$, the curves

For the polymeric case with $N_\mathrm{pol} = 20$, the Gaussian Process Regression method (see Appendix A) was used to obtain the fitting curves.

By comparing the results obtained under simple shear and planar extension for $N_\mathrm{pol} = 1$ and $N_\mathrm{pol} = 20$, we find that the former fluid displays almost no rate dependence and very little flow-type dependence of the stress coefficients, while the latter features rate dependence and a strong flow-type dependence (Figure 1). The statistics of interaction between monomers in the two types of flow does not differ significantly, leading to a quasi-Newtonian rheology. On the contrary, the fact that polymer chains are kept in an elongated state by extensional flows, while they tend to tumble and keep a spherical shape in simple shear, generates a richer phenomenology. In fact, the shear-thinning trend displayed in simple shear is qualitatively similar to that of monomeric fluids, while we find an opposite trend---a shear-thickening behavior—in extensional flows at small deformation rates, followed by a non-monotonic behavior of the viscosity. The elastic effects generated by the bonds are crucial in this case and give rise, in simple shear, to the normal stress differences associated with the coefficient $\lambda_3$, also featuring a non-monotonic curve.

The direct connection between the statistics of chain conformation in the polymeric case and the observed rheology can be evinced from the probability distribution, presented in Figure 2, of the normalized radius of gyration $R_g$ and asphericity $\alpha$ defined as follows. Let the components of the $3\times3$ gyration tensor be

with $m_i$ the mass of the $i$-th monomer, $\boldsymbol{r_i}$ its position, $M = \sum_im_i$, and $r_\mathrm{cm}$ the center of mass of the polymeric chain. With $R_0$ the effective monomer radius, we set

Moreover, the asphericity is given by

where $a_1\leq a_2\leq a_3$ are the eigenvalues of the gyration tensor $\boldsymbol{\mathsf{R}}$.

At low deformation rate, the distributions obtained in simple shear and extensional flows are almost identical (blue curves in Figure 2) because the chains are equally and mildly stretched by the imposed flow. As $\dot{\varepsilon}$ grows larger, the extensional flow distributions feature a single peak that is progressively shifted rightward, toward more stretched configurations. On the other hand, in simple shear we notice a broadening of the distribution eventually leading to a doubly peaked shape (red curves in Figure 2). This indicates that chains spend roughly half of the time in a stretched configuration while frequently going back to a coiled configuration. In this very different microscopic behavior we find the origin of the strong flow-type dependence of the viscosity at larger $\dot{\varepsilon}$.

3.3.   Data-driven material functions

Ideally, we would like to perform micro-scale simulations with many values of $\dot{\varepsilon}$ and $\beta_3$, since flows in generic geometries can easily feature variations of the flow type with $\beta_3$ ranging from minus to plus infinity and values of $\dot{\varepsilon}$ that cover many orders of magnitude. The implementation of simulation algorithms for mixed flows different from simple shear and planar extension is a nontrivial task that will be addressed in future works, but we still need, for our macro-scale simulations, a suitable definition of material functions that covers the whole range of kinematic parameters. We thus apply simple extrapolation strategies to build material functions out of the available computational data.

Specifically, we extend the fitted curves $\eta_\mathrm{sh}(\dot{\varepsilon})$, $\eta_\mathrm{ext}(\dot{\varepsilon})$, $\lambda_\mathrm{3, sh}(\dot{\varepsilon})$ and $\lambda_\mathrm{3, ext}(\dot{\varepsilon})$ out of their natural domains by taking constant values. Then we set

that is an even function of $\beta_3$ with $\eta_\mathrm{ext}(\dot{\varepsilon}) = \tilde{\eta}(\dot{\varepsilon}, 0)$ and $\eta_\mathrm{sh}(\dot{\varepsilon}) = \tilde{\eta}(\dot{\varepsilon}, 1)$, and

that is an odd function of $\beta_3$ and such that $\lambda_\mathrm{3, ext}(\dot{\varepsilon}) = \tilde{\lambda}_3(\dot{\varepsilon}, 0) \approx 0$ and $\lambda_\mathrm{3, sh}(\dot{\varepsilon}) = \tilde{\lambda}_3(\dot{\varepsilon}, 1)$. These operations in combination with (2.4) allow us to build a complete set of data for the stress tensor as a function of $\dot{\varepsilon}$ and $\beta_3$.

4.   Macro-scale simulation

We consider the flow of an incompressible fluid with constant density $\rho$ and velocity vector $\boldsymbol{{u}}$. From the macroscopic point of view the continuum description leads to the conservation of mass and the balance of linear momentum. In the incompressible case, they translate into the following partial differential equations on the domain $\Omega$:

Here $\boldsymbol{\mathsf{{\boldsymbol{\sigma}}}}$ denotes the Cauchy stress tensor and $p$ the pressure field. Equation (4.1) expresses the incompressibility constraint. We partition the boundary of the domain as disjoint union of inlet $\Gamma_\mathrm{in}$, outlet $\Gamma_\mathrm{out}$, and solid walls $\Gamma_\mathrm{w}$, so that $\partial \Omega = \Gamma_\mathrm{in}\cup \Gamma_\mathrm{out}\cup \Gamma_\mathrm{w}$. As usual, $\boldsymbol{n}$ is the outward unit normal to $\partial \Omega$.

A weak formulation of the problem is retrieved by multiplying Eq (4.2) by a test function $\boldsymbol{v} \in H^1_{\Gamma_\mathrm w}(\Omega) = \big\{\boldsymbol{f} \in H^1(\Omega) : \boldsymbol{f}\vert_{\Gamma_\mathrm{w}} = \boldsymbol{0}\big\}$ and integrating over $\Omega$. By applying Green's formula, we obtain the integral equation

Moreover, (4.1) can be multiplied by a scalar test function $q \in L^2_0(\Omega) = \big\{g\in L^2(\Omega):\int_\Omega g = 0\big\}$ and integrated over $\Omega$ to give

We decompose the Cauchy stress tensor in spherical and deviatoric parts by introducing the traceless extra stress $\boldsymbol{\tau}$ such that $\boldsymbol{\sigma} = -(p+\bar{p}){\boldsymbol{\mathsf{I}}} + \boldsymbol{\tau}$, where $\bar{p}\in L^2_0(\Omega)$ is a given pressure field used to impose a chosen pressure gradient from inlet to outlet. On top of the Dirichlet boundary conditions for $\boldsymbol{u}$ on $\Gamma_\mathrm{w}$ we assume a form of periodicity for $\boldsymbol{u}$ and $p$, requiring that they take the same values on $\Gamma_\mathrm{in}$ and $\Gamma_\mathrm{out}$. Thanks to the independence of the test functions $\boldsymbol{v}$ and $q$, we can take the sum of (4.3) and (4.4) and substitute the expression for $\boldsymbol{\sigma}$ to obtain the complete weak formulation of the problem:

Find $\boldsymbol{u} \in L^2(\mathbb{R}^+; H^1_{\Gamma_\mathrm w}(\Omega))$ and $p \in L^2(\mathbb{R}^+; L_0^2(\Omega)),$ such that for any suitable test functions $\boldsymbol{v}$ and $q$ we have

Equation (4.5) must be completed with suitable initial conditions for the velocity field, that we will assume to vanish identically at time $t = 0$. To exploit the reflection symmetry of some domains, we will consider a fictitious boundary $\Gamma_\mathrm{c}$ corresponding to the symmetry axis. On this boundary we will assume that the normal component of the velocity field and the tangential component of the traction $\boldsymbol{\tau}\boldsymbol{n}$ vanish.

For the time integration of (4.5) we apply a semi-implicit Euler scheme. In particular, the nonlinear convective term is approximated explicitly, the rest is approximated implicitly. This approximation is suitable to simulate flows at sufficiently low Reynolds number. Consequently, we have

with bilinear forms $b(\boldsymbol{{u}}, q): = -\int_{\Omega} q \nabla\cdot{\boldsymbol{u}}$ and $a(\boldsymbol{{u}}, \boldsymbol{{v}}): = \int_{\Omega}{\boldsymbol{\tau}(\boldsymbol{u})}:\nabla{\boldsymbol{v}}.$ Note that the latter is already linearized since the (nonlinear) material functions in (2.4) are computed using $\boldsymbol{{u}}^n.$

For each time step, Eq (4.6) is discretized in space and solved, in a standard way, by means of mixed finite elements $\mathsf{P}_2$-$\mathsf{P}_0$ for the approximation of the velocity-pressure pair, which are known to satisfy the Ladyzhenskaya–Babuška–Brezzi inf-sup condition for a stable solution of the associated saddle-point problem ^[3]. The numerical method is implemented within the Python library FEniCS ^[25] and some details pertaining the micro-to-macro coupling through the dynamic reconstruction of the extra stress are given in Appendix B.

5.   Numerical results for different geometries

The aim of this section is to demonstrate the robustness and effectiveness of our heterogeneous multi-scale method by presenting the results of simulations in three different planar geometries: flows in a straight channel, through a 4:1 contraction, and past a deep hole. These are paradigmatic geometries frequently used to test the non-Newtonian behavior of fluid models ^[4].

5.1.   Channel flow

The multi-scale method is firstly tested on the classical case of a planar channel flow. The total length of the channel is $L = {5}\;{\rm{m}}$ along the $x$-axis and the system width is $W = {2}\;{\rm{m}}$ along the $y$-axis. Exploiting the reflection symmetry of the problem, we discretized only the upper half of the channel on a mesh with $100$ triangular elements along the $x$-axis and $20$ elements along the $y$-axis, for a total of $4000$ triangles in the entire computational domain. On the fictitious boundary that corresponds to the center of the channel we impose a vanishing vertical velocity and a vanishing normal traction to respect the reflection symmetry of the problem. We consider a fluid with unit density and drive the flow by imposing a constant horizontal pressure gradient $C = {1}\;{\rm{Pa/m}}$. As for the velocity field $\boldsymbol{u}$, we assume periodic boundary conditions so that the velocity at inlet and outlet is the same. The time step is $\Delta t = {1 \times 10^{-3}}\;{\rm{s}}$. We let the system evolve up to the time $T = {2}\;{\rm{s}}$ to reach a steady state.

The flow type is everywhere equivalent to that of a simple shear. Indeed, this is a viscometric flow with $\beta_3 = 1$ in the top half of the domain and $\beta_3 = -1$ in the other half. A discontinuity would appear at the midline of the channel consistent with the fact that $\beta_3$ is not defined for $\dot{\varepsilon} = 0$. The homogeneity of the flow type makes the flow-type dependence of the stress tensor irrelevant. To assess the importance of non-Newtonian effects, we compare the solution obtained from the multi-scale approach with the analytical solution of a reference Newtonian model, represented by the classical stress-strain relation

with $\bar{\eta}$ constant. We take $\bar{\eta} = {1.7}\;{\rm{Pa}} \;{\rm{s}}$ for $N_\mathrm{pol} = 1$ and $\bar{\eta} = {10}\;{\rm{Pa}} \;{\rm{s}}$ for $N_\mathrm{pol} = 20$, corresponding to the zero-rate limit of the interpolated shear viscosity from our NEMD simulations.

The comparison of velocity, deformation rate, and viscosity profiles is reported in Figure 3. In the case $N_\mathrm{pol} = 1$ the non-Newtonian solution is extremely close to the Newtonian one, as expected from the mild rate dependence of the viscosity and the small values of $\lambda_3$. In the case $N_\mathrm{pol} = 20$, we observe a much larger deviation that is small only around the center of the channel, where the shear rate vanishes and the viscosity approaches that of the reference Newtonian model.

Due to the non-uniformity of the viscosity in the development of the flow and in its steady state, the characterization of the flow by means of dimensionless quantities such as the Reynolds number is not straightforward. Nevertheless, we can identify the range of local parameter values as follows. We take the channel width $W$ as reference length and $P = CW$ as reference pressure drop. The local Reynolds number is then given by

Another important dimensionless quantity is the reorientation angle $\varphi$, such that

which measures the rotation of the stress eigenvectors with respect to the eigenvectors of $\boldsymbol{\mathsf{D}}$. This geometric information is precisely the meaning of the first normal stress difference in simple shear flows, extended to generic local flows by the definition of $\lambda_3$ and $\eta$ ^[17].

In the monomeric case, we have $\eta_\mathrm{min} = {1.6}\;{\rm{Pa}} \;{\rm{s}}$ and $\eta_\mathrm{max} = {1.8}\;{\rm{Pa}} \;{\rm{s}}$ so that $\mathrm{Re}\in[1.57, 1.77]$, while $\varphi$ is negligible. In the polymeric case, we have $\eta_\mathrm{min} = {2}\;{\rm{Pa}} \;{\rm{s}}$ and $\eta_\mathrm{max} = {10}\;{\rm{Pa}} \;{\rm{s}}$ so that $\mathrm{Re}\in[0.28, 1.41]$. We can then consider our examples to be at low Reynolds number. Since $\dot{\varepsilon}_\mathrm{max} = {0.24}\;{{\rm{s}}^{-1}}$, we have $\varphi\in[-22, 22]$ degrees, that includes quite significant values.

5.2.   Flow through a contraction

As a first example of flow in a complex geometry we consider the 4:1 contraction. This flow is characterized by the presence of different local flow types: simple shear in regions that are sufficiently far from the ends of the contraction, extensional and mixed motion in proximity of entrance and exit of the contraction, and mixed rotational flow near the corners of the domain. The spatial non-uniformity of the flow type must be taken into account to obtain a precise macro-scale description of the fluid motion. To show this, we compare the prediction of the full non-Newtonian model with that of a modified non-Newtonian model, wherein the flow-type dependence is artificially suppressed. Specifically, $\beta_3$ is replaced by the function $\mathrm{sign}(\beta_3)$. In this way, only the micro-scale data obtained in simple shear (i.e., $\beta_3 = \pm 1$) are used, but the results cannot properly reflect the fluid behavior.

We discuss results obtained with a channel of length $L = {5}\;{\rm{m}}$ and maximum width $W = {2}\;{\rm{m}}$, only in the polymeric case $N_\mathrm{pol} = 20$, with time step $\Delta t = 10^{-3}$, letting the system evolve to reach a steady state for $T = {0.2}\;{\rm{s}}$. We can appreciate the time convergence to the steady solution by looking at the evolution of the flow rate $Q$ through vertical sections of the domain, reported in Figure 4, left panel. The discretization employs an unstructured triangular mesh with average diameter of the triangles $h\approx{0.03}\;{\rm{m}}$. The mesh convergence of the multi-scale method is assessed by considering the $L^2$ norm of the difference between solutions obtained with mesh parameter $h$ and $2h$. From the data reported in Figure 4, right panel, we see that the incremental correction scales linearly in $h$, thus showing convergence of our simulation.

The pressure gradient that drives the flow is $C = {0.3}\;{\rm{Pa/m}}$. The velocity field $\boldsymbol{u}$ is assumed periodic so that it is matched at inlet and outlet and, to exploit symmetry and compute the numerical solution only in the upper half of the domain, we use the same conditions as in the straight channel simulation, namely, vanishing vertical velocity and normal traction at the center of the channel. The flow is laminar with the characteristic increase of horizontal velocity $u_x$ in correspondence of the contraction. The pressure gradient is concentrated along the narrower portion of the domain and features oblique isolines due to the presence of normal stress differences in that region (Figure 5). The non-uniformity of the deformation rate $\dot{\varepsilon}$ and of the flow-type parameter $\beta_3$ has a strong influence on the local values of $\eta$, that range from ${3.8}\;{\rm{Pa}} \;{\rm{s}}$ to ${13}\;{\rm{Pa}} \;{\rm{s}}$ in the contracting region, and $\lambda_3$, that range from ${-4}\;{\rm{Pa}} \;{\rm{s}}$ to ${4}\;{\rm{Pa}} \;{\rm{s}}$, considering also the lower half of the domain (see Figure 6). From these values and taking now $W/4$ as reference length, we can estimate $\mathrm{Re}\in[0.015, 0.05]$ and $\varphi\in[-23.5, 23.5]$ degrees.

In this complex flow, the flow type measured by the parameter $\beta_3$ varies significantly through the domain (Figure 6b), ranging from pure extension ($\beta_3 = 0$) to simple shear ($|\beta_3| = 1$) to more rotational flows ($|\beta_3| > 1$). Regions of mixed and extensional motion cover a major part of the domain, going well within the contraction. Since the viscosity in the full non-Newtonian model is a function of $\dot{\varepsilon}$ and $\beta_3$, the effect of the flow type is clearly dominant in determining its value. The importance of correctly embedding the flow-type dependence in the multi-scale model can be highlighted by comparing the prediction of the full non-Newtonian model with that obtained from the modified non-Newtonian model (Figure 7). The viscosity in the modified model presents a very different profile in the extensional and mixed flow regions, where it happens to decrease instead of increasing, as predicted by the correct method and in line with the extensional rheology of the polymeric case. This has an immediate effect on the velocity profile, that features a more rapid flow and a clear asymmetry, due to the overemphasized role of the normal stress differences related to $\lambda_3$ in the modified model.

5.3.   Flow past a deep hole

In this section we analyze a second example of flow in a complex geometry. We consider a channel of length $L = {5}\;{\rm{m}}$ and width $W = {1}\;{\rm{m}}$. At a distance of ${1.75}\;{\rm{m}}$ from the inlet we find the hole of width $W_\mathrm{hole} = {1.5}\;{\rm{m}}$ and depth $H = {5}\;{\rm{m}}$. The discretization employs an unstructured triangular mesh with average diameter of the triangles $h\approx{0.06}\;{\rm{m}}$. The time step is $\Delta t = 10^{-3}$ and we let the system evolve to reach a steady state for $T = {0.4}\;{\rm{s}}$ (see Figure 4, left panel). We present simulation results for the polymeric case $N_\mathrm{pol} = 20$, driving the flow with an outlet-to-inlet pressure gradient $C = {0.3}\;{\rm{Pa/m}}$. The velocity field is again assumed to be periodic, and hence equal at inlet an outlet.

The flow past a deep hole features a complex distribution of local flow types with a consequent variation of the material response. Also in this case, a comparison of the viscosity and the horizontal velocity obtained with the full and the modified non-Newtonian models (Figure 8) confirms the need for a correct treatment of flow-type-dependent rheologies. The viscosity appears to be significantly different right above the deep hole. This has a mild but noticeable influence on the flow. The velocity depression is shifted rightward and, since the viscosity is lower, the flow is globally faster with the modified model, as seen also from the values of the flow rate reported in Figure 4.

From the component $u_y$ of the velocity (Figure 9a) we can clearly see the region where the deep hole modifies the flow in the upper channel. Moreover, the presence of a clockwise-rotating vortex can be inferred at the center of the shown portion of the hole. The deformation rate $\dot{\varepsilon}$ and the flow-type parameter $\beta_3$ highlight the presence of mixed flows and high-vorticity regions (Figure 9b and 9c). Specifically, $\beta_3$ confirms the presence of a vortex in the hole. The sharp contrast between yellow and green regions nearby the openings of the upper channel indicates the viscometric nature of the flow far from the hole. The values of $\lambda_3$ across the domain (Figure 9d) are a direct consequence of the non-Newtonian rheology and the complexity of the flow. The same applies to the viscosity shown in Figure 7c. For this flow, the relevant ranges of dimensionless parameters are $\mathrm{Re}\in[0.05, 0.12]$ and $\varphi\in[-20.8, 20.8]$ degrees. Hence, we have significant non-Newtonian effects while turbulent motions are prevented.

6.   Conclusions and outlook

We have developed and tested a heterogeneous multi-scale method for the simulation of flows in arbitrary geometries for complex fluids by means of a data-driven reconstruction of the stress tensor. This is based on employing the results of micro-scale simulations according to a general decomposition of the stress ^[17] that allows to correctly align its different components in flows that display a broad spectrum of local flow types. We have shown that a proper treatment of the flow-type dependence, that was missing in previous multi-scale methods exclusively focused on simple shear rheology, is essential to capture the macroscopic dynamics.

In our case, the micro-scale simulations are non-equilibrium molecular dynamics simulations of polymeric chains with internal FENE bonds and an effective Weeks-Chandler-Andersen interaction potential. Nevertheless, our coupling scheme does not rely on a specific simulation strategy neither at the micro-scale nor at the macro-scale. Indeed, we have chosen to implement the macroscopic continuum simulation with a well-established mixed finite element method based on $\mathsf{P}_2$-$\mathsf{P}_0$ elements for the velocity-pressure pair to deal with the incompressibility constraint and a semi-implicit time integration scheme, but any other computational method to solve for the continuum evolution can in principle be employed.

Polymeric fluids are well known to display a significant degree of flow-type dependence in their rheological properties. In fact, the different average conformation of the long molecular chains in simple shear and extensional flows causes the same material to show opposite trends in the dependence of the viscosity on the rate of deformation. Our micro-scale simulations reproduce this feature. Hence, we were able to assess the relevance of properly taking into account the flow-type dependence of the stress in the macro-scale simulation. From the results obtained with our multi-scale approach in paradigmatic geometries, it is clear that the predictions of a model that suppresses the flow-type dependence are not reliable.

We performed micro-scale simulations imposing simple shear and planar extension as average kinematic conditions under which we evaluate the stress response. From these, we constructed the rheological functions in mixed flows by means of a simple interpolation and extrapolation procedure. In view of the impact on macroscopic flows in non-viscometric geometries of the local flow conditions, we plan to expand the micro-scale simulations to obtain data under mixed and rotational flow conditions. This will require the implementation of a suitable generalization of the Kraynik–Reinelt boundary conditions.

In fluids composed by molecules with longer microscopic relaxation times, it may happen that the local flow type experienced by a fluid element changes fast along streamlines. The method that we presented can in principle be extended to take this into account by devising MD simulations for which the flow type of the imposed background motion is slowly varying. The macroscopic material functions would then depend also on the derivative along streamlines of the flow-type parameter. Other directions for future work are the testing of our method in parameter ranges where elastic instabilities may be observed, that may lead to the presence of highly oscillatory material coefficients in space and time, and the implementation of an on-the-fly coupling between macro-scale and micro-scale simulations, to be able to sample the space of kinematic conditions in a problem-driven fashion.

Acknowledgments

G. G. G. and F. T. acknowledge the support of the Italian National Group of Mathematical Physics (GNFM-INdAM) through the funding scheme GNFM Young Researchers' Projects 2020. The research of M. L.-M. and L. Y. was supported by the German Science Foundation (DFG) within the Project number 233630050-TRR 146, C5 project. They also gratefully acknowledge the computing time granted on the HPC cluster Mogon at Johannes Gutenberg University Mainz. M. L.-M. is grateful to the Gutenberg Research College of the Johannes Gutenberg University Mainz for the research support.

Conflict of interest

The authors declare no conflict of interest.

A.   Fitting procedure for the micro-scale NEMD data

We report here the details about the fitting of micro-scale data that we performed using Gaussian Process Regression (GPR). Since some choices are involved in the procedure, it seems appropriate to give a brief account of ours. GPR is a Bayesian approach to regression based on machine learning techniques. The fitting is done through a Gaussian stochastic process $f(x; {\bf{w}})$ that evolves along the variable $x$ and depends on a list of hyper-parameters ${\bf{w}}$. The observed data are a finite number $N$ of pairs $(x_i, y_i)_{i = 1}^N$. The functional form of the Gaussian process is specified by its mean and covariance functions, $m({x}; {\bf w}) = \mathbb{E}[f({x}; {\bf w})]$ and $K({x}, { x'}; {\bf w}) = \mathbb{E}\big[\big(f({x}; {\bf w})-m({x}; {\bf w})\big)\big(f({x'}; {\bf w})-m({x'}; {\bf w})\big)\big]$, where $\mathbb E$ denotes the expected value. A fundamental feature of GPR is that it does not produce a definite value of the fitting parameters ${\bf{w}}$, but considers them random variables and seeks to determine suitable mean and variance for their distribution. This is achieved by an iterative optimization that takes into account the likelihood of the parameters distributions given the knowledge of the observed data ^[33].

What one needs to specify is the type of mean and covariance functions to be used for the GPR. Typical choices for the former are polynomial functions of $x$ with coefficients given by ${\bf{w}}$, while for the latter are Radial Basis Functions

with a set of two hyper-parameters ${\bf w} = (\sigma, l)$, or the Matérn covariance of degree $\nu$ given by

where $\Gamma$ is the Euler gamma function and $K_\nu$ is the modified Bessel function of the second kind. In our treatement, we took as $x$ variable the quantity $\dot{\varepsilon}$ or $\log_{10}(\dot{\varepsilon})$, as $y$ variable we took $\eta$, $\log_{10}(\eta)$, or $\lambda_3$. We tested regressions with mean function zero, constant, or polynomial up to degree four and, for the covariance function, RBF or Matérn kernel of degree up to 5 and picked the combinations that give the best results after optimization.

B.   Multi-scale coupling in the code

The purpose of this section is to present the portion of code that implements the multi-scale coupling in the context of an otherwise standard Finite Element simulation of the flow equations (see ^[23] for an introduction to the library FEniCS). The peculiarity of our method concerns the integration of NEMD data with the continuum solver and the construction of the stress tensor based on the general representation (2.4).

In our code, $\verb|numpy|$ arrays contain the interpolated curves $\eta_{\mathrm{sh}}$, $\lambda_{3, \mathrm{sh}}$ and $\eta_{\mathrm{ext}}$, while $\lambda_{3, \mathrm{ext}}$ is simply zero. In particular, the first column of the array $\verb|S|$ contains the values of $\dot{\varepsilon}$ on which the Carreau function or GPR is sampled and the second column the corresponding expected values of $\eta_\mathrm{sh}$. For an efficient use of those data in the stress computation, we transform them in projections on a Finite Element space defined on a one-dimensional mesh, the nodes of which are given by the sampled values of $\dot{\varepsilon}$. Considering, for instance, the curve $\eta_\mathrm{sh}$, it becomes the function $\verb|etas|$ in the space $\verb|E|$ of $\mathsf{P}_1$ elements defined on the mesh $\verb|mesh_etas|$. Its construction is performed with the following lines of code. Similar commands are used to define the functions $\verb|etae|$ and $\verb|lambda3s|$ that represent $\eta_\mathrm{ext}$ and $\lambda_{3, \mathrm{sh}}$.

The second portion of code relevant to the multi-scale coupling is the local computation of the kinematic parameters $\dot{\varepsilon}$ and $\beta_3$ and the reconstruction, at each time step, of the spatial profile of the stress coefficients $\eta$ and $\lambda_3$ by means of the material functions $\tilde{\eta}(\dot{\varepsilon}, \beta_3)$ and $\tilde{\lambda}_3(\dot{\varepsilon}, \beta_3)$. This is accomplished through the following definitions, where $\verb|mesh|$ is the two-dimensional triangular mesh defined on the flow domain, and $\verb|u|$ is a three-component variable containing the velocity and pressure fields.

At each time step, we need to compute the new stress coefficients based on the velocity field $\verb|u_n|$ obtained in the previous step. These will be used to update the stress. To this end, we apply the interpolation between simple shear and extensional data according to (3.1)–(3.2).