Abstract

In this study, we unveil a three-dimensional flow solver designed to simulate viscoelastic two-phase flows using the Oldroyd-B formulation. Acknowledging the challenges that researchers encounter in this dynamic field, we have integrated the three-dimensional Log conformation approach into the open-source flow solver basilisk, significantly enhancing its capabilities beyond its two-dimensional predecessors. Our solver stands as a testament to rigorous testing against a wide range of three-dimensional viscoelastic flow challenges, encompassing both single and two-phase scenarios drawn from established literature. True to its two-dimensional roots, it exhibits extraordinary robustness, adeptly managing viscoelastic flows, even at high Weissenberg numbers. By offering this powerful solver as an open-source resource, we aspire to empower the computational fluid dynamics community. We believe it will become an invaluable tool for researchers delving into the complexities of viscoelastic flows, fostering innovation and inspiring new progress in the field.

Introduction

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.

Governing Equations

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

Continuity and Navier–Stokes Equations.

The mass and momentum conservation equations are as follows:
(1)
(2)
For the case of two-phase flows, the advection equation for the interface is given by
(3)
In this context, ρ, u(u,v,w), p, and γ refer to density, velocity, pressure, and surface tension, respectively. The term κ denotes the curvature of the free surface, while n̂ represents the unit normal vector to the interface. Additionally, δs is the Dirac delta function, and α signifies the volume fraction of the viscoelastic phase. The variable α is stored at grid cell centers and the velocity vectors used while solving the advection equation are stored at grid cell faces. The total shear stress in the viscoelastic fluid, denoted as τm, comprises both the solvent stress (τs), and the polymeric stress (τp)
(4)
The stress tensor in the solvent can be described by the Newtonian flow equation
(5)

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.

Constitutive Model for the Viscoelastic Fluid.

For this study, the Oldroyd-B model is used with constant polymeric viscosity which means there is no shear thinning during the flow of viscoelastic fluid. For the Oldroyd-B model, the polymeric stress τp is expressed as follows [23]:
(6)
where λ and μp are the relaxation time and polymeric viscosity, respectively. λ is associated with the Weissenberg number (Wi), which is a dimensionless number given by Wi=λU/l. Here, U and l are, respectively, the velocity scale and characteristic length scale in the problem. τp is the Upper convected time derivative of polymeric stresses and is given by
(7)
By using equation Eqs. 6 and (7), we get
(8)
In Eq. (8), the second and third terms on the left-hand side imply convection and deformation , respectively. The right-hand side terms denote the relaxation and the source term. The boundary condition used for the normal component of polymeric stresses at a wall is Dirichlet boundary condition with zero value. For shear components of the polymeric stresses, we used homogeneous Neumann boundary condition. In the Oldroyd-B model, polymeric stress tensor, τp can also be expressed in terms of conformation tensor A as
(9)
where I is the identity matrix. Using Eqs. (8) and (9), the constitutive equation in terms of the conformation tensor A becomes
(10)

Log Conformation Method.

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.

In their work, Fattal and Kupferman [11] proposed a strategy for addressing the HWNP by transforming the conformation tensor (A) into a logarithmic conformation tensor (Ψ). This approach aims to stabilize the numerical simulations. Their recommendation is to utilize the LCM as a means of tackling the challenges posed by high Weissenberg numbers. In the context of applying the LCM reformulation, the transpose of the velocity gradient can be effectively decomposed in a constructive manner
(11)
In this context, we consider the divergence-free velocity field denoted as u. The tensors N and Ω are characterized as antisymmetric, while B is recognized as a symmetrical, traceless tensor. Within this framework, we progress in time using the logarithm of the conformation tensor, expressed as Ψ=logA. It is important to note that Ψ represents the log conformation tensor, and A is a symmetric positive definite (SPD) matrix that can be diagonalized accordingly
(12)
where R is an orthogonal tensor which is formed by the three eigenvectors of A. Λ is a diagonal tensor which is formed by the three eigenvalues of A. Now, N, Ω, and B are decomposed by the tensor R and its transpose RT
(13)
(14)
(15)
The velocity gradient (u) can be derived using equations Eqs. (13)(15). In a similar manner, the Log-conformation constitutive equation is obtained by solving equations Eqs. (10)(12). By following the methodology proposed by Fattal and Kupferman [11], we arrive at the following equation:
(16)
The conformation tensor (A) and log-conformation tensor (Ψ) are linked with each other by the eigenvectors of the conformation tensor (A). Moreover, log-conformation tensor Ψ has homogeneous Neumann boundary conditions by default. We have
(17)
(18)

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.

Numerical Implementation

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 1: To advance the volume fraction α in time, we employ the nondiffusive, geometric Volume-of-Fluid method, represented by the following equation:
(19)
Step 2: Next, we update the fluid properties using the equation
(20)

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

Step 4: The “guessed” velocity u* is then calculated according to the following equation:
(21)
In basilisk, the above equation is solved by decomposing it into advection, diffusion, and acceleration steps. In the advection step, the equation subpart being solved is as follows:
(22)
This step is solved with the help of Bell–Collela–Glaz second-order accurate scheme, Bell et al. [25]. In the diffusion step, the following equation is solved with the help of multigrid methods, Popinet [26]
(23)
In the acceleration step, the effect of forcing term F due to surface tension, viscoelastic stresses, electric stresses etc. is applied. The velocity field under different forcing terms is computed as
(24)

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

Step 5: The last step is the projection step. The intermediate velocities ucc*** are not divergence free. Imposing ·un+1=0 yields a Poisson equation for the pressure
(25)
where ufc*** is evaluated as
(26)

Here cf denotes the face-centered averaging process of cell-centered quantities.

The Poisson equation for pressure is solved for pccn+1 using multigrid-based solvers. Finally, the cell-centered and the face-centered velocities are updated using
(27)
and
(28)

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.

The polymeric stress tensor (τp) is calculated by a time-split procedure [29]. This procedure decomposes Eq. (16) as follows:
(29)
(30)
(31)

Assuming polymeric stresses are available at midstep n12, τpn12, then time stepping advances as follows:

Step 1: Conformation tensor at midstep n12 is calculated from Eq. (9)
(32)

Step 2: The conformation tensor is diagonalized, A=RΛRT, to obtain its eigenvalues and eigenvectors, Λn12 and Rn12.

Step 3: The logarithm of the conformation tensor is calculated
(33)

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.

Step 5: The BCG scheme is used for the advection of log-conformation tensor
(34)
Step 6: Equation (30) is solved explicitly
(35)
Step 7: Log conformation tensor is diagonalized
(36)

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

Step 8: Converting Eq. (31) in terms of conformation tensor A and later integrated analytically,
(37)
Step 9: Then, An+12 is calculated by
(38)
Step 10: Finally
(39)

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].

Validations

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

Lid Cavity Flow.

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(8t4)]x2(1x)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×104, while for the M2 grid, it is t=5×105.

Fig. 1
Schematic of viscoelastic fluid flow in a three-dimensional lid driven cavity
Fig. 1
Schematic of viscoelastic fluid flow in a three-dimensional lid driven cavity
Close modal
Table 1

Grid refinement levels

RefinementLevel 6 (M1)Level 7 (M2)
Grid size64×64×64128×128×128
RefinementLevel 6 (M1)Level 7 (M2)
Grid size64×64×64128×128×128

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.

Fig. 2
Horizontal velocity (ux) profile at Wi=1.0 for different grid sizes
Fig. 2
Horizontal velocity (ux) profile at Wi=1.0 for different grid sizes
Close modal
Fig. 3
Vertical velocity (uy) profile at Wi=1.0 for different grid sizes
Fig. 3
Vertical velocity (uy) profile at Wi=1.0 for different grid sizes
Close modal

The components of the polymeric stresses (τp) in the xy 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.

Fig. 4
Polymeric stresses τp in x−y plane on M2 grid at Wi=1.0 in Lid driven cavity
Fig. 4
Polymeric stresses τp in x−y plane on M2 grid at Wi=1.0 in Lid driven cavity
Close modal

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].

Fig. 5
Streamlines and Velocity magnitude contours on M2 grid at Wi=1.0 in Lid driven cavity
Fig. 5
Streamlines and Velocity magnitude contours on M2 grid at Wi=1.0 in Lid driven cavity
Close modal

Drop Impact Case.

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.

Fig. 6
Three-dimensional sketch of the computational domain for the drop impact case
Fig. 6
Three-dimensional sketch of the computational domain for the drop impact case
Close modal

In this test case, we consider the initial velocity of the drop as uy=1m/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.81m/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×104. 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.

Fig. 7
Time variation of the width of the viscoelastic drop on different grids
Fig. 7
Time variation of the width of the viscoelastic drop on different grids
Close modal

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.02.2, the droplet continues to spread out and it become more flattened with a broader base. At t=3.55.0, the droplet reaches to the maximum base area and tries to regain its initial state till its shape looks like a cookie.

Fig. 8
Isometric view of the three-dimensional droplet spreading obtained on M2 grid for Wi=1 and β=0.1
Fig. 8
Isometric view of the three-dimensional droplet spreading obtained on M2 grid for Wi=1 and β=0.1
Close modal

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 xy, yz, and zx planes and the white line represents the interface between drop and surrounding fluid.

Fig. 9
Contours of polymeric stresses (τp,xx, τp,xy, τp,xz, τp,yy, τp,yz, and τp,zz) after the drop impact on the bottom wall
Fig. 9
Contours of polymeric stresses (τp,xx, τp,xy, τp,xz, τp,yy, τp,yz, and τp,zz) after the drop impact on the bottom wall
Close modal

Drop Deformation Due to Shear.

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 γ˙=(VTVB)/H. For this investigation, we utilize two types of uniform grids: M1(64×64×64) and M2(128×128×128).

Fig. 10
Initial configuration of a three-dimensional viscoelastic drop under shear deformation
Fig. 10
Initial configuration of a three-dimensional viscoelastic drop under shear deformation
Close modal

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 Φ=(RmaxRmin)/(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].

Fig. 11
Time evolution for the deformation parameter (Φ) for the viscoelastic drop deforming under shear
Fig. 11
Time evolution for the deformation parameter (Φ) for the viscoelastic drop deforming under shear
Close modal
Fig. 12
Deformation of the droplet under shear at a nondimensional time of t=10
Fig. 12
Deformation of the droplet under shear at a nondimensional time of t=10
Close modal

Application: Viscoelastic Filament Thinning

In this study, we present an investigation involving 3D numerical simulation of a viscoelastic filament undergoing thinning within a domain defined by 0x,y,z2π. The initial configuration of the viscoelastic filament takes the form of a cylinder with a surface perturbation of ϵ=0.2 and is surrounded by an air environment. The initial cylindrical shape, oriented parallel to the x-axis, is characterized by its radius as follows:
(40)

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.

Table 2

The nondimensional parameters for BOAS

DeOhβρgρsμgμsϵ
603.160.250.010.010.2
DeOhβρgρsμgμsϵ
603.160.250.010.010.2

Figure 13 illustrates the thinning process of the Oldroyd-B filament within a three-dimensional solution domain defined by 0x,y,z2π, 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.

Fig. 13
Thinning of a viscoelastic filament under capillary forces. The actual solution domain is the three-dimensional geometry shown with the blue rectangle.
Fig. 13
Thinning of a viscoelastic filament under capillary forces. The actual solution domain is the three-dimensional geometry shown with the blue rectangle.
Close modal

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.

Fig. 14
(a) Evolution of the filament interface between t=26.04 and t=42.42, (b) variation of the drop diameter (D) with nondimensional time, (c) variation of the thread diameter, and (d) with nondimensional time
Fig. 14
(a) Evolution of the filament interface between t=26.04 and t=42.42, (b) variation of the drop diameter (D) with nondimensional time, (c) variation of the thread diameter, and (d) with nondimensional time
Close modal

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.

Fig. 15
(a) Contours of the dimensionless axial velocity, ux, at t=34.72 and (b) plot of axial velocity through the center of the filament at different times
Fig. 15
(a) Contours of the dimensionless axial velocity, ux, at t=34.72 and (b) plot of axial velocity through the center of the filament at different times
Close modal

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.

Polymeric Stresses During Filament Thinning.

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.

Fig. 16
(a) Contour plot of the axial polymeric stress component, τp,xx, at t=34.72 and (b) plot of τp,xx as a function of the axial location passing through the center of the filament at different nondimensional times
Fig. 16
(a) Contour plot of the axial polymeric stress component, τp,xx, at t=34.72 and (b) plot of τp,xx as a function of the axial location passing through the center of the filament at different nondimensional times
Close modal

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.

Fig. 17
The distribution of axial stress, τp,xx, as a function of radius of the filament at different nondimensional times. The stress values are from midaxial location, x=3.14.
Fig. 17
The distribution of axial stress, τp,xx, as a function of radius of the filament at different nondimensional times. The stress values are from midaxial location, x=3.14.
Close modal

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.

Fig. 18
(a) Contour plot of the polymeric stress component, τp,yy, at t=34.72 and (b) variation of τp,yy along the axial location passing through the center of the filament at different nondimensional times
Fig. 18
(a) Contour plot of the polymeric stress component, τp,yy, at t=34.72 and (b) variation of τp,yy along the axial location passing through the center of the filament at different nondimensional times
Close modal

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 xy 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.

Fig. 19
(a) Contour plot of the polymeric stress component, τp,zz, at t=34.72 and (b) a plot of τp,zz as a function of the axial location passing through the center of the filament and at different nondimensional times
Fig. 19
(a) Contour plot of the polymeric stress component, τp,zz, at t=34.72 and (b) a plot of τp,zz as a function of the axial location passing through the center of the filament and at different nondimensional times
Close modal

Conclusions

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.

Acknowledgment

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.

Nomenclature

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

References

1.
Carré
,
A.
,
Gastel
,
J.-C.
, and
Shanahan
,
M. E.
,
1996
, “
Viscoelastic Effects in the Spreading of Liquids
,”
Nature
,
379
(
6564
), pp.
432
434
.10.1038/379432a0
2.
Brummer
,
R.
,
2006
, “
Skin and its Care
,”
Rheology Essentials of Cosmetic and Food Emulsions
,
Springer
,
Berlin, Heidelberg
, pp.
15
16
.10.1007/3-540-29087-7_3
3.
Ofei
,
T. N.
,
Lund
,
B.
,
Saasen
,
A.
, and
Sangesland
,
S.
,
2022
, “
The Effect of Oil–Water Ratio on Rheological Properties and Sag Stability of Oil-Based Drilling Fluids
,”
ASME J. Energy Resour. Technol.
,
144
(
7
), p.
073008
.10.1115/1.4052033
4.
Tian
,
F.
,
Feng
,
Q.
,
Chen
,
Q.
,
Liu
,
C.
,
Li
,
T.
, and
Sun
,
J.
,
2019
, “
Manipulation of Bio-Micro/Nanoparticles in non-Newtonian Microflows
,”
Microfluid. Nanofluid.
,
23
(
5
), pp.
1
9
.10.1007/s10404-019-2232-z
5.
Makris
,
N.
,
Constantinou
,
M.
, and
Dargush
,
G.
,
1993
, “
Analytical Model of Viscoelastic Fluid Dampers
,”
J. Struct. Eng.
,
119
(
11
), pp.
3310
3325
.10.1061/(ASCE)0733-9445(1993)119:11(3310)
6.
Misra
,
J.
,
Sinha
,
A.
, and
Shit
,
G.
,
2010
, “
Flow of a Biomagnetic Viscoelastic Fluid: Application to Estimation of Blood Flow in Arteries During Electromagnetic Hyperthermia, a Therapeutic Procedure for Cancer Treatment
,”
Appl. Math. Mech.
,
31
(
11
), pp.
1405
1420
.10.1007/s10483-010-1371-6
7.
Oliveira
,
P. J.
,
Pinho
,
F. D.
, and
Pinto
,
G.
,
1998
, “
Numerical Simulation of Non-Linear Elastic Flows With a General Collocated Finite Volume Method
,”
J. Non-Newtonian Fluid Mech.
,
79
(
1
), pp.
1
43
.10.1016/S0377-0257(98)00082-2
8.
Phillips
,
T.
, and
Williams
,
A.
,
1999
, “
Viscoelastic Flow Through a Planar Contraction Using a semi-Lagrangian Finite Volume Method
,”
J. Non-Newtonian Fluid Mech.
,
87
(
2–3
), pp.
215
246
.10.1016/S0377-0257(99)00065-8
9.
Aboubacar
,
M.
,
Phillips
,
T. N.
,
Tamaddon-Jahromi
,
H.
,
Snigerev
,
B.
, and
Webster
,
M.
,
2004
, “
High-Order Finite Volume Methods for Viscoelastic Flow Problems
,”
J. Comput. Phys.
,
199
(
1
), pp.
16
40
.10.1016/j.jcp.2004.01.033
10.
Tomé
,
M. F.
,
Grossi
,
L.
,
Castelo
,
A.
,
Cuminato
,
J. A.
,
Mangiavacchi
,
N.
,
Ferreira
,
V. G.
,
De Sousa
,
F.
, and
McKee
,
S.
,
2004
, “
A Numerical Method for Solving Three-Dimensional Generalized Newtonian Free Surface Flows
,”
J. Non-Newtonian Fluid Mech.
,
123
(
2–3
), pp.
85
103
.10.1016/j.jnnfm.2004.06.012
11.
Fattal
,
R.
, and
Kupferman
,
R.
,
2004
, “
Constitutive Laws for the Matrix-Logarithm of the Conformation Tensor
,”
J. Non-Newtonian Fluid Mech.
,
123
(
2–3
), pp.
281
285
.10.1016/j.jnnfm.2004.08.008
12.
Afonso
,
A.
,
Oliveira
,
P. J.
,
Pinho
,
F.
, and
Alves
,
M.
,
2009
, “
The Log-Conformation Tensor Approach in the Finite Volume Method Framework
,”
J. Non-Newtonian Fluid Mech.
,
157
(
1–2
), pp.
55
65
.10.1016/j.jnnfm.2008.09.007
13.
Trebotich
,
D.
,
Colella
,
P.
, and
Miller
,
G.
,
2005
, “
A Stable and Convergent Scheme for Viscoelastic Flow in Contraction Channels
,”
J. Comput. Phys.
,
205
(
1
), pp.
315
342
.10.1016/j.jcp.2004.11.007
14.
Hulsen
,
M. A.
,
Fattal
,
R.
, and
Kupferman
,
R.
,
2005
, “
Flow of Viscoelastic Fluids Past a Cylinder at High Weissenberg Number: Stabilized Simulations Using Matrix Logarithms
,”
J. Non-Newtonian Fluid Mech.
,
127
(
1
), pp.
27
39
.10.1016/j.jnnfm.2005.01.002
15.
Hao
,
J.
, and
Pan
,
T.-W.
,
2007
, “
Simulation for High Weissenberg Number: Viscoelastic Flow by a Finite Element Method
,”
Appl. Math. Lett.
,
20
(
9
), pp.
988
993
.10.1016/j.aml.2006.12.003
16.
Oliveira
,
M. S.
,
Oliveira
,
P. J.
,
Pinho
,
F. T.
, and
Alves
,
M. A.
,
2007
, “
Effect of Contraction Ratio Upon Viscoelastic Flow in Contractions: The Axisymmetric Case
,”
J. Non-Newtonian Fluid Mech.
,
147
(
1–2
), pp.
92
108
.10.1016/j.jnnfm.2007.07.009
17.
Kane
,
A.
,
Guénette
,
R.
, and
Fortin
,
A.
,
2009
, “
A Comparison of Four Implementations of the Log-Conformation Formulation for Viscoelastic Fluid Flows
,”
J. Non-Newtonian Fluid Mech.
,
164
(
1–3
), pp.
45
50
.10.1016/j.jnnfm.2009.08.003
18.
Mompean
,
G.
,
Thais
,
L.
,
Tomé
,
M. F.
, and
Castelo
,
A.
,
2011
, “
Numerical Prediction of Three-Dimensional Time-Dependent Viscoelastic Extrudate Swell Using Differential and Algebraic Models
,”
Comput. Fluids
,
44
(
1
), pp.
68
78
.10.1016/j.compfluid.2010.12.010
19.
Habla
,
F.
,
Marschall
,
H.
,
Hinrichsen
,
O.
,
Dietsche
,
L.
,
Jasak
,
H.
, and
Favero
,
J. L.
,
2011
, “
Numerical Simulation of Viscoelastic Two-Phase Flows Using openFOAM®
,”
Chem. Eng. Sci.
,
66
(
22
), pp.
5487
5496
.10.1016/j.ces.2011.06.076
20.
Balci
,
N.
,
Thomases
,
B.
,
Renardy
,
M.
, and
Doering
,
C. R.
,
2011
, “
Symmetric Factorization of the Conformation Tensor in Viscoelastic Fluid Models
,”
J. Non-Newtonian Fluid Mech.
,
166
(
11
), pp.
546
553
.10.1016/j.jnnfm.2011.02.008
21.
Oishi
,
C. M.
,
Martins
,
F.
,
Tome
,
M. F.
,
Cuminato
,
J. A.
, and
Mckee
,
S.
,
2011
, “
Numerical Solution of the eXtended Pom-Pom Model for Viscoelastic Free Surface Flows
,”
J. Non-Newtonian Fluid Mech.
,
166
(
3–4
), pp.
165
179
.10.1016/j.jnnfm.2010.11.001
22.
Figueiredo
,
R.
,
Oishi
,
C.
,
Cuminato
,
J. A.
, and
Alves
,
M. A.
,
2013
, “
Three-Dimensional Transient Complex Free Surface Flows: Numerical Simulation of XPP Fluid
,”
J. Non-Newtonian Fluid Mech.
,
195
, pp.
88
98
.10.1016/j.jnnfm.2013.01.004
23.
López-Herrera
,
J.-M.
,
Popinet
,
S.
, and
Castrejón-Pita
,
A.-A.
,
2019
, “
An Adaptive Solver for Viscoelastic Incompressible Two-Phase Problems Applied to the Study of the Splashing of Weakly Viscoelastic Droplets
,”
J. Non-Newtonian Fluid Mech.
,
264
, pp.
144
158
.10.1016/j.jnnfm.2018.10.012
24.
Popinet
,
S.
,
2020
, “
Basilisk Flow Solver and PDE Library
,” accessed Apr. 9, 2025, http://basilisk.fr/
25.
Bell
,
J. B.
,
Colella
,
P.
, and
Glaz
,
H. M.
,
1989
, “
A Second-Order Projection Method for the Incompressible Navier-Stokes Equations
,”
J. Comput. Phys.
,
85
(
2
), pp.
257
283
.10.1016/0021-9991(89)90151-4
26.
Popinet
,
S.
,
2015
, “
A Quadtree-Adaptive Multigrid Solver for the Serre–Green–Naghdi Equations
,”
J. Comput. Phys.
,
302
, pp.
336
358
.10.1016/j.jcp.2015.09.009
27.
Popinet
,
S.
,
2003
, “
Gerris: A Tree-Based Adaptive Solver for the Incompressible Euler Equations in Complex Geometries
,”
J. Comput. Phys.
,
190
(
2
), pp.
572
600
.10.1016/S0021-9991(03)00298-5
28.
Afkhami
,
S.
, and
Bussmann
,
M.
,
2008
, “
Height Functions for Applying Contact Angles to 2D VOF Simulations
,”
Int. J. Numer. Methods Fluids
,
57
(
4
), pp.
453
472
.10.1002/fld.1651
29.
Fattal
,
R.
, and
Kupferman
,
R.
,
2005
, “
Time-Dependent Simulation of Viscoelastic Flows at High Weissenberg Number Using the Log-Conformation Representation
,”
J. Non-Newtonian Fluid Mech.
,
126
(
1
), pp.
23
37
.10.1016/j.jnnfm.2004.12.003
30.
Popinet
,
S.
,
2009
, “
An Accurate Adaptive Solver for Surface-Tension-Driven Interfacial Flows
,”
J. Comput. Phys.
,
228
(
16
), pp.
5838
5866
.10.1016/j.jcp.2009.04.042
31.
Kopp
,
J.
,
2008
, “
Efficient Numerical Diagonalization of Hermitian 3× 3 Matrices
,”
Int. J. Mod. Phys. C
,
19
(
3
), pp.
523
548
.10.1142/S0129183108012303
32.
Max Planck Society,
2024
, “
Numerical Diagonalization of 3x3 Matrices
,”
Max Planck Society
,
Munich, Germany
, accessed May 12, 2025, https://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html.
33.
Habla
,
F.
,
Tan
,
M. W.
,
Haßlberger
,
J.
, and
Hinrichsen
,
O.
,
2014
, “
Numerical Simulation of the Viscoelastic Flow in a Three-Dimensional Lid-Driven Cavity Using the Log-Conformation Reformulation in OpenFOAM®
,”
J. Non-Newtonian Fluid Mech.
,
212
, pp.
47
62
.10.1016/j.jnnfm.2014.08.005
34.
Figueiredo
,
R.
,
Oishi
,
C. M.
,
Afonso
,
A.
,
Tasso
,
I.
, and
Cuminato
,
J. A.
,
2016
, “
A Two-Phase Solver for Complex Fluids: Studies of the Weissenberg Effect
,”
Int. J. Multiphase Flow
,
84
, pp.
98
115
.10.1016/j.ijmultiphaseflow.2016.04.014