A thermodynamically consistent Johnson-Segalman-Giesekus model: numerical simulation of the rod climbing effect

TL;DR

This work addresses rod climbing in viscoelastic fluids by deriving a thermodynamically consistent Johnson--Segalman--Giesekus (JSG) model within the Rajagopal--Srinivasa framework and implementing it via a high-order ALE finite-element method in Firedrake. It contrasts the thermodynamically consistent Model I with the engineering Model II, showing that only Model I guarantees nonnegative dissipation and yields physically plausible free-surface behavior that agrees with experimental data, especially Beavers and Joseph's rod-climbing measurements. The authors provide a robust numerical pipeline, including a monolithic formulation, high-order elements, and open-source code, enabling accurate predictions of climbing heights and surface shapes in axisymmetric rotating-rod configurations. The results highlight the importance of thermodynamic consistency in constitutive modeling for free-surface viscoelastic flows and establish a foundation for future three-dimensional extensions and parameter-matching with experiments.

Abstract

Viscoelastic rate-type fluids represent a popular class of non-Newtonian fluid models due to their ability to describe phenomena such as stress relaxation, non-linear creep, and normal stress differences. The presence of normal stress differences in a simple shear flow gives rise to forces acting in directions orthogonal to the primary flow direction. The rod climbing effect, i.e. the rise of a fluid along a rod rotating about its axis, is associated with this phenomenon. Within the class of viscoelastic rate-type fluids that includes the Oldroyd-B and Giesekus models with Gordon--Schowalter convected derivatives, we show -- by means of thermodynamical analysis and numerical simulations -- that a thermodynamically consistent variant of the Johnson--Segalman model captures experimental data exceedingly well and is therefore superior to other models in this class, including the standard Johnson--Segalman model, which is widely used in engineering applications but is shown here to be incompatible with the second law of thermodynamics. We release a robust and computationally efficient higher-order finite-element implementation as open-source software on GitHub. The implementation is based on an arbitrary Lagrangian--Eulerian (ALE) formulation of the governing equations and is developed using the Firedrake library.

Figures (9)

• Figure 1: Sketch of the reference configuration $\kappa_R(\mathcal{B})$, the current configuration $\kappa_t(\mathcal{B})$ and the natural configuration $\kappa_{p(t)}(\mathcal{B})$. The deformation gradient $\mathbb{F}_{\kappa_R}$ is multiplicatively decomposed. The natural configuration is defined as the configuration that the body in the current configuration would take if the external stimuli were removed. Hence the natural configuration $\kappa_{p(t)}(\mathcal{B})$ is associated with the current configuration $\kappa_t(\mathcal{B})$ and it evolves with the body as the body produces entropy. This allows us to split the total deformation $\mathbb{F}_{\kappa_R}$ into a purely elastic (reversible) part $\mathbb{F}_{\kappa_{p(t)}} : \kappa_{p(t)}(\mathcal{B}) \longrightarrow \kappa_{t}(\mathcal{B})$ and the rest (dissipation) $\mathbb{G} : \kappa_{R}(\mathcal{B}) \longrightarrow \kappa_{p(t)}(\mathcal{B})$ so that: $\mathbb{F}_{\kappa_R} = \mathbb{F}_{\kappa_{p(t)}} \mathbb{G}$.
• Figure 2: Sketch of the axi-symmetric rod climbing configuration used in all simulations. The problem is formulated in the meridional plane of cylindrical coordinates, with a rotating rod, a deformable free surface, and no-slip container walls.
• Figure 3: Spatial discretization convergence is examined for $\omega = 2.6$ rev/s. The meshes are labelled M0, …, M4 from coarsest to finest. Mesh M1 with $p=4$ (approximately $1.3\times 10^{5}$ DoF) produces results essentially identical in the region near the rod to those obtained on mesh M4 with $p=2$ (approximately $1.5\times 10^{6}$ DoF), demonstrating that higher-order discretizations achieve same accuracy at much lower computational cost. The convergence is not strictly monotone, with small oscillations around what appears to be the mesh-converged solution.
• Figure 4: Comparison of the Oldroyd-B results obtained in the present work with the reference results FIGUEIREDO201698. The observed differences may partly be attributed to spatial under-resolution in the reference simulations, as suggested by the coarse structure of the reported climbing vortex.
• Figure 5: Visualization of the resulting steady state for $\omega = 2.9$ rev/s in ParaView. The domain is colored and warped by the computed mesh displacement $\hat{\mathbf{u}} = (0,0,\hat{u}_z)$ and shown in a three-dimensional perspective obtained by a $270^\circ$ rotation about the rotation axis. The right face cut shows the computational mesh (M1), while the left face cut displays a line integral convolution (LIC) representation of the velocity field $\mathbf{v} = (v_r,0,v_z)$ in the meridional plane, showing the secondary flow. The rotating rod, shown for reference, is not a part of the simulation; its motion is prescribed via boundary conditions.

Theorems & Definitions (1)

• Remark : Remarks on the Johnson--Segalman Model II in the rheological literature