A Three-Dimensional Open-Source Solver for Incompressible Viscoelastic Two-Component Flows Open Access

Viscoelastic fluids play a vital role in numerous practical applications, including paints and coatings [1], food and cosmetics [2], oil and gas drilling [3], microfluidics [4], vibration dampening [5], and biomedical applications [6]. These fluids exhibit both fluid-like and solid-like behaviors. However, in the realm of simulation, the numerical community has encountered the challenge known as the high Weissenberg number problem (HWNP), which affects viscoelastic fluid flow. This issue arises when the numerical solution diverges beyond a certain value of the nondimensional Weissenberg number. In response to the HWNP, various numerical methods, computational techniques, and analytical strategies have been developed over the years. To provide a comprehensive overview of this area of research, we will conduct an extensive literature review of the previous work related to this topic.

Oliveira et al. [7] made significant contributions by developing a collocated finite volume method on unstructured grids for simulating nonlinear elastic flows. Their investigations encompassed the channel entry flow problem, as well as flow around bounded and unbounded cylinders. In another notable study, Phillips and Williams [8] explored viscoelastic flow through a planar contraction using a staggered grid-based finite volume method. Their discretization approach for the convection term in the Navier–Stokes equations and the constitutive equations was conducted in a semi-Lagrangian manner, revealing interesting insights into the effect of the Weissenberg number on vortex behavior during contraction.

Additionally, Aboubacar et al. [9] introduced innovative hybrid cell-vertex and cell-centered finite volume approaches for analyzing viscoelastic flows. They employed a semi-implicit formulation using the Taylor–Galerkin method, complemented by a pressure-correction step in their hybrid case, while the cell-centered formulation included a semi-Lagrangian step for convection terms in both the momentum and constitutive equations.

Tome et al. [10] successfully developed a numerical technique for simulating the free surface flow of Generalized Newtonian fluids. Their approach utilized a finite difference method solver on staggered grids, incorporating marker particles to effectively capture the free surface dynamics. In addressing the challenges associated with high Weissenberg numbers, Fattal and Kupferman [11] proposed a matrix logarithm of the conformation tensor for various viscoelastic constitutive models, demonstrating its utility by simulating the lid-driven cavity problem at a Weissenberg number of five. They further expanded their techniques by developing a log conformation method (LCM) using a finite difference method formulation, which proved effective in mitigating high Weissenberg number problems.

Afonso et al. [12] conducted a thorough investigation into the log conformation method's performance with a finite volume-based formulation for steady and unsteady creeping flows of viscoelastic fluids, finding it to be robust and stable even at high Deborah numbers. Trebotich et al. [13] presented an impressive stable and convergent hyperbolic scheme for simulating viscoelastic flows in channels with sudden contractions, enabling the simulation of a diverse range of flows from highly viscous to highly elastic.

Employing a discontinuous Galerkin method-based finite element formulation alongside discrete elastic viscous stress splitting (DEVSS) and the log conformation approach, Hulsen et al. [14] effectively simulated viscoelastic flow around a cylinder, achieving converged results even at elevated Weissenberg numbers. Hao and Pan [15] also utilized the log conformation approach in conjunction with finite element methods to study viscoelastic flow within a lid-driven cavity, attaining converged results at high Weissenberg numbers.

Further explorations by Oliveira et al. [16] focused on the influence of contraction ratios on the flow behavior of Oldroyd-B and Phan–Thien–Tanner fluids in axisymmetric contractions, analyzing vortex patterns across various contraction ratios and Deborah numbers. Kane et al. [17] contributed by comparing the performance of four distinct finite element method implementations of the log conformation approach in the context of viscoelastic flow around a cylinder, highlighting differences in their treatment of advection terms in the governing equations.

In the realm of the extrudate swell problem, Mompean et al. [18] introduced an algebraic model for elastic stresses, noting a significant reduction in computational demands with this approach compared to more traditional differential equations. Habla et al. [19] took strides in developing a viscoelastic two-phase flow solver using the open-source framework openfoam, implementing a range of differential constitutive equations and validating their findings against literature data.

Lastly, Balci et al. [20] proposed a square root-based symmetric factorization of the conformation tensor, demonstrating its efficiency and stability at high Weissenberg numbers. Oishi et al. [21] developed a finite difference method-based viscoelastic free surface solver employing the extended Pom–Pom model as the constitutive equation, effectively investigating the die swell effect and achieving converged results for Weissenberg numbers up to 20. Based on a previous formulation, Figueiredo et al. [22] extended a two-dimensional viscoelastic flow solver for the extended Pom-Pom model to three dimensions, investigating both the three-dimensional extrudate swell and jet buckling problems.

Following a comprehensive review of the existing literature, it is evident that the log-conformation approach stands out as a robust technique for addressing the challenges associated with high Weissenberg numbers. Additionally, it is noteworthy that the majority of research to date has utilized two-dimensional formulations developed in-house. In this paper, we introduce a three-dimensional implementation of the log-conformation approach, leveraging the open-source flow solver, basilisk. This solver extends the two-dimensional framework established by López-Herrera et al. [23]. The paper is structured as follows: the forthcoming section outlines the governing equations pertinent to our study; the subsequent section details the numerical formulation; we then provide an overview of the results from various test cases; and finally, we offer concluding remarks on our findings.

In this section, the governing equations for incompressible, isothermal, single, and two-phase viscoelastic flows are presented.

In this context, $μs$ represents the viscosity of the solvent, while D denotes the rate of deformation tensor. Additionally, $∇u$ refers to the velocity gradient tensor, and $∇uT$ is the transpose of this gradient tensor.

The positive definiteness of conformation tensor $(A)$ is essential for ensuring the stability of the solution according to Eq. (10). However, in high elasticity flows, the conformation tensor $(A)$ may not meet the criteria for strict positive definiteness (SPD), potentially leading to numerical instability. In regions characterized by high deformation rates, the stretching and relaxation terms can exhibit exponential growth. The convection term plays a crucial role in balancing this growth; however, since it relies on polynomial interpolations, it may struggle to effectively counteract the exponential amplification. This leads to what is known as the HWNP. When the Weissenberg number is high, the numerical solution may experience significant growth over time, marked by an exponential increase in deformation and relaxation terms.

where $ΛΨ$ and $ΛA$ are tensors formed by eigenvalues of log conformation tensor and conformation tensor, respectively. The Log conformation method helps the conformation tensor to be strictly positive-definite. The positive-definiteness can be verified by checking $det(A)≥0$⁠. Hence, the high Weissenberg number instability is removed.

We developed the three-dimensional, viscoelastic two-phase solver based on the open-source code basilisk [24] which provides several libraries and solvers for compressible, incompressible, multiphase flows. The incompressible flow solver uses a second-order accurate projection method. The steps involved in the time stepping of the Navier–Stokes equation is described as follows:

Step 3: Following this, we proceed to advance the viscoelastic stresses to the midstep $n+12$⁠, resulting in $τpn+12$⁠.

where the subscripts cc and fc, respectively, denote the cell-centered and the face centered quantities. $f→c$ refers to the cell-centered averaging process of face-centered values.

Here $c→f$ denotes the face-centered averaging process of cell-centered quantities.

Further details on the discretization process can be found in Popinet [27] and Afkhami and Bussmann [28]. Here, the time-step, $△t$ is calculated based on the Courant–Friedrich–Levy condition.

Assuming polymeric stresses are available at midstep $n−12$⁠, $τpn−12$⁠, then time stepping advances as follows:

Step 2: The conformation tensor is diagonalized, $A=RΛRT$⁠, to obtain its eigenvalues and eigenvectors, $Λn−12$ and $Rn−12$⁠.

Step 4: The velocity field is calculated at time-step n, $un$⁠. Thus, velocity gradient is decomposed according to Eq. (11) to calculate $Bn$ and $Ωn$⁠.

to obtain $Λ**$⁠, $R**$ and $A**=RΛRT|**$⁠.

Further details on the discretization process can be found in the literature, check, for instance, Refs. [23] and [30]. The eigenvalues and eigenvectors of the various matrices are calculated using the routines provided by Kopp [31,32].

In this section, we perform various tests for the numerical validation of the viscoelastic single and two-phase solvers.

In this test case, we explore the behavior of a viscoelastic (Oldroyd-B) fluid within a three-dimensional lid-driven cavity shown in Fig. 1. The upper boundary of the cavity is designed to move with a tangential velocity that varies with both time and position, as described by the equation: $ux(x,t)=8[1+tanh(8t−4)]x2(1−x)2$⁠. The other walls of the cavity remain stationary, adhering to a no-slip boundary condition for velocity. For this numerical simulation, we have selected specific parameters: a Reynolds number of $Re=0.01$⁠, a solvent viscosity ratio of $β=0.5$⁠, and a Weissenberg number of $Wi=1$⁠. To achieve accurate results, a uniform grid as detailed in Table 1 was utilized. The maximum time-step for the M1 grid is set at $△t=5×10−4$⁠, while for the M2 grid, it is $△t=5×10−5$⁠.

We present two velocity profiles: the horizontal velocity $ux$ as a function of the y coordinate at $x=0.5$⁠, and the vertical velocity $uy$ as a function of the x coordinate at $y=0.5$⁠. In Figs. 2 and 3, we illustrate the horizontal velocity (⁠$ux$⁠) and vertical velocity (⁠$uy$⁠) profiles on two distinct grids, M1 (dashed line) and M2 (solid line), alongside the results from Habla et al. [33] (dotted line). The results indicate a commendable convergence of the velocity profiles with those reported by Habla et al. [33], suggesting that the grid size has minimal impact on both $ux$ and $uy$⁠.

The components of the polymeric stresses (⁠$τp$⁠) in the $x−y$ plane using the M2 grid are depicted in Figs. 4(a)–4(f). The stress contours for $τp,xx$⁠, $τp,xy$⁠, and $τp,yy$ demonstrate a noteworthy agreement with the findings of Habla [33]. Specifically, $τp,xx$ reveals a boundary layer at the top wall of the lid, while the contours for $τp,xy$ and $τp,yy$ are more pronounced in the upper right corner. Additionally, the magnitudes of $τp,xz$ and $τp,yz$ are relatively smaller compared to the other components, with both developing near the upper left and upper right corners of the lid. The behavior of $τp,zz$ varies from $−1.6$ to 1.6, as illustrated in Fig. 4(f). It is noteworthy that the contour lines are more pronounced in the upper half of the lid, with minimal stress evident in the bottom half.

The streamlines and contours of velocity magnitude in the driven cavity at $Wi=1$ are illustrated in Figs. 5(a) and 5(b). Notably, a primary vortex is present, accompanied by a smaller secondary vortex in the bottom left corner, as observed in Fig. 5(a). Additionally, Fig. 5(b) displays the contours of velocity magnitude at $Wi=1$⁠, demonstrating findings that are consistent with those reported by Habla [33].

In this section, we evaluate the performance of our two-phase solver using the three-dimensional drop impact scenario. As depicted in Fig. 6, we consider the initial configuration of a viscoelastic drop with a diameter D, which is released from a height of $1.5×D$⁠. The fluid surrounding the drop within the computational domain is a Newtonian fluid. In this test, we focus on analyzing the time-dependent evolution of the drop's shape as it falls under the influence of gravity. Our findings are compared with those documented in the literature by Figueiredo et al. [34], providing a comprehensive perspective on the problem.

In this test case, we consider the initial velocity of the drop as $uy=−1 m/s$ with diameter $D=0.02m$⁠. The distance from the center of the drop to bottom wall is 0.04 m. The gravity magnitude is $g=−9.81 m/s2$⁠. We have used Oldroyd-B model in our implementation and transformed all the dimensional data into nondimensional numbers as follows: $Re=5.0$⁠, $Fr=2.26$⁠, $Wi=1.0$ and $β=0.1$⁠. No slip boundary condition were applied on the bottom wall of the domain. Numerical simulation has been implemented with two kind of grid size M1 and M2 with a maximum time-step of $△t=1.0×10−4$⁠. Figure 7 represents the numerical prediction of the variation of the dimensionless width, $w=W/D$⁠, of the viscoelastic drop versus the dimensionless time, $t=t*U/D$⁠, on M1 and M2 grid. These results are compared with those given in Ref. [34]. The results obtained for the two grids are consistent and in good agreement with the data from literature.

The behavior of the viscoelastic fluid after the impact on surface depends on the interplay between surface tension and polymeric stresses. Figure 8 shows the isometric view of 3D droplet spreading obtained on M2 grid for $Wi=1$ and $β=0.1$⁠. The figure represents the evolution in the shape of droplet at different time instants, denoted as $t=0$⁠, $t=1.7$⁠, $t=2.0$⁠, $t=2.2$⁠, $t=3.5$⁠, $t=5.0$⁠. At $t=0$⁠, the droplet is in perfect spherical shape before the impact with the surface. This shape is the reference point for measuring the deformation and spreading over time. At $t=1.7$⁠, the droplet contacts with the surface and starts to spread out. As a result, the height of the droplet decreases while its base area increases. Between $t=2.0−2.2$⁠, the droplet continues to spread out and it become more flattened with a broader base. At $t=3.5−5.0$⁠, the droplet reaches to the maximum base area and tries to regain its initial state till its shape looks like a cookie.

The polymeric stress components $τp,xx$⁠, $τp,xy$⁠, $τp,xz$⁠, $τp,yy$⁠, $τp,yz$⁠, and $τp,zz$ after impact of the drop on the bottom wall are shown in Figs. 9(a)–9(f). The stresses in these figures are represented on $x−y$⁠, $y−z$⁠, and $z−x$ planes and the white line represents the interface between drop and surrounding fluid.

In this section, we aim to validate our formulation through a specific test case involving a 3D viscoelastic drop submerged in a Newtonian fluid, which deforms under a pure Couette shear flow. In Fig. 10, we present the initial configuration employed for the numerical simulation of the viscoelastic drop in response to the pure Couette shear flow. The properties of the fluids are characterized by their respective density and viscosity, denoted as $ρ1$⁠, $μ1$ for the viscoelastic fluid and $ρ2$⁠, $μ2$ for the surrounding Newtonian fluid. The radius of the viscoelastic drop is represented as $Rd$⁠, with the overall domain size set to $16Rd×8Rd$⁠. The interfacial tension between the two fluids is indicated by $γ$⁠. Under the influence of the Couette shear flow, the top and bottom walls are set into motion in opposite directions, with velocities $VT$ and $VB$⁠, resulting in an effective constant shear rate for the droplet, calculated as $γ˙=(VT−VB)/H$⁠. For this investigation, we utilize two types of uniform grids: $M1(64×64×64)$ and $M2(128×128×128)$⁠.

The results obtained in this study were compared with those presented by Figueiredo et al. [34]. The nondimensional parameters considered include $Re=0.3$⁠, $We=0.18$⁠, $μr=ρr=1$⁠, $De=0.4$⁠, and $β=0.5$⁠. In Fig. 11, we illustrate the time evolution of the deformation parameter (⁠$Φ$⁠) related to the viscoelastic drop deformation within a Newtonian fluid. The deformation parameter $Φ$ is defined as the ratio of the difference between the maximum and minimum distances of the droplet interface from its center, expressed as $Φ=(Rmax−Rmin)/(Rmax+Rmin)$⁠. The findings demonstrate a commendable alignment with existing literature data. Furthermore, Fig. 12 presents the interface shape of the deformed droplet at $t=10$⁠, showcasing the three-dimensional form of the droplet in comparison to the two-dimensional results reported by Figueiredo [34].

In this context, we consider $ϵ=0.2$ as the perturbation parameter, which serves to initiate the capillary thinning of the filament. The parameter $R0$ represents the initial radius of the filament. Several nondimensional parameters are utilized in the simulation, including the Deborah number (⁠$De=λρR03/γ$⁠), Ohnsorge number (⁠$Oh=μ0ρR0γ$⁠), and viscosity ratio $β=μsμ0$⁠, and are detailed in Table 2. A higher value of the Deborah number, specifically $(De=60)$⁠, suggests that the elastic properties of the fluid are becoming increasingly significant compared to its viscous characteristics.

Figure 13 illustrates the thinning process of the Oldroyd-B filament within a three-dimensional solution domain defined by $0≤x,y,z≤2π$⁠, represented by a blue square. The images capturing the filament's thinning are derived through the application of symmetry and reflection across the x-axis at various nondimensional time instants: (a) t = 0, (b) t = 26.04, (c) t = 30.38, (d) t = 34.72, and (e) t = 42.42. In the initial phases, from (a) to (c), the morphological changes in the filament occur gradually, as viscous forces exert a greater influence during this interval until t = 26.04. Subsequently, during phases (c) and (e), the thinning process advances further, with surface tension forces (or capillary forces) becoming increasingly significant in driving the dynamics. These forces act to minimize the surface area of the thread, leading to a notable decrease in its diameter over time. This phenomenon is recognized as capillary-driven thinning.

The evolution of the viscoelastic filament interface is illustrated in Fig. 14(a) for the time interval between $t=26.04$ and $t=42.42$⁠. During this period, it can be observed that the neck region connecting the slender thread to the drop becomes increasingly steep as time progresses. Additionally, Figs. 14(b) and 14(c), respectively, present the changes in drop diameter (D) and thread diameter (d) at various moments. Notably, the thinning of the filament follows a faster trend after $t=26.04$⁠, indicating a transition as the filament evolves into a drop linked by a slender thread. This process contributes to an increase in the drop radius and a reduction in the thread radius over time.

Understanding the axial component of velocity is crucial for comprehending the viscoelastic filament thinning process. In Fig. 15(a), we observe the contours of the dimensionless axial velocity component, $ux$⁠, at time t = 34.72. It is important to note that capillary forces play a significant role in the thinning process by facilitating the movement of fluid from the thinning filament to the drop, where the capillary pressure is comparatively lower. This pressure gradient effectively drives fluid flow from the center of the thinning thread toward the drop, resulting in an axial flux directed toward the drop. This flux contributes to the growth of the drop while simultaneously promoting the thinning of the connected thread.

In Fig. 15(b), we present a plot illustrating the axial velocity component through the center of the filament over different times. Notably, the axial velocity experiences an increase in magnitude near the drop-filament junction, as indicated by the two parallel black dotted lines, before gradually decreasing along the remaining length of the thread. The plot demonstrates that, as time progresses, the axial velocity continues to evolve while maintaining an almost linear profile in the filament, signifying a consistent thinning of the filament.

In Fig. 16(a), we observe the contour map representing the nondimensional axial polymeric stress component, $τp,xx$⁠, at time t = 34.72. This stress component acts in the axial direction and appears minimal within the drop, while showing considerable significance within the filament. It plays a crucial role in resisting thread collapse, which is induced by capillary forces. Examination of two cross-sectional planes at $x=3.1$ and $x=5$ indicates that $τp,xx$ tends to increase along the length of the filament.

Furthermore, Fig. 16(b) illustrates the variation of $τp,xx$ along the centerline of the filament over time. Notably, while there is little variation of $τp,xx$ within the drop, a significant rise is observed at the junction of the drop and filament, with more gradual changes occurring along the length of the filament. It is important to highlight that $τp,xx$ along the centerline shows a steady increase as time progresses. The distribution of axial stress $τp,xx$ as a function of radius of the filament at different nondimensional times is shown in Fig. 17. It can be noticed in the figure that the radius of the filament decreases while the stress values increase as time progresses. Also, the location of peak stress shifts toward the center of the filament.

In Fig. 18(a), we observe the contour map distribution of the polymeric stress component, $τp,yy$ at t = 34.72. Unlike $τp,xx$⁠, $τp,yy$ demonstrates a more uniform distribution within the thread region; however, it exhibits increased stress near the drop, attributed to the elevated extensional strain rates. This stress component plays a crucial role in resisting the radial deformation of the thread as fluid is expelled toward the drop. The cross-sectional plane at $x=0.75$ reveals a heightened level of $τp,yy$ near the center of the drop, while the plane at $x=3.1$ illustrates an increasing trend of $τp,yy$ along the neck region. Figure 18(b) presents the variation of $τp,yy$ along the filament centerline over different time intervals. It is noteworthy that $τp,yy$ rises within the drop region but begins to decline starting from $x=0.75$⁠, ultimately reaching a stable value beyond the neck region.

Another component of polymeric stress, $τp,zz$⁠, is presented in Fig. 19(a) at t = 34.72. Much like $τp,yy$⁠, $τp,zz$ reflects the radial deformation of the filament, specifically within the $x−y$ plane. Notably, we observed higher values of $τp,zz$ within the drop, indicating a region of significant extensional deformation. Conversely, lower values are noted in the thread, attributable to comparatively lesser extensional deformation. Figure 19(b) illustrates the variation of $τp,zz$ along the filament centerline at various time points. As expected, the results show an increase in $τp,zz$ within the drop and a decrease along the thread. In the next section, we conclude this work.

In this work, we embarked on an innovative journey to develop, test, and validate a three-dimensional flow solver for viscoelastic two-phase flows. Built on the open-source flow solver basilisk, our creation embodies both robustness and accuracy. We applied the solver to standard single and two-phase test cases from the literature, tackling challenges like the three-dimensional driven cavity problem, the drop fall problem, the drop under shear problem, and the bead on a string problem. The results achieved by our solver across these diverse test cases match well with findings from previous literature, showcasing the powerful potential of our advancements. The role of polymeric stresses is highly significant in the process of viscoelastic filament thinning. In particular, the axial stress component, $τp,xx$⁠, plays a dominant role in the thread region due to the axial stretching as the fluid moves toward the drop from the filament. Additionally, the stresses $τp,yy$ and $τp,zz$ arise from the radial expansion of the drop and are notably pronounced in the drop region.

Thanks go to Professor Stéphane Popinet for developing the open-source software basilisk which is a wonderful piece of scientific software. Also, the original, two-dimensional implementation of Log conformation approach in basilisk by Professor JM Lopez Herrera is gratefully acknowledged.

A =

conformation tensor

B =

symmetric tensor

D =

rate of deformation tensor

D =

diameter

g =

gravity vector

h =

initial filament radius

I =

identity tensor

N =

antisymmetric tensor

p =

pressure

R =

orthogonal tensor

R =

radius

r =

thread radius

t =

time

u =

velocity vector

W =

width

$n̂$ =

unit normal vector

$α$ =

volume fraction

$β$ =

ratio of solvent viscosity with total viscosity

$γ$ =

surface tension

$γ˙$ =

shear rate

$δs$ =

Dirac delta function

$ϵ$ =

perturbation parameter

$θ$ =

any fluid property

$κ$ =

curvature

$λ$ =

relaxation time

$Λ$ =

diagonal matrix of eigenvalues

$μ$ =

viscosity

$μp$ =

polymer dynamic viscosity

$μr$ =

viscosity ratio

$μs$ =

solvent dynamic viscosity

$ρ$ =

density

$ρr$ =

density ratio

$τm$ =

total shear stress

$τs$ =

solvent stress

$τp$ =

polymeric stress

$Φ$ =

deformation parameter

$Ψ$ =

logarithm conformation tensor

$Ω$ =

antisymmetric tensors

De =

Deborah number

Fr =

Froude number

Oh =

Ohnsorge number

Re =

Reynolds number

Wi =

Weissenberg number

p =

polymer

s =

solvent

$n$ =

time instant

$▽$ =

upper convected time derivative

DEVSS =

discrete elastic viscous stress splitting

HWNP =

high Weissenberg number problem

LCM =

log conformation method