A high-resolution finite-difference method for simulating two-fluid, viscoelastic gel dynamics

Abstract

An important class of gels are those composed of a polymer network and fluid solvent. The mechanical and rheological properties of these two-fluid gels can change dramatically in response to temperature, stress, and chemical stimulus. Because of their adaptivity, these gels are important in many biological systems, e.g. gels make up the cytoplasm of cells and the mucus in the respiratory and digestive systems, and they are involved in the formation of blood clots. In this study we consider a mathematical model for gels that treats the network phase as a viscoelastic fluid with spatially and temporally varying material parameters and treats the solvent phase as a viscous Newtonian fluid. The dynamics are governed by a coupled system of time-dependent partial differential equations which consist of transport equations for the two phases, constitutive equations for the viscoelastic stresses, two coupled momentum equations for the velocity fields of the two fluids, and a volume-averaged incompressibility constraint. We present a numerical method based on a staggered grid, second order finite-difference discretization of the momentum equations and a high-resolution unsplit Godunov method for the transport equations. The momentum and incompressibility equations are solved in a coupled manner with the Generalized Minimum Residual (GMRES) method using a multigrid preconditioner based on box-relaxation. We present results on the accuracy and robustness of the method together with an illustration of the interesting behavior of this gel model for the four-roll mill problem.

Keywords: Krylov subspace; Mixture theory; Multigrid; Multiphase flow; Transient network model; Viscoelastic flow simulations.

Fig. 1.

Location of the unknowns in the space- and time-staggered grid for the 2-D gel model: ▹ = network/solvent horizontal velocity, Δ = network/solvent vertical velocity, •∙= pressure, ◯ = network/solvent volume fractions, □ = components of the viscoelastic stress tensor and $z$.

Fig. 2.

Test 1: Left column of (a), (b), and (c) displays the respective network volume fraction, network velocity, and solvent velocity at $t=4$ using the Set 1 parameters. Right column displays the relative errors in the solutions for the corresponding variable in the left column using various values of the grid spacing $h$.

Fig. 3.

Test 1: Continuation of Fig. 2, but for the (a) elastic modulus, (b) trace of the viscoelastic stress, and (c) viscoelastic shear stress.

Fig. 4. Test 2: Left column of (a), (b), and (c) displays the respective network volume fraction, network velocity, and solvent velocity at $t=4$ using the Set 2 parameters. Right column displays the relative errors in the solutions for the corresponding variable in the left column using various grid spacing $h$.

Fig. 5.

Test 2: Continuation of Fig. 4, but for the (a) elastic modulus, (b) trace of the viscoelastic stress, and (c) viscoelastic shear stress.

Fig. 6.

Initial volume fraction (57) for the third refinement test.

Fig. 7.

Test 3: Left column of (a), (b), and (c) displays the respective network volume fraction, network velocity, and solvent velocity at $t=2.5$ using the Set 3 parameters. Right column displays the relative errors in the solutions for the corresponding variable in the left column using various grid spacing $h$.

Fig. 8.

Test 3: Continuation of Fig. 7, but for the (a) elastic modulus, (b) trace of the viscoelastic stress, and (c) viscoelastic shear stress. Note that only the region $-0.4⩽x,y⩽0.4$ is displayed for the figures in the left column to better illustrate the structure of the variables near the concentration of network.

Fig. 9.

Stills of ${\theta }_{\mathrm{n}},{\mathbf{u}}_{\mathrm{n}}$, and ${\mathbf{u}}_{\mathrm{s}}$ from the first snapback test using the Set 3 parameters listed in Table 3. Red (dark) curves in velocity field plots correspond to the contour ${\theta }_{\mathrm{n}}(x,y,t)=0.1$, while the yellow (light) curves correspond to the initial contour ${\theta }_{\mathrm{n}}(x,y,0)=0.1$. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10.

Stills of ${\theta }_{\mathrm{n}},{\mathbf{u}}_{\mathrm{n}}$, and ${\mathbf{u}}_{\mathrm{s}}$ from the second snapback test using the Set 4 parameters listed in Table 3. Red (dark) curves in velocity field plots correspond to the contour ${\theta }_{\mathrm{n}}(x,y,t)=0.1$, while the yellow (light) curves correspond to the initial contour ${\theta }_{\mathrm{n}}(x,y,0)=0.1$. (For interpretation of the references in colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 11.

Time traces of the number of iterations required for the iterative solver to reach a tolerance of 10−8) for the two snapback tests. The jump in iteration count from 3 to 6 in both tests corresponds to time $t=2.5$ when the background force is turned off.