Active- and transfer-learning applied to microscale-macroscale coupling to simulate viscoelastic flows

Relations between the stress tensor and kinematic measures of deformation for non-Newtonian polymeric fluids are generally complicated, and reflect the collective dynamical response of polymer micro-structures to macroscopic deformation [11]. In the present study, the DPD method is employed to explicitly simulate the dynamics of polymer micro-structures at a coarse-grained level. Similar to molecular dynamics (MD) systems, a DPD model consists of many interactive particles governed by the Newton’s equation of motion [12],

$$m_{i}\frac{{\rm d}^{2}\mathbf{r}_{i}}{{\rm d}t^{2}}=m_{i}\frac{{\rm d}\mathbf{v}_{i}}{{\rm d}t}=\mathbf{F}_{i}=\sum_{j\neq i}\left(\mathbf{F}^{C}_{ij}+\mathbf{F}^{D}_{ij}+\mathbf{F}^{R}_{ij}\right),$$

(1)

where $m_{i}$ is the mass of a particle $i$ and $t$ denotes time. The bold symbols $\mathbf{r}_{i}$ and $\mathbf{v}_{i}$ represent position and velocity vectors, and $\mathbf{F}_{i}$ is the total force acting on the particle $i$ due to the presence of neighboring particles. The summation for computing $\mathbf{F}_{i}$ is carried out over all other particles within a cutoff radius beyond which the pairwise interactions are negligible. The pairwise force is comprised of a conservative force $\mathbf{F}^{C}_{ij}$, a dissipative force $\mathbf{F}^{D}_{ij}$ and a random force $\mathbf{F}^{R}_{ij}$, which respectively are of the forms:

$$\begin{split}&\mathbf{F}^{C}_{ij}=a_{ij}\omega_{C}(r_{ij})\mathbf{e}_{ij},\\ &\mathbf{F}^{D}_{ij}=-\gamma_{ij}\omega_{D}(r_{ij})(\mathbf{v}_{ij}\mathbf{e}_{ij})\mathbf{e}_{ij},\\ &\mathbf{F}^{R}_{ij}\cdot{\rm d}t=\sigma_{ij}\omega_{R}(r_{ij}){\rm d}\tilde{W}_{ij}\mathbf{e}_{ij},\end{split}$$

(2)

where $r_{ij}=|\mathbf{r}_{ij}|=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ is the distance between two particles $i$ and $j$, $\mathbf{e}_{ij}=\mathbf{r}_{ij}/r_{ij}$ is the unit vector, and $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$ is the difference of their velocities; ${\rm d}\tilde{W}_{ij}$ is an independent increment of the Wiener process [13]. The coefficients $a_{ij}$, $\gamma_{ij}$ and $\sigma_{ij}$ determine the strength of each force respectively. To satisfy the fluctuation-dissipation theorem (FDT) [13], and to maintain the DPD system at a constant temperature [12], the dissipative force and the random force are related by $\sigma_{ij}^{2}=2\gamma_{ij}k_{B}T$ and $\omega_{D}(r_{ij})=\omega^{2}_{R}(r_{ij})$. In the classical DPD model, all the forces in Eq. (2) have the same finite interaction range $r_{c}$ and their amplitudes decay with the distance $r_{ij}$ according to corresponding weight functions $\omega_{C}(r_{ij})=1-r_{ij}/r_{c}$, $\omega_{D}(r_{ij})=\omega^{2}_{R}(r_{ij})=(1-r_{ij}/r_{c})^{1/2}$ with $r_{c}$ being the cutoff radius.

The micro-structures of polymer molecules are represented as bead-spring chains so that their dynamics can be explicitly simulated when their suspensions are subject to flows. The chain structure of a polymer is modeled as $N$ beads connected by a finitely extensible nonlinear elastic (FENE) spring [14]. The FENE potential is given by

$$U_{B}(r)=\left\{\begin{matrix}-\frac{1}{2}k_{s}r_{0}^{2}\ln\left[1-(r/r_{0})^{2}\right]\ ,&r\leq r_{0}\\ \infty\ ,&r>r_{0}\\ \end{matrix}\right.$$

(3)

where $k_{s}$ is the spring constant and $r_{0}$ the maximum bond extension of the FENE spring. When chosen to be short and stiff enough the nonlinear FENE bond has been shown to minimize phantom bond-crossings [15] during chain entanglements.

The polymer solution is then represented by many bead-spring chains suspended in a solvent of free DPD particles, a widely used model in studies of polymer solutions [16, 17, 18, 19]. DPD is a coarse-grained method whose input parameters must be prescribed, and whose macroscopic rheological properties, such as viscosity and diffusivity, are outputs. However, it is worth noting that the DPD model is grounded in atomistic dynamics and can be derived directly from MD systems using the rigorous mathematical framework of the Mori-Zwanzig projection [20]. In our demonstration we calibrate the DPD model of polymer solution with experimental values measured at various shear strain rates by Purnode and Crochet [5] for a 0.5% aqueous polyacrylamide solution. We start with the classical DPD model with constant cutoff radius for pairwise interactions. Specifically, we consider a DPD system with a number density of $\rho_{n}=3.0$ at $k_{B}T=1.0$, and set the parameters in Eq. (2) as $a_{ss}=a_{pp}=25.0,~{}a_{ps}=12.5,~{}\gamma_{ij}=4.5,\sigma_{ij}=3.0,r_{c}=1.0$, where the subscripts ‘$s$’ represents solvent and ‘$p$’ for polymer. The polymer solution consists of $756$ polymer chains suspended in $88\,200$ DPD solvent particles in a computational box of $60\times 70\times 10$ (in reduced DPD units). Simple-shear DPD simulations at various shear strain rates are performed to obtain the shear stress and the first normal stress difference computed by the Irving-Kirkwood formula [21].

The DPD results are mapped to the corresponding measured physical ones in the experimental data by choosing appropriate characteristic length, mass and time units. For the standard DPD model with constant cutoff radius, this mapping yields $[L]=1.138\times 10^{-3}~{}m$, $[M]=4.92\times 10^{-7}~{}kg$ and $[T]=5.43\times 10^{-3}~{}s$ as the characteristic length, mass and time units, respectively. Figure 2 presents the rheological properties of the polymer solution obtained by the standard DPD model of bead-spring chains suspended in a solvent of free particles with $r_{c}=1.0$ compared to experimental data for a 0.5% aqueous polyacrylamide solution [5]. It is obvious from Fig. 2 that the DPD model with constant cutoff radius $r_{c}=1.0$ correctly captures the shear-thinning feature of polymer solution, however, the magnitude of viscosity variation in the DPD model is one-order smaller than the experimental data, as shown in Fig. 2(a). The apparent reason appears to be that the DPD model is a coarse-grained representation of the polymer molecules which does not model their atomistic structures. Each bead in a DPD chain represents a cluster of many monomers, where the changes of unresolved polymer structures due to shear are not captured well by the DPD model. Polymer coils become stretched out under increasing shear strain rate. The range of shear strain rates for stretch-out is far wider for real polymer coils than that found in simulations with the coarse-grained polymer model. This may explain why the latter fails to capture correct shear-thinning effect. To remedy this, we take the cutoff radius to be a function of the second invariant of the strain rate tensor ${\bm{\epsilon}}$, in the form of

$$r_{c}(\dot{\gamma})=r_{c,0}+\kappa\log\left(\frac{\chi}{\dot{\gamma}}\right),$$

(4)

where $\dot{\gamma}=\sqrt{2{\rm{tr}}({\bm{\epsilon}}^{2})}$ represents the second invariant of the strain rate tensor; $r_{c,0}$, $\kappa$ and $\chi$ are undetermined parameters. Eq. (4) assumes that the interaction range of polymer beads decreases with increasing shear strain rates. Because the transport properties such as viscosity and diffusivity of a DPD system are determined mainly by the dissipative and random forces, we set the $r_{c}$ in $\mathbf{F}^{D}_{ij}$ and in $\mathbf{F}^{R}_{ij}$ as a function of $\dot{\gamma}$ in Eq. (4), while we maintain the $r_{c}$ in $\mathbf{F}^{C}_{ij}$ constant as $r_{c}=r_{c,0}$. In particular, the variable cutoff radius $r_{c}(\dot{\gamma})$ is bounded by the region [1,2] to avoid unreasonable large or small cutoff radii, i.e., $r_{c}(\dot{\gamma})=\max[1,\min(2,r_{c}(\dot{\gamma}))]$. We set the parameters in Eq. (4) as $r_{c,0}=1.0$, $\kappa=0.265$ and $\chi=10.0$ in reduced DPD units to generate rheological properties consistent with the experimental data.

For the DPD model with variable cutoff radius, the reduced DPD units are mapped to physical units and compared to experimental data by setting $[L]=1.44\times 10^{-4}~{}{\rm m}$, $[M]=9.99\times 10^{-10}~{}{\rm kg}$ and $[T]=1.1\times 10^{-3}~{}{\rm s}$ as the characteristic length, mass and time units, respectively. Then, all other physical units of flow and rheological variables can be mapped to reduced DPD units by dimensional analysis using the three basic units $[T]$, $[L]$ and $[M]$, i.e., the shear strain rate is converted by the strain rate unit $[\dot{\gamma}]=1/[T]=903.22~{}{\rm s}^{-1}$. Also, the dynamic viscosity and the shear stress are converted by the factors $[\eta]=[M][L]^{-1}[T]^{-1}=6.304\times 10^{-3}~{}{\rm Pa\,s}$, $[\tau]=[M][L]^{-1}[T]^{-2}=5.694~{}{\rm Pa}$, and $[N_{1}]=2[M][L]^{-1}[T]^{-2}=11.39~{}{\rm Pa}$, respectively. Figure 2 shows the dependence of the non-Newtonian viscosity $\eta$ and the first normal stress difference $N_{1}$ on shear strain rate $\dot{\gamma}$ in steady shear flow. The DPD results with variable cutoff radius $r_{c}(\dot{\gamma})$ given by Eq. (4) are now consistent with the experimental data for a 0.5% aqueous polyacrylamide solution [5].

In a steady shear flow, analytical expressions for the FENE-P fluid were derived by Purnode and Crochet [6] for the strain rate dependent viscosity $\eta$, and the first normal stress difference $N_{1}$, as

$$\begin{split}\eta&=\eta_{s}+\eta_{p}\left(\frac{L^{6}}{4\lambda^{2}\dot{\gamma}^{2}\left(L^{2}-3\right)^{2}}\right)^{1/3}\left(\Delta_{1}^{1/3}+\Delta_{2}^{1/3}\right),\\ N_{1}&=\frac{\eta_{p}L^{4}}{\lambda\left(L^{2}-3\right)}\left(\frac{2L^{6}}{\lambda^{2}\dot{\gamma}^{2}\left(L^{2}-3\right)^{2}}\right)^{-1/3}\left(\Delta_{1}^{1/3}+\Delta_{2}^{1/3}\right)^{2},\end{split}$$

(5)

where $\eta_{s}$ and $\eta_{p}$ are the solvent and polymer viscosities respectively, $\lambda$ is the relaxation time, and $L$ is the ratio of the extended length of the spring to its equilibrium length. Hence $L^{2}$ is a measure of the extensibility of the FENE dumb-bells. $\Delta_{1}=1+[1+2L^{6}/(27\lambda^{2}\dot{\gamma}^{2}(L^{2}-3)^{2})]^{1/2}$, and $\Delta_{2}=2-\Delta_{1}$. Further details of the FENE-P fluid will be introduced in Section 2.2.

Purnode and Crochet [6] fitted their experimental data plotted in Fig. 2 using a FENE-P parameter set $\eta_{s}=0.012~{}{\rm Pa\,s},~{}\eta_{p}=0.378~{}{\rm Pa\,s},~{}L^{2}=4,~{}\lambda=2~{}s$. Although their non-Newtonian viscosities are well-fitted, the normal stress differences deviate significantly from the FENE-P curve, i.e., the opposite sign of the second derivative of the curves, which suggests that the 0.5% aqueous polyacrylamide solution is not a perfect FENE-P fluid. Likewise it also shows that the standard bead-spring DPD model does not yield the FENE-P stress even when chain lengths are varied substantially. The improved fits of both properties with the introduction of a variable cutoff radius $r_{c}(\dot{\gamma})$ dependent on a macro level measure of deformation rate indicates that modeling mesoscopic inter-molecular forces is a crucial factor in obtaining realistic macroscopic viscoelastic constitutive equations.

The flow of an incompressible fluid is governed by the conservation of mass and momentum [22],

$$\begin{split}\nabla\cdot\mathbf{u}&=0,\\ \rho\frac{\partial\mathbf{u}}{\partial t}+\rho\mathbf{u}\cdot\nabla\mathbf{u}=&-\nabla p+\nabla\cdot\bm{\tau},\end{split}$$

(6)

where $t$, $\rho$ and $p$ denote time, density and pressure, respectively; $\mathbf{u}$ is the velocity vector, and $\bm{\tau}$ the stress tensor. The momentum equation has to be closed by a specified constitutive equation. This work is motivated by the flows of polymer solutions, hence the chosen constitutive equation is FENE-P  [23, 24] given by,

$$\bm{\tau}=\bm{\tau^{s}}+\bm{\tau^{p}},~{}~{}~{}\bm{\tau}^{s}=\eta_{s}\cdot 2{\bm{\epsilon}},~{}~{}~{}\bm{\tau}^{p}=\frac{\eta_{p}}{\lambda}\frac{L^{2}}{L^{2}-3}\left(\frac{L^{2}-3}{L^{2}-tr\left(\bm{C}\right)}\bm{C}-\bm{I}\right),$$

(7)

where $\bm{\tau}^{s}$ is the stress generated by the solvent with $\eta_{s}$ being the solvent viscosity, ${\bm{\epsilon}}=[\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}]/2$ is the rate of deformation tensor, $\bm{\tau}^{p}$ is the stress generated by the polymer with $\eta_{p}$ being the polymer viscosity, $\lambda$ is the relaxation time, $L$ is the extensibility of the FENE dumbbell, and $\bm{C}$ is the conformation tensor computed from the following equation,

$$\overset{\triangledown}{\mathbf{C}}=-\frac{1}{\lambda}\frac{L^{2}}{L^{2}-3}\left(g\mathbf{C}-\mathbf{I}\right),~{}~{}~{}g=\frac{L^{2}-3}{L^{2}-tr\left(\mathbf{C}\right)},$$

(8)

in which $\overset{\triangledown}{\mathbf{C}}$ is the Oldroyd upper convected derivative defined by

$${\overset{\triangledown}{\mathbf{C}}=\frac{\partial\mathbf{C}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{C}-\mathbf{C}\cdot(\nabla\mathbf{u})-(\nabla\mathbf{u})^{T}\cdot\mathbf{C}.}$$

(9)

In order to improve the performance of the continuum model in reproducing the experimental observations, instead of using a classical FENE-P model with constant relaxation time and polymer viscosity, we propose a modified FENE-P model with polymer viscosity and relaxation-time dependent on the second invariant of strain rate by drawing an analogy to the White-Metzner model [25] with variable viscosity and relaxation time, i.e., $\eta_{p}(\dot{\gamma})$ and $\lambda(\dot{\gamma})$, which can be obtained from simple shear DPD simulations of polymer solution at various shear strain rates.

Eqs. (6)-(8) are non-dimensionalized by the characteristic length $L_{0}$, velocity $U_{0}$ and viscosity $\eta_{0}$. All physical variables are scaled accordingly as: $x^{*}=x/L_{0}$, $u^{*}=u/U_{0}$, $\eta^{*}=\eta/\eta_{0}$, $t^{*}=tU_{0}/L_{0}$, $p^{*}=p/(\rho U_{0}^{2})$ and $\tau^{*}=\tau L_{0}/(\eta_{0}U_{0})$. The non-dimensionalized Eq. (6)-(8) were solved using the SEM method [22] with a high-order time-splitting method with stiffly-stable scheme [26] for temporal discretization. Moreover, we employed the entropy viscosity method (EVM) [27] that introduces a nonlinear diffusion term to the FENE-P constitutive equation to stabilize the SEM simulation at high Weissenberg numbers (see details in Appendix A).

A polymer solution can exhibit non-Newtonian (shear-thinning) behavior, i.e., the shear viscosity and normal stress differences change with shear strain rate or shear strain rate history. The reason is that the polymer molecules suspended in solvent have chain like micro-structures, namely the length along the polymer chain is much larger than the other dimensions of the molecule. More importantly, the micro-structure of polymer is deformable and can be significantly changed by the flow conditions. The dissolved polymer molecules are generally in a random coil state in an equilibrium solution [16, 28]. However, under shear flows, the polymer chains trend to be aligned in the direction of shear and generate a reduced flow resistance as the shear strain rate increases. Consistently, the apparent viscosity of a polymer solution will decrease with increasing shear strain rate and exhibits a shear thinning behavior, while the normal stress differences increase as the shear strain rate increases. Using multiscale simulations, we couple the microscale and macroscale descriptions and directly connect the dynamics of polymer chains and the continuum equation for bulk rheology. Hence, a constitutive closure model can be on-the-fly constructed from the microscale dynamics using an active-learning algorithm rather than using an ad hoc empirical formula for the constitutive relation.

At the microscale level, the polymer solution is modeled by bead-spring chains suspended in a solvent of free DPD particles. Specifically, the DPD system is constructed with $756$ chains in a computational box of $60\times 70\times 10$ in reduced DPD units. Each chain consists of $50$ DPD beads connected by the FENE spring given by Eq. (3) with spring constant $k_{s}=30$ and $r_{0}=1.5$. The overall particle number density of the system is $\rho_{n}=3$ and the temperature is $k_{B}T=1.0$. The non-bonded interactions between DPD particles are computed by Eq. (2). For the conservative force, we use $a_{ss}=a_{pp}=25,a_{ps}=12.5$, and $\omega_{C}(r)=1-r/r_{c}$ with a cutoff radius $r_{c}=1.0$, where the subscript ‘$s$’ represents solvent and ‘$p$’ for polymer. For the dissipative force, we use $\gamma_{ss}=\gamma_{pp}=\gamma_{ps}=4.5$ and $\omega_{D}(r)=[1-r/r_{c}(\dot{\gamma})]^{1/2}$ with a variable cutoff radius $r_{c}(\dot{\gamma})$ given by Eq. (4). The parameters in the random force are given by $\sigma_{ij}^{2}=2\gamma_{ij}k_{B}T$ and $\omega^{2}_{R}(r)=\omega_{D}(r)$. We will perform simple shear DPD simulations at a few selected shear strain rates $\dot{\gamma}$ to obtain the stress tensor $\bm{\tau}$ at various shear strain rates, and we will compute the shear tensor is computed by the Irving-Kirkwood formula [21]. Let $\tau_{xy}$ be the shear stress and $N_{1}=\tau_{xx}-\tau_{yy}$ be the first normal stress difference, the apparent dynamic viscosity can be computed by $\eta=\tau_{xy}/\dot{\gamma}$. The shear strain rate dependent relaxation time of polymer solution is calculated from $\lambda=N_{1}/(2\dot{\gamma}^{2}\eta_{p})$ with $\eta_{p}$ being the polymer viscosity.

At the macroscale level, the bulk rheology of the polymer solution is modeled by the continuum equations given by Eqs. (6). We consider a 0.5% aqueous polyacrylamide solution according to Purnode’s experiments [5]. We take the dynamic viscosity $\eta_{0}=0.378~{}{\rm Pa\,s}$ at a shear strain rate of $\dot{\gamma}_{0}=0.504~{}{\rm s^{-1}}$ in the experimental data as the zero shear viscosity, and the normal stress difference at a median shear strain rate $\dot{\gamma}=17.625~{}{\rm s^{-1}}$ is $N_{1}=3.925~{}{\rm Pa}$. Based on these physical quantities, we can determine the length unit $[L]=1.44\times 10^{-4}~{}{\rm m}$, the mass unit $[M]=9.99\times 10^{-10}~{}{\rm kg}$, and the time unit $[T]=1.1\times 10^{-3}~{}{\rm s}$ of the DPD system. For hydrodynamic and rheology problems, all the physical quantities with physical units can be mapped to reduced DPD units by dimensional analysis using the three basic units $[L]$, $[M]$ and $[T]$. For example, the shear strain rate is converted by the shear strain rate unit $[\dot{\gamma}]=1/[T]=903.22~{}{\rm s}^{-1}$, the dynamic viscosity is converted by the viscosity unit $[\eta]=[M][L]^{-1}[T]^{-1}=6.304\times 10^{-3}~{}{\rm Pa\,s}$, the shear stress is converted by the stress unit $[\tau]=[M][L]^{-1}[T]^{-2}=5.694~{}{\rm Pa}$, and the first normal stress difference is $[N_{1}]=2[M][L]^{-1}[T]^{-2}=11.39~{}{\rm Pa}$.

Flow between two parallel plates, as shown in Fig. 1, is used to demonstrate the application of the active-learning scheme to microscale-macroscale coupling for modeling viscoelastic flows. At the macroscopic level, the bulk rheology of the polymer solution is modeled by the continuum equations in the form of Eq. (6), which are discretized and solved by SEM [26] with an entropy-viscosity method (EVM) [27] for numerical stabilization at high Weissenberg numbers. Let coordinates $x$ and $y$ denote the streamwise and transverse (wall-normal) directions, and $u$ and $v$ represent the corresponding velocity components, respectively. We perform the SEM simulation in a two-dimensional rectangular computational region of $20L_{0}\times 2L_{0}$, with two flat plates fixed at $y=\pm L_{0}$ as the wall boundaries, as shown in Fig. 1. The computational domain is discretized by $8\times 6$ spectral elements with a polynomial order of $p=7$. The inlet is at $x=-10L_{0}$, the outlet is at $x=10L_{0}$.

The channel flow is initialized with a zero velocity field. A pressure gradient is applied to drive the polymer solution to flow through the channel. The periodic boundary condition is applied in the $x$-direction, while the surfaces of the upper and lower plates are set to be solid walls with a no-slip ($u=0,v=0$) boundary condition. The Reynolds number is defined by $Re=UL/\nu_{0}=2\tilde{R}U/(2.14\tilde{U})$ with the channel width $L=2L_{0}$, the flow velocity $U$ and the zero-shear-rate viscosity $\nu_{0}=2.14\tilde{\nu}$. Flows at different $Re$ can be generated by applying various pressure gradients to drive the flow, or by changing the zero-shear-rate viscosity of the fluid. The elasticity of the polymer solution is described by the Weissenberg number, which is defined as the ratio of the microscopic time scale to the macroscopic strain rate, i.e., $We=\lambda_{0}U/L=0.5\tilde{W}U/\tilde{U}$.

In the microscale-macroscale coupled simulation, a FENE-P model, modified to have shear-thinning polymer viscosity and relaxation time, is used to close the continuum equations. They are discretized and solved by SEM on rectangular elements, while the microscale dynamics of polymer chains is simulated with DPD. The continuum solver acts as a master solver which actively calls the slave DPD simulator on-the-fly as necessary. A GPR model is employed to decide when and where to do the next sampling. To initialize the GPR model for active-learning, we first set up the SEM system and perform three DPD simulations of simple shear flow at low strain rates. The effective dynamic viscosity and relaxation time of the polymer solution can be derived from the computed stress tensor by $\eta=\tau_{xy}/\dot{\gamma}$ and $\lambda=(\tau_{xx}-\tau_{yy})/(2\dot{\gamma}^{2}\eta_{p})$. Three DPD simulations of the polymer solution in simple shear flow are performed to generate the initial training data for GPR, namely at strain rates in reduced units: $\dot{\gamma}=0.5~{}{\rm s^{-1}}$ ($5.5\times 10^{-4}$ in reduced unit), $\dot{\gamma}=1.0~{}{\rm s^{-1}}$ ($1.1\times 10^{-3}$ in reduced unit) and $\dot{\gamma}=1.5~{}{\rm s^{-1}}$ ($1.65\times 10^{-3}$ in reduced unit).

The GPR model predicts not only the mean apparent viscosity $\eta(\dot{\gamma})$ and normal stress difference $N_{1}(\dot{\gamma})$ at various strain rates, but also their corresponding uncertainties in terms of standard deviations $\sigma_{\eta}(\dot{\gamma})$ and $\sigma_{N_{1}}(\dot{\gamma})$. Because the range of the shear strain rate during the transient channel flow could change by several orders of magnitude, we perform the GPR on logarithmic scales, e.g., $\log\eta$ and $\log N_{1}$ versus $\log\dot{\gamma}$. Given the GPR prediction uncertainties, we can define an acquisition function $\alpha_{\eta}=\sigma_{\eta}(\dot{\gamma})$ for viscosity $\eta$ and an acquisition function $\alpha_{N_{1}}=\sigma_{N_{1}}(\dot{\gamma})$ for $N_{1}$ to indicate the location of next optimal sampling point. We note that the acquisition function $\alpha$ can also be defined by other functions of prediction uncertainty as well as its magnitude and curvature, wherein the selection of acquisition function would depend on the intrinsic smoothness of target functions in different problems. Let $\alpha=\max\{\alpha_{\eta},\alpha_{N_{1}}\}$ represent the maximum prediction uncertainty, we choose an acceptance criterion $\max(\alpha)\leq\delta_{\rm tol}$ for the GPR-informed effective viscosity and normal stress difference, where $\delta_{\rm tol}=1\%$ is a predefined tolerance.

The results of DPD simulation performed at three low shear strain rates ($\dot{\gamma}_{O_{1}}=0.5~{}{\rm s^{-1}}$, $\dot{\gamma}_{O_{2}}=1~{}{\rm s^{-1}}$ and $\dot{\gamma}_{O_{3}}=1.5~{}{\rm s^{-1}}$) are taken as the initial training points of a GPR model. Fig. 3(a1) presents the GPR predictions of dynamical viscosity $\eta$ and normal stress difference $N_{1}$ based on the three training points, where the solid line is the GPR prediction of $\eta$ and the dotted line is the prediction of $N_{1}$, and the training data for $\eta$ and $N_{1}$ are represented by circles and squares, respectively. The predicted uncertainties are visualized by the $95\%$ confidence interval shown as the shaded area. Fig. 3(a2) plots the magnitude of the corresponding acquisition function $\alpha_{\eta}(\dot{\gamma})$ and $\alpha_{N_{1}}(\dot{\gamma})$, where the predefined tolerance $\delta_{\rm tol}=1\%$ is shown by the horizontal dot-dashed line. It is obvious in Fig. 3(a2) that the magnitude of the acquisition function $\alpha=\max\{\alpha_{\eta},\alpha_{N_{1}}\}$ reaches its maximum value $\alpha=0.10$ at the shear strain rate $\dot{\gamma}=0.71~{}{\rm s^{-1}}$, which is far beyond the designed tolerance $\delta_{\rm tol}$ and breaks the acceptance criterion. Therefore, we need additional training points by running DPD simulations to reduce the magnitude of $\alpha$. In practice, we take the shear strain rate $\dot{\gamma}_{p_{1}}$ at the maximum $\alpha$ as the next sampling point. Then, we perform a simple shear DPD simulation of the polymer solution at $\dot{\gamma}_{p_{1}}=0.71~{}{\rm s^{-1}}$ and obtain the corresponding $\eta$ and $N_{1}$ derived from the computed stress tensor. With only one additional data point, we observe a better GPR prediction of $\eta$ and $N_{1}$ with smaller uncertainties as shown in Fig. 3(b1). Although the magnitude of $\alpha$ is significantly reduced by adding one data point at $\dot{\gamma}_{p_{1}}$, the maximum of $\alpha$ is still larger than the tolerance $\delta_{\rm tol}$. Similarly, we need to perform another simple shear DPD simulation at the next sampling point $\dot{\gamma}_{p_{2}}=1.26~{}{\rm s^{-1}}$, which is indicated by the location of maximum $\alpha$ in Fig. 3(b2). This process of adding new points can be repeated until the GPR predicted uncertainties in terms of $\alpha$ are smaller than the prescribed tolerance $\delta_{\rm tol}$. Fig. 5(a1) presents the GPR predictions of $\eta$ and $N_{1}$ based on the five data points, whose uncertainty described by $\alpha$ is smaller than the tolerance $\delta_{\rm tol}$ over the range of shear strain rate up to $\dot{\gamma}=1.86~{}{\rm s^{-1}}$.

Because the flow field of the SEM system is initialized with zero velocity, we can advance the time integration of the SEM solver using the GPR predictions of $\eta$ and $N_{1}$. The GPR model based on the five data points shown in Fig. 5(a1) is valid until the global maximum of $\dot{\gamma}_{\max}$ in the SEM simulation exceeds a value that breaks the acceptance criterion $\max(\alpha)\leq\delta_{\rm tol}$. In the start-up channel flow, the global maximum of the strain rate $\dot{\gamma}_{\max}$ in the SEM system changes with time before the flow reaches a steady state. Therefore, we need to monitor the magnitude of the acquisition function $\alpha(\dot{\gamma})$ (up to $\dot{\gamma}_{\max}$) in order to decide when we need to run additional DPD simulations. We can see in Fig. 5(a2) that the acquisition functions $\alpha_{\eta}$ and $\alpha_{N_{1}}$ increase sharply with the shear strain rate $\dot{\gamma}$. Let $\dot{\gamma}_{\rm gpr}$ be the maximum strain rate covered by a valid GPR-prediction $\eta$ and $N_{1}$ with $\alpha\leq\delta_{\rm tol}$. When the global maximum of strain rate $\dot{\gamma}_{\max}$ in the channel flow is larger than $\dot{\gamma}_{\rm gpr}$, we need additional training points for extrapolation to reduce the magnitude of $\alpha$. Unlike interpolation of new points in GPR, the extrapolation of the next sampling point $\dot{\gamma}_{p}$ requires the size of the jump to be determined. To that end, a tunable parameter $\beta$ is introduced to determine jump size for extrapolation, so that the extrapolation point is found from $\alpha(\dot{\gamma}_{p})=\beta\delta_{\rm tol}$. In general, too small a $\beta$ will fail to reduce the number of training points while a too large one may break the acceptance criterion, and require extra interpolating points. The optimum choice for $\beta$ achieves $\max(\alpha)\simeq\delta_{\rm tol}$ with only a single extrapolation sampling point.

To find a proper value of $\beta$ for extrapolation in the GPR, we perform a few numerical experiments using different values of $\beta$ as shown in Fig. 4. In the GPR of the effective viscosity $\eta$, we first try $\beta_{\eta}=5$ and obtain the next sampling point at $\dot{\gamma}_{p_{3}}=2.44~{}{\rm s^{-1}}$ by setting $\alpha(\dot{\gamma}_{p_{3}})=5\delta_{\rm tol}$. After adding one sampling point at $\dot{\gamma}_{p_{3}}=2.44~{}{\rm s^{-1}}$, we observe in Fig. 4(a) that $\max(\alpha)=0.47\%<\delta_{\rm tol}/2$, which means that $\beta_{\eta}=5$ results in an overcautious jump for extrapolation. Subsequently, we try a larger value of $\beta_{\eta}=10$, leading to the next sampling point at $\dot{\gamma}_{p_{3}}=3.23~{}{\rm s^{-1}}$. This leads to $\max(\alpha)=2.34\%>\delta_{\rm tol}$ and breaks the acceptance criterion, as shown in Fig. 4(a). Consequently, using $\beta_{\eta}=10$ may be too aggressive for extrapolation and we will need to add extra interpolating points to satisfy $\max(\alpha)|_{\dot{\gamma}<\dot{\gamma}_{p_{3}}}<\delta_{\rm tol}$. Based on the tests on $\beta_{\eta}=5$ and $\beta_{\eta}=10$ in Fig. 4(a), using $85\%$ of the value obtained by a linear interpolation, we estimate that $\beta_{\eta}=6.41$ could be a good option for extrapolation with the prescribed tolerance $\delta_{\rm tol}=1\%$. Using the same method, we test $\beta_{N_{1}}=10$ and $\beta_{N_{1}}=20$ for the GPR of the first normal stress difference $N_{1}$. We observe in Fig. 4(b) that $\beta_{N_{1}}=10$ results in $\max(\alpha)=0.56\%$ and $\beta_{N_{1}}=20$ results in $\max(\alpha)=1.25\%$, from which an optimal value for extrapolation of $N_{1}$ is $\beta_{N_{1}}=14.19$. To reduce the number of DPD simulations, we selected a proper value of $\beta$ to avoid performing extra interpolating samplings after each extrapolation, with the next sampling point determined by

$$\dot{\gamma}_{p}=\min(\dot{\gamma}_{\eta,p},\dot{\gamma}_{N_{1},p}),$$

(10)

where the magnitudes of $\dot{\gamma}_{\eta,p}$ and $\dot{\gamma}_{N_{1},p}$ are computed by $\alpha_{\eta}(\dot{\gamma}_{\eta,p})=\beta_{\eta}\delta_{\rm tol}$ and $\alpha_{N_{1}}(\dot{\gamma}_{N_{1},p})=\beta_{N_{1}}\delta_{\rm tol}$, respectively. Fig. 5 shows the extrapolating process in the GPR model to construct the effective viscosity $\eta$ and the first normal stress difference $N_{1}$. More specifically, Fig. 5(a2) and (b2) plot the acquisition functions $\alpha_{\eta}(\dot{\gamma})$ and $\alpha_{N_{1}}(\dot{\gamma})$, in which the next sampling points $\dot{\gamma}_{p3}=2.61~{}s^{-1}$ in Fig. 5(a2) and $\dot{\gamma}_{p4}=3.55~{}s^{-1}$ in Fig. 5(b2) are determined by $\alpha(\dot{\gamma})=\beta_{\eta}\delta_{\rm tol}$ using $\beta_{\eta}=6.41$. Fig. 5(c1) shows that the macroscale-microscale coupled system of a channel flow at $\tilde{R}=10$ and $\tilde{W}=1$ is eventually closed with 7 sampling points in total, and its inset shows the learned relaxation time $\lambda$ as a function of shear strain rate $\dot{\gamma}$. Fig. 8(a) plots the distribution of the shear strain rate $\dot{\gamma}$, the dynamic viscosity $\eta$ and relaxation time $\lambda$ in the channel flow at $\tilde{R}=10$ and $\tilde{W}=1$, showing that both the effective viscosity $\eta$ and the relaxation time $\lambda$ are decreased near the wall surfaces because of increased shear strain rate $\dot{\gamma}$. The inset of Fig. 8(a) shows that the global maximum of strain rate in the channel flow increases with time before reaching $\dot{\gamma}_{\max}=3.20~{}{\rm s^{-1}}$ at the steady state, while the GPR model for $\eta$ and $N_{1}$ based on seven training points is valid up to $\dot{\gamma}_{\rm gpr}=3.93~{}{\rm s^{-1}}$. Consequently, the GPR-informed $\eta$ and $N_{1}$ computed from the microscopic dynamics are sufficient to close the macroscopic equations with the modified FENE-P model in this channel flow.

Using the same driving force, we performed simulations using the classical FENE-P model (FENE-P) and also the macroscale-microscale coupled model (FENE-P+DPD). Fig. 6 shows different profiles of the polymer solutions flowing through a channel, including the velocity profile $u(y)$ in Fig. 6(a), the polymer normal stress profile $\tau^{p}_{xx}(y)$ in Fig. 6(b), and the polymer shear stress profile $\tau^{p}_{xy}(y)$ in Fig. 6(c). For the channel flow with $\tilde{R}=10$ and $\tilde{W}=1$, our SEM solutions for a classical FENE-P model with constant parameters are represented by the squares, which is consistent with the analytic solution of a FENE-P model (solid line) and verifies the SEM solver. The modified FENE-P model with shear-thinning viscosities $\eta(\dot{\gamma})$ and relaxation times $\lambda(\dot{\gamma})$ results in the macroscale-microscale coupling solutions represented by circles. It is obvious in Fig. 6(a) that the FENE-P+DPD model has a smaller flow rate than the classical FENE-P model in the channel flow. The FENE-P+DPD model has smaller viscosity than the classical FENE-P model because of shear thinning effects, as shown in Fig. 5(c1); the observation in Fig. 6(a) would contradict the physical intuition about inelastic fluids that a less viscous fluid should have a higher flow rate in the channel flow. After the flow reaches its steady state, the maximum flow velocity in the classical FENE-P model is $U=56.31~{}cm/s$ ($18.71$ in reduced unit) resulting in $Re=174.29$ and $We=9.34$, while the flow velocity in the macroscale-microscale coupled model is $U=16.63~{}cm/s$ ($5.52$ in reduced unit) resulting in $Re=51.47$ and $We=2.76$.

With straight path lines steady channel flow velocity profiles are determined solely by the viscosity function. Hence the fluids elasticity affects the flow only in so far as a change of relaxation time alters the viscosity function. Fig. 7 shows a parametric sensitivity analysis of the FENE-P model in a channel flow, affecting the velocity profile $u(y)$, the polymer normal stress profile $\tau^{p}_{xx}(y)$ and the polymer shear stress profile $\tau^{p}_{xy}(y)$. In Fig. 7(a1-a3), we decrease the polymer viscosity $\eta_{p}$ from $2.0$ to $1.0$ while keeping the relaxation time $\lambda$ constant. We observe in Fig. 7(a1) that the FENE-P fluid with a decreased viscosity has a higher flow rate, which is consistent with the physical intuition of viscous flows. However, when we decrease the relaxation time $\lambda$ from $1.0$ to $0.2$ while keeping the polymer viscosity $\eta_{p}$ constant in Fig. 7(b1-b3), we find that the FENE-P fluid with a decreased relaxation time has a lower velocity profile and a decreased flow rate, which is an opposite trend as the effects of decreasing the viscosity. Therefore, the flow profiles of the microscale-macroscale coupled simulation (FENE-P+DPD) shown in Fig. 6 are determined by the competition between viscous effects and elastic effects in the channel flow.

For a shear-thinning fluid flowing though a straight channel, a higher Weissenberg number can yield a more flattened velocity profile around the channel centerline and a sharper layer in the vicinity of the wall, which will increase the maximum strain rate in the system. For simplicity, we call the case of channel flow at $\tilde{R}=10$ and $\tilde{W}=1$ as channel #1, and another case of channel flow at $\tilde{R}=1$ and $\tilde{W}=5$ as channel #2. The GPR model for the effective viscosity $\eta$ and the first normal stress difference $N_{1}$ (computed for the relaxation time) constructed in the channel #1 can be transferred to the channel #2. As shown in Fig. 5(c2). the GPR model for $\eta$ and $N_{1}$ based on the seven training points in the channel #1 is valid up to $\dot{\gamma}_{\rm gpr}=3.93~{}{\rm s^{-1}}$. Fig. 8(b) plots the profiles of the shear strain rate $\dot{\gamma}$, the dynamic viscosity $\eta$ and relaxation time $\lambda$ in channel flow #2. The inset of Fig. 8(b) shows the time evolution of the global maximum strain rate of $\dot{\gamma}_{\max}$, where an overshoot response induced by the increased $\lambda$ is observed, leading to $\dot{\gamma}_{\max}=5.08~{}{\rm s^{-1}}$ beyond the point $\dot{\gamma}_{\rm gpr}=3.93~{}{\rm s^{-1}}$ in Fig. 5(c2). Then, we need to perform additional DPD simulations to include new data points once the acquisition function $\alpha$ breaks the acceptance criterion $\max(\alpha)<\delta_{\rm tol}$.

Fig. 9 presents the transfer-learning of $\eta$ and $N_{1}$ from channel #1 ($\tilde{R}=10,~{}\tilde{W}=1$) to channel #2 ($\tilde{R}=1,~{}\tilde{W}=5$). Using the same method for extrapolation in Fig. 5 with $\beta_{\eta}=6.41$ and $\beta_{N_{1}}=14.19$, the next sampling point determined by Eq. (10) for extrapolation is $\dot{\gamma}_{p_{5}}=4.97~{}{\rm s^{-1}}$, as shown in Fig. 9(a2). Subsequently, we perform a DPD simulation of a simple shear flow at $\dot{\gamma}_{p_{5}}=4.97~{}{\rm s^{-1}}$ and have a GPR model for $\eta$ and $N_{1}$ valid up to $\dot{\gamma}_{\rm gpr}=5.63~{}{\rm s^{-1}}$ as shown in Fig. 9(b2), which is sufficient to close the macroscopic SEM solver for channel #2 whose maximum strain rate is $\dot{\gamma}_{\max}=5.08~{}{\rm s^{-1}}$. After the flow reaches its steady state, the maximum flow velocity in the classical FENE-P model is $U=7.97~{}cm/s$ ($3.74$ in reduced unit) resulting in $Re=3.50$ and $We=9.31$, while the flow velocity in the macroscale-microscale coupled model is $U=2.24~{}cm/s$ ($1.05$ in reduced unit) resulting in $Re=0.98$ and $We=2.62$.

The active learning method is now applied to pressure-driven flow past a circular cylinder within a straight channel similar to the previous cases, the channel flow $\#1$ and $\#2$, with walls located at $y=\pm 2L_{0}$. The inlet is at $x=-10L_{0}$, the outlet is at $x=15L_{0}$. A circular cylinder with diameter $D=2L_{0}$ placed at the origin, as shown in Fig. 1. The flow is driven by a uniform pressure gradient, $f={dp}/{dx}$, along the channel. Boundary conditions are imposed on the walls as: the no-slip for velocity, and the Neumann condition for the extra-stress tensor. Periodic conditions are imposed in the streamwise direction. The cylinder Reynolds number is $Re=UD/\nu_{0}=2\tilde{R}U/(2.14\tilde{U})$ where the diameter $D=2L_{0}$, $U$ and $\nu_{0}=2.14\tilde{\nu}$ are the mean flow velocity and the zero-shear-rate viscosity, respectively. The Weissenberg number is defined as $We=\lambda_{0}U/D=0.5\tilde{W}U/\tilde{U}$ where $\lambda_{0}$ is the relaxation time. Flow past the cylinder is initiated by direct transfer of the constitutive equation as previously learned in the channel flow $\#2$.

The constitutive closure shown in Fig. 9(b) or Fig. 10(a) is valid only up to $\dot{\gamma}_{\rm gpr}=5.63~{}s^{-1}$ to satisfy $\max{\alpha}<\delta_{\rm tol}$. As the velocity field around the cylinder develops, the active-learning scheme automatically detects an acceptance violation since the GPR model was constructed for the channel flow $\#2$. This initializes new DPD simulations at extrapolation sampling points $\dot{\gamma}_{p6}=7.23~{}s^{-1}$ and $\dot{\gamma}_{p7}=9.34~{}s^{-1}$, which are determined by Eq. (10) and $\alpha_{\eta}(\dot{\gamma}_{p})=\beta_{\eta}\delta_{\rm tol}$ with $\beta_{\eta}=6.41$. Fig. 10(b1) shows how the GPR-informed method improves the closure of previous channel flow data $\#2$ by adding sampling points at $\dot{\gamma}_{p6}$ and $\dot{\gamma}_{p7}$, thus yielding an updated constitutive closure valid up to $\dot{\gamma}_{\rm gpr}=10.43s^{-1}$. The filled square in Fig. 10(b2) shows the global maximum shear strain rate to be $\dot{\gamma}_{\rm max}=8.52~{}s^{-1}<\dot{\gamma}_{\rm gpr}$ for cylinder flow at $\tilde{R}=1,~{}\tilde{W}=2$ and $\tilde{f}=3$. This example demonstrates how the learned effective closure constructed in the channel flow can be transferred directly to flow past a cylinder with very few additional expensive DPD simulations needed for a satisfactory closure.

With a constant uniform driving force $\tilde{f}=3$, the flow past a cylinder in the channel eventually reaches a steady state, where the maximum shear strain rate for the microscale-macroscale coupled system is $\dot{\gamma}_{\max}=7.12s^{-1}$. The effective viscosity of this shear-thinning fluid is $\nu=1.53\times 10^{-4}~{}m^{2}/s$ at shear strain rate of $7.12s^{-1}$. With the same driving force $\tilde{f}=3$, we simulate flow past a cylinder confined in a channel for two FENE-P fluids (FENE-P) and a microscale-macroscale coupled fluid (FENE-P+DPD). At their steady states, the maximum flow velocity at the inlet boundary is $U=3.54~{}cm/s$ for the classical FENE-P fluid resulting in $Re=2.45$ and $We=2.62$, $U=1.69~{}cm/s$ for the microscale-macroscale coupled fluid resulting in $Re=1.17$ and $We=1.25$. Fig. 11 compares stress component $\tau_{xx}^{p}$ for the FENE-P fluids and a macroscale-microscale coupled fluid along different sections of flow past a cylinder under the same driving force $\tilde{f}=3$. The three selected sections for extracting the stress component $\tau_{xx}^{p}$ are indicated by the blue solid lines in the insets of Fig. 11, i.e., section (a) is along the horizontal centerline $\left(y=0\right)$ and on the periphery of the cylinder, section (b) is along the vertical centerline $\left(x=0\right)$ and on the periphery of the cylinder, and section (c) is a vertical line $\left(x=1\right)$ tangent to the cylinder. Results show that the classical FENE-P fluid generates much higher stress $\tau_{xx}^{p}$ on the solid-wall surfaces than the macroscale-microscale coupled fluid. For a shear-thinning polymer solution, the relaxation time decreases as the shear strain rate increases. The classical FENE-P model that uses a constant relaxation time independent of flow conditions may overestimate the magnitude of stress in the flow past a cylinder. Because the shear strain rate can be large on the no-slip wall surfaces, we observe significant difference of stress on the wall surface and also on the cylinder surface between the classical FENE-P models and the macroscale-microscale coupled model, as shown in Fig. 11.

Upon increasing the Reynolds number by changing the flow conditions, the active-learning scheme is able to automatically detect the inaccuracy of the learned constitutive closure, and subsequently it initiates additional DPD simulations for the extra data needed to once again close the microscale-macroscale coupled system. Fig. 12 shows the process of constructing a new GPR-informed constitutive closure model when the flow is changed from $\tilde{R}=1,~{}\tilde{W}=2,~{}\tilde{f}=3$ to another flow past a cylinder at $\tilde{R}=20,~{}\tilde{W}=1,~{}\tilde{f}=0.4$ in the simulation, wherein Fig. 12(a1) plots the previously learned constitutive relations in the case $\tilde{R}=1,~{}\tilde{W}=2,~{}\tilde{f}=3$, and Fig. 12(b1) show an updated constitutive relations by adding one sampling point at $\dot{\gamma}_{p8}=12.92~{}s^{-1}$ according to Eq. (10). However, the learned constitutive relations are still insufficient to close the microscale-macroscale coupled system, hence the active-learning scheme adaptively initiates two more DPD simulations at sampling points $\dot{\gamma}_{p9}=16.56~{}s^{-1}$ and $\dot{\gamma}_{p10}=0.57~{}s^{-1}$ to construct a constitutive closure model valid up to $\dot{\gamma}_{\rm gpr}=18.42s^{-1}$, as shown in Fig. 12(d1). Therefore, three more training data are required to once again close the microscale-macroscale coupled system when we change to the flow past a cylinder with $\tilde{R}=20,~{}\tilde{W}=1,~{}\tilde{f}=0.4$.

The flow past a cylinder confined in the channel driven by a constant pressure gradient $\tilde{f}=0.4$ with $\tilde{R}=20$ and $\tilde{W}=1$ will eventually reach a steady state, when the maximum shear strain rate for the microscale-macroscale coupled system is $\dot{\gamma}_{\max}=15.01~{}s^{-1}$. Considering the shear-thinning effect, the effective viscosity of the fluid is $\nu=1.10\times 10^{-4}~{}m^{2}/s$ at the shear strain rate of $15.01~{}s^{-1}$. Using the same driving force $\tilde{f}=0.4$, we carried out simulations of flow past a cylinder in a channel for two classical FENE-P fluids (FENE-P) and a microscale-macroscale coupled fluid (FENE-P+DPD). After the flows reach their steady states, the maximum flow velocity at the inlet boundary is $U=16.10~{}cm/s$ for the classical FENE-P fluid resulting in $Re=70.45$ and $We=1.89$, $U=12.27~{}cm/s$ for the microscale-macroscale coupled fluid resulting in $Re=53.69$ and $We=1.44$.

Fig. 13 shows the streamline contours for (a, b) the classical FENE-P viscoelastic model and (c) the coupled macroscale-microscale fluid model for flow past a circular cylinder at $\tilde{R}=20,~{}\tilde{W}=1$ and $\tilde{f}=0.4$, where a stationary vortex pair behind each cylinder is shown in Figs. 13(a-c). The recirculation length $L_{w}$ is used to quantify the size of a steady separation bubble behind the cylinder. The values of $L_{w}$ in Figs. 13(a), 13(b) and 13(c) are $2.45D$, $1.20D$ and $0.76D$, respectively, which indicates that under the same driving pressure gradient the wake length in the macroscale-microscale coupled system generates a shorter wake than the classical FENE-P viscoelastic model. This is because the polymer solution possesses a shear-thinning behavior and has larger viscosity in the low shear strain rate regions. The contours of the effective viscosity $\eta$ in the flow past a circular cylinder in steady state are shown in Fig. 13(d), and the contours of the relaxation time $\lambda$ are shown in Fig. 13(e). We observe that both $\eta$ and $\lambda$ in the macroscale-microscale coupled system are significantly reduced in the vicinity of the cylinder surface and channel walls. Because the friction force on the channel wall surfaces can create an effective confinement suppressing flow separation, we observe that the macroscale-microscale coupled system produces a shorter wake behind the cylinder in Fig. 13 compared to the classical FENE-P models.

The drag per unit length on the cylinder is computed from [29],

$$F_{D}=\int_{0}^{2\pi}\left[\left(-p+2\eta_{s}\frac{\partial u}{\partial x}+\tau_{xx}\right)\cos\theta+\left(\eta_{s}\left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)+\tau_{xy}\right)\sin\theta\right]_{r=R}Rd\theta.$$

(11)

Fig. 14 shows the drag coefficient ${F_{D}}/{(1/2\rho UD^{2})}$ for the FENE-P fluids and the microscale-macroscale coupled fluid. The results show that under the same driving force the microscale-macroscale coupled fluid produces larger drag coefficient than the classical FENE-P models, which is consistent with our observations in Fig. 13 that the classical FENE-P models with constant viscosity and relaxation time produce wakes with longer recirculation length.

Fig. 15 compares the normal stress component $\tau_{xx}^{p}$ for the classical FENE-P fluids and the coupled macroscale-microscale fluid with $\tilde{R}=20,~{}\tilde{W}=1$ along different sections of flow past a cylinder driven by the same force $\tilde{f}=0.4$. The three selected sections for extracting the stress component $\tau_{xx}^{p}$ are indicated by the blue solid lines in the insets of Fig. 15, i.e., section (a) is along the horizontal centerline $\left(y=0\right)$ and on the periphery of the cylinder, section (b) is along the vertical centerline $\left(x=0\right)$ and on the periphery of the cylinder, and section (c) is a vertical line $\left(x=1\right)$ tangent to the cylinder. Similar to the previous case of flow past a cylinder, the classical FENE-P fluids generate higher surfaces stress $\tau_{xx}^{p}$ on the solid-wall than the macroscale-microscale coupled fluid. The decrease of relaxation time with shear-rate of the latter renders it less elastic than the FENE-P model with its constant relaxation time which overestimates stress levels in the flow past a cylinder. Fig. 15 shows significant differences of normal stress $\tau_{xx}^{p}$ on both the wall surfaces and on the cylinder surface between the standard FENE-P model and the macroscale-microscale coupled model.