Lagrangian Heterogeneous Multiscale Method (LHMM) for Simulating Polymer Solutions/Melts Behavior under Complex Flows using DPD-SPH

TL;DR

This work introduces LHMM, a fully Lagrangian multiscale framework that couples microscale DPD with FENE chains to a GENERIC-consistent SPH macroscale for simulating viscoelastic polymer melts in complex flows. Through bidirectional macro–micro coupling via Irving–Kirkwood averaging and a continuous intermittent time stepping scheme, the method preserves deformation history and yields invariant Weissenberg numbers across scales. Microscopic rheology is characterized with virtual rheometry and Carreau–Yasuda fits, providing input parameters such as zero-shear viscosity and relaxation time that drive macro predictions. Benchmark tests (Reverse Poiseuille Flow, periodic cylinder arrays) and flow in porous media demonstrate accurate transient and steady-state viscoelastic behavior, scalability to hundreds of millions of particles, and robustness across flow regimes, underscoring the method’s potential for predictive, constitutive-free multiscale modeling in polymeric systems.

Abstract

We present a Lagrangian Heterogeneous Multiscale Method (LHMM) for simulating the non-Newtonian rheology of polymer melts in complex two-dimensional flows. The method couples Dissipative Particle Dynamics (DPD) at the microscale with a GENERIC-compliant Smoothed Particle Hydrodynamics (SPH) at the macroscale, in a concurrent framework, overcoming the limitations of traditional Eulerian-based methods in capturing long-memory and history-dependent effects. At the microscale, DPD serves as a virtual rheometer, employing FENE (Finitely Extensible Nonlinear Elastic) bead-spring polymer chains. This approach provides key rheological properties, including shear-thinning and zero-shear-rate viscosities, relaxation times, and viscoelastic dynamics, which are quantified via Carreau-Yasuda fitting and spectral analysis. The LHMM couples SPH-derived strain rates with microscopic stress responses using the Irving-Kirkwood formalism. This approach enables a concurrent interaction between macroscopic strain rates and microscopic stress tensors, ensuring a consistent viscoelastic response across scales. The method is validated against benchmark flows, including Reverse Poiseuille Flow and flow through a Periodic Array of Cylinders, across Weissenberg numbers $0.5 < \text{Wi} < 30$ and low Reynolds numbers ($\text{Re} < 1$). A final demonstration of flow in a 2D porous medium highlights LHMM's capability to handle highly heterogeneous geometries. The LHMM is implemented in LAMMPS, making it suitable for integrating multiple models to describe microscales. In contrast, large-scale simulations efficiently utilize GPU and CPU resources, managing multiple coupling and time-scaling levels to maintain numerical stability and accuracy. The framework offers a predictive, constitutive-free tool that links microscopic polymer dynamics to macroscopic flow behavior, making it suitable for multiscale applications.

Figures (19)

• Figure 1: (a) Sketch of the Arbitrary Flow Boundary Conditions and their definition of buffer (), Boundary Condition (), and Core () regions for the / simulation. This diagram also illustrates how each micro-domain is delimited within the framework for monodisperse bead-spring chains (polymer melts). Here, is the total length of the arbitrary flow boundary setup, is the length of the core region, is the length of the boundary condition region, and is the length of the buffer region. (b) Details of the initial configuration within the core regions for polymer melts with different chain lengths, $\in \{8, 16, 32\}$.
• Figure 2: (a) Carreau-Yasuda model validation for polymer melts under steady shear flow. Lines represent fitted viscosity using the Carreau-Yasuda model, and symbols correspond to simulation data for various polymer chain lengths ($N_{\rm d}$). (b) Details of the distribution of bead-chains within the regions for polymer melts with different $N_{\rm d}=[8;16;32]$ after reaching the shear rate steady-state of $\dot{\gamma}_{xy}=0.02$. (c) Normal stress coefficients ($\psi_1$) for polymers under steady shear flow. Lines represent fitted values based on the Carreau-Yasuda model, and symbols indicate simulation data for different polymer chain lengths ($N_{\rm d}$). (d) Validation of relaxation time using constant shear rate for $= 8$ to $= 32$. The asymptotic values of $\lambda$ were determined as the shear rate ($\dot{\gamma}_{xy} \to 0$) approaches zero. (e) Relaxation Time ($\lambda$) as a function of polymer chain length (). The fitted curve shows a power-law relationship $\lambda \sim N_{\rm d}^{1.90}$, consistent with the Rouse model scaling Goto2021Murashima2011.
• Figure 3: Comparison of viscoelastic properties for $=16$ (a) and $=32$ (b). (I) Storage modulus ($G'$) and loss modulus ($G"$) derived from and , normalized by their crossover point value $\text{G}_\text{c}$ (=). (II) Complex viscosity from and (See Section \ref{['SM']} in the \ref{['SM']}), normalized by the maximum complex viscosity $\eta'_{\text{max}}$. $1/\lambda$ is the inverse of the characteristic relaxation time from the model, represented in Panels (I) as a vertical line, and CP is the crossover point.
• Figure 4: Micro–macro coupling in for $=1.0\cdot10^{-3}$ and $=1.0\cdot10^{-4}$ with $=16$ ($=0.46$, $=0.83$, and $=2.7$). (a) Spatial distribution of key variables at $t=\lambda(16)$. Sub-panel (a.I) shows the local Weissenberg number, $=(\underline{\underline{\dot{\boldsymbol{\gamma}}}}:\underline{\underline{\dot{\boldsymbol{\gamma}}}})^{1/2}\lambda$, and sub-panel (a.II) presents the dimensionless polymeric stress magnitude, $=$$($:$)^{1/2}/\tau_{\rm max}^{\rm RPF}$. Streamlines of five representative particles are overlaid in both panels to illustrate their trajectories. s(b) Micro-domain (polymer melt ) distributions at $t/\lambda(16)= 5.5$ for different chain lengths $=\{8,16,32\}$. Sub-panel (b.I) corresponds to the micro-domain linked to particle 5, while sub-panel (b.II) shows the one linked to particle 2. Both provide insight into the local polymeric structure, its dependence on chain length, and its deformation relative to the initial monomer distribution (Fig. \ref{['Fig:1']}(b)). Videos showing the time evolution of the micro-domain deformation for each chain length are provided in the \ref{['SM']}.
• Figure 5: Variation of $=\{8, 16, 32\}$ and their effects on the time evolution of the average velocity in the $x$-direction ($v_{x}$) and the shear stress of , compared with Newtonian fluid reference simulations (-only). (a) and (b) with $\alpha = 1.0\cdot10^{-4}$ and $\beta=\{0.15, 0.46, 0.55\}$ for $N_{\rm d}=\{8,16,32\}$, and (c) with $\alpha = 2.0\cdot10^{-3}$ and $\beta=\{0.77, 0.94, 0.96\}$ for $N_{\rm d}=\{8,16,32\}$. Results are shown for: (a) $F=1.8\cdot10^{-4}$ with ${\rm Wi}=\{0.13, 0.45, 1.52\}$. (b) $F=1.0\cdot10^{-3}$ with ${\rm Wi}=\{0.7, 2.7, 9.9\}$. (c) $F=1.0\cdot10^{-2}$ with ${\rm Wi}=\{4.15, 15.35, 53.43\}$. (I) Zoomed-in view of the velocity plateau, with vertical lines marking $\lambda(N_{\rm d})$. (II) Velocity response over time for different $N_{\rm d}$. (III) Shear stress response over time, showing the dimensionless total shear stress $=\bar{\pi}_{xy}/\tau_{\rm max}^{\rm RPF}$.