Numerical Simulation of Three-Dimensional Free Surface Flows Using the K–BKZ–PSM Integral Constitutive Equation

This work introduces a novel numerical method designed to address three-dimensional unsteady free surface flows incorporating integral viscoelastic constitutive equations, specifically the K–BKZ–PSM (Kaye–Bernstein, Kearsley, Zapas–Papanastasiou, Scriven, Macosko) model. The new proposed methodology employs a second-order finite difference approach along with the deformation fields method to solve the integral constitutive equation and the marker particle method (known as marker-and-cell) to accurately capture the evolution of the fluid’s free surface. The newly developed numerical method has proven its effectiveness in handling complex fluid flow scenarios, including confined flows and extrudate swell simulations of Boger fluids. Furthermore, a new semi-analytical solution for velocity and stress fields is derived, considering fully developed flows of a K–BKZ–PSM fluid in a pipe.


Introduction
Discovered and developed in the late twentieth century, viscoelastic materials have been used in a number of different applications (polymer industry, biomedicine, automotive industry, food, paints, etc.).Their use is often based on trial and error procedures, resulting in wasting raw material and time before a good end product is achieved.To mitigate this problem, numerical simulations are often used to predict the best material processing conditions.Usually, the simulations are based on the finite element, finite volume and finite difference methods, and the constitutive equations are, in most cases, defined by rheological differential models, such as Oldroyd-B [1], UCM (Upper-convected Maxwell) [2,3], PTT (Phan-Thien-Tanner) [4,5], FENE-P (Finite Extensible Nonlinear Elastic Peterlin) [6,7] and Giesekus models [7].Simulations of three-dimensional real-world applications require a great deal of computational effort, making the convergence of the algorithms in a reasonable amount of time a difficult task [4].However, recent technological advances in scientific computing and the development of faster computers have led researchers to perform simulations in more complex geometries and use more sophisticated rheological models (that use integral equations instead of partial differential equations).
It is known that integral constitutive equations can better model various viscoelastic fluids, such as high-density polyethylene (HDPE) [8,9] and low-density polyethylene (LDPE) [10] (used in the injection molding industry), and one of the most successful integral models is the K-BKZ-PSM [11][12][13] (see also [14,15]).Therefore, there is significant interest among research groups worldwide in developing numerical methods to deal with the K-BKZ model, particularly with an emphasis on its application to polymer flows.Many studies have focused on simulating data and phenomena associated with polymer melt flows in rheology and polymer processing; however, there is still a need for further efforts to tackle numerical solutions of the K-BKZ-PSM for three-dimensional, time-dependent, free surface flows.
The vast majority of problems studied in the literature (considering integral models) are about confined flows, such as entry flows [9,16,17] and flows in abrupt contractions [18][19][20][21][22]. Regarding free surface flows, Mitsoulis and Malamataris [20] extended the implementation of the free boundary condition (FBC) method to viscoelastic fluids governed by integral constitutive equations.Specifically, they focused on the K-BKZ-PSM model.To validate their numerical approach, they used the finite element method (FEM) to obtain results in simple test cases, including planar flow at an angle and Poiseuille flow in a tube, where analytical solutions are available for comparison.Furthermore, they have applied the FBC method to the K-BKZ-PSM model using data from a benchmark polymer melt, specifically the IUPAC-LDPE melt.Some other researchers have also considered flows with free surfaces [8,14,[23][24][25][26][27].Ganvir et al. [25] developed a novel approach for simulating extrudate swell using an Arbitrary Lagrangian Eulerian (ALE) technique in conjunction with a finite element formulation.The constitutive behavior of the melt was modeled using a differential exponential Phan-Thien-Tanner (PTT) viscoelastic model.With the proposed method, they have conducted simulations of extrudate swell in both planar and axisymmetric extrusion scenarios, which involve an abrupt contraction ahead of the die exit.Regarding three-dimensional (3D) flows, Rasmussen [28] developed a Galerkin finite element method for simulating three-dimensional transient viscoelastic flows.The method used a Lagrangian kinematic description and integral constitutive models.The numerical implementation was validated with the calculation of various transient and steady drag correction factors for the motion of a sphere in a cylinder containing an upper convected Maxwell fluid.Later, Marín and Rasmussen [29] extended the Galerkin finite element method for simulating three-dimensional transient and non-isothermal flows of K-BKZ type fluids.Tomé et al. [27] proposed a novel numerical approach to tackle the simulation of 3D viscoelastic unsteady free surface flows governed by the Oldroyd-B differential constitutive equation.The numerical method involves solving the governing equations using a finite difference approach on a 3D-staggered grid.To validate the accuracy and reliability of the proposed technique, an exact solution of the flow of an Oldroyd-B fluid inside a 3D-pipe was employed.The results obtained through numerical simulations included the analysis of transient extrudate swell and jet buckling.
Summarizing, previous studies in the field of free surface flows have predominantly centered around two-dimensional (2D) scenarios and used the finite element method.These investigations primarily revolve around the extrudate swell problem, considering both steady and unsteady flows; however, it is worth noting that these studies relied on differential viscoelastic constitutive equations.Therefore, this work introduces a novel numerical method specifically designed to address 3D unsteady free surface flows incorporating integral viscoelastic constitutive equations, specifically, the K-BKZ-PSM model.The key innovation lies in the development of a robust numerical method for integral models using the finite difference method on a staggered grid, which enables accurate predictions of extrudate swell phenomena.We also derive a new semi-analytical solution for the fully developed flow of a K-BKZ-PSM viscoelastic fluid, which can serve for other authors to verify their own numerical implementations of the K-BKZ-PSM integral viscoelastic model.
It is worth noting that the FEM prominently features in the limited body of work concerning this subject.Nevertheless, both FEM and FDM stand as extensively used numerical approaches for tackling partial differential equations (PDEs).When properly employed within suitable conditions, both techniques exhibit stability.Within our research group, a longstanding tradition exists regarding leveraging the finite difference method [2,4,27,[30][31][32], resulting in a profound mastery of its implementation.Furthermore, the group has made innovative strides in enhancing the fundamental finite difference methodology.This progression equips us to adeptly handle different grid structures and a range of discretization choices.Consequently, this method takes precedence in our current work.
The paper is structured as follows.In Section 2, we introduce the governing equations for isothermal and incompressible viscoelastic flows modelled by the K-BKZ-PSM constitutive equation.Section 3 is devoted to the numerical method, where the variant of the marker particle method that employs the finite difference method on a staggered grid is described for 3D flows using the K-BKZ-PSM viscoelastic integral model.In Section 4, we derive a new semi-analytical solution for the fully developed flow of a K-BKZ-PSM viscoelastic fluid.For validation of the newly developed numerical method, two case studies are analyzed in Section 5, the confined pipe flow and the extrudate swell free surface flow of a Boger fluid.The paper ends with the conclusions in Section 6.

Governing Equations
The isothermal and incompressible fluid flow considered in this work is governed by the dimensionless continuity and linear momentum equations [27], together with a constitutive equation for the stress.Φ is a stress tensor given by where τ is a non-Newtonian stress tensor, v(u, v, w) is the velocity field, p is the kinematic pressure, g is the gravity acceleration vector and t is the time.In these equations, Fr = U/ Lg is the Froude number, Re = ρ 0 UL/η 0 is the Reynolds number, η 0 is the zero-shear-rate viscosity, ρ 0 is the fluid density, g the magnitude of the gravity acceleration vector and U and L are the characteristic velocity and length scales, respectively.Note that all variables are dimensionless, with: x = x/L, v = v/U, t = t U/L, p = p/(ρU 2 ) and Φ = Φ/(ρU 2 ).
The constitutive equation for the non-Newtonian stress tensor is given by the K-BKZ-PSM model [11], where is the memory function, λ k is the relaxation time of the fluid, a k is a model parameter and m 1 is the number of modes.H(I 1 , I 2 ) is the Papanastasiou-Scriven-Macosko decay function, being given by B t (t) is the Finger tensor, and ] are the first and second invariants of B t (t), respectively.The parameters a k , λ k , α and β are obtained from a fit to rheological data.Wi = λ re f U/L is the Weissenberg number, the viscosity is given by η is the mean relaxation time [14].
In this work, the method of deformation fields [33] is used to update the Finger tensor as the fluid flows.In this methodology, (N + 1)-deformation instants (t ) are defined in the interval [0, t] where the history of deformation is stored.This deformation is updated by solving the transport equation for B t (t), The governing equations are solved in a Cartesian 3D system (x, y, z, t) where p = p(x, y, z, t), v = (u(x, y, z, t), v(x, y, z, t), w(x, y, z, t)) T , This results in the following system of equations that need to be solved for the pressure, velocity and stress: linear momentum equations: where g x , g y , g z are the Cartesian components of the gravity vector.

Numerical Method
The governing equations are solved by a variant of the marker particle method [4,27], which employs the finite difference method on a staggered grid.This methodology is implemented in the FREEFLOW-3D code developed by researchers from the Institute of Mathematical and Computing Sciences (ICMC) at the University of São Paulo (USP) in Brazil.Code details can be found in [4,27,30,31] considering 2D, 3D and radial symmetry flows.The precision of the numerical technique and its validation for three-dimensional viscoelastic flows with a free surface is presented in the works of Tomé et al. [27,30] (which only uses differential constitutive equations).The use of integral models in three-dimensional flows (considering free surface problems) has not yet been tested on this system (see Tomé et al. [4,27,30,31]).The novelty of this work is to incorporate equations for viscoelastic fluids using integral models (more complex than the differential type models) in the FREEFLOW-3D system.
In this methodology, the velocity field is approximated in the face of the cells, and the other variables, denoted by ζ, are evaluated in the center of the computational cells (see Figure 1a).The technique adopted here is presented by Tomé et al. [30,31] for differential models, where the free surface is defined by marker particles that move with the local fluid velocity.In addition, the computational cells are defined as (see also Figure 1b):  To solve Equations ( 1) and ( 2), one must specify boundary conditions for the velocity field.A velocity field (V in f ) is prescribed in the fluid inlet cells (inflows), and a fully developed flow is assumed in the outflow (a homogeneous Neumann boundary condition ∂v/∂n = 0 is assumed, where n is the normal direction to the contour).In the fluid inlet cells (inflows), we assume that ∂p ∂n = 0 and the Finger tensor B t is the identity matrix.In outflow regions, we assume Neumann conditions for Finger tensor ∂B ∂n = 0 , and we assume p = 0. We also take v = 0 in rigid boundaries.Details of the boundary conditions adopted in this work can be found in Tomé et al. [30] or Castelo et al. [31].
The solutions v(x, t n+1 ), p(x, t n+1 ) and τ(x, t n+1 ) at time step t n+1 = t + ∆t are obtained in the following way: first, using the values of τ(x, t n ), the velocity and pressure fields at time t n+1 are calculated.Then, v(x, t n+1 ) is used to calculate the tensor τ(x, t n+1 ) by the method of deformation fields, and, lastly, the free surface is updated.Specifically, the following steps are performed: Step 1 -Calculation of v(x, t n+1 ) and p(x, t n+1 ) It is assumed that, at time t, the variables v(x, t) = v (n) , p(x, t) = p (n) , τ(x, t) = τ (n) and the marker's positions x(t) = x (n) are known.Then, v(x, t n+1 ) and p(x, t n+1 ) are obtained as follows: , and, from the EVSS [34] transformation, Calculate an intermediate velocity field v (n+1) using the ideas of the projection method [30,31] to uncouple the conservation of mass and momentum equations.An intermediate velocity field v (n+1) is obtained from Equation (2) using explicit Euler Methods, where p (n) is an approximation to p (n+1) .The boundary conditions for v (n+1) are the same as those for the final velocity v (n+1) .Details of boundary conditions for full cells (F), outflow cells (O) and free surface cells (S) are provided in detail in Tomé et al. [4,30] and will not be presented here for the sake of conciseness.It can be shown that v (n+1) possesses the correct vorticity at time t n+1 , but it does not conserve mass in general.Therefore, there is a potential function ψ (n+1) so that, 3. Solve the Poisson equation for the potential function ψ for every F-cell in the domain, The boundary conditions required for solving this Poisson equation are the homogeneous Neumann conditions for rigid walls and inflows, while homogeneous Dirichlet conditions are used at outflows.4. Compute the final velocity field from Equation (13); 5. Compute the final pressure field (see [4]) by Details of the discretization of the equations (temporal and spatial) considering all types of cells (see Figure 1b) are given in [4,27,30].
Step 2 -Calculation of the extra stress tensor τ(x, t n+1 ) and free surface update To calculate the extra stress tensor τ(x, t n+1 ), initially, the Finger tensor is updated at t n+1 for every full cell (F) and surface cell (S) for every intermediate time t using Equation (7).Details of the calculation of the Finger B tensor (in two dimensions) can be found in [32] and will not be presented here because the extension to three dimensions is straightforward.Note that, for each computational cell (F and S), Equation ( 7) is solved N times (for each t ), considering each of the components of the Finger tensor.Thus, considering three-dimensional flows, the computational cost to obtain the deformation history in each cell is high, demanding a great deal of memory and simulation time (since, for each cell, it is necessary to calculate the Finger tensor N times for each component of the deformation matrix).For inflow cells (I), boundary cells (B) and empty cells (E), the Finger tensor is the identity tensor.In the outflow, the Neumann condition is assumed.
The definition of the points t for the calculation of the components of the Finger tensor and the tensor τ are given as follows.Let t j , j = 0, 1, • • • , N, be (N + 1)-points in the interval [0, t n+1 ].Then, the constitutive equation can be written in the form where an even N is assumed.For t < 0, B t (t n+1 ) = B 0 (t n+1 ), and, therefore, the first integral becomes 0 and can be solved exactly.
Each integral under the summation operator is approximated by the 3-points quadrature formula The coefficients A 0 , A 1 , A 2 are obtained by solving the (3 × 3) linear system and are found to be One of the key issues of the deformation fields method is how the integration nodes 0 = t 0 < t 1 < • • • < t N = t n+1 are distributed over the interval [0, t n+1 ] because such distribution can affect the accuracy of the results when solving complex flows.In this work, we used an ad hoc methodology for the discretization of the interval [0, t n+1 ], where a geometric progression is employed to calculate the integration nodes.Note that we consider time-dependent flows, and, therefore, the integration nodes are calculated at time t n+1 as follows: (a) Set t 0 = 0 and t The last step in the calculation is to update the position of the moving free surface (the S-cell in Figure 1b).The fluid surface is represented by a piecewise linear surface composed of triangles and quadrilaterals having marker particles on their vertices (see [30]).The particle coordinates, stored at each time step, are updated, solving the equation by Euler's method.With the new coordinates of each marker particle, a reclassification of the free surface cells is performed.A free surface cell can become an empty cell (E-cell in Figure 1b) or a full cell (F-cell in Figure 1b) or remain an S-type cell.Details on the marker particles that define the free surface and the steps for inserting and removing particles will not be shown here, but the reader can consult Tomé et al. [30] or Castelo et al. [31].

Semi-Analytical Solution
We will now derive a semi-analytical solution for a fully developed three-dimensional tube flow of a K-BKZ-PSM fluid to validate the numerical implementation.Due to the complexity of the integral model, some simplifications need to be assumed to develop the analytical solution, which is only possible for some types of domains.We consider cylindrical coordinates (see Figure 2) for simplicity and assume pure shear flow.After finding the semi-analytical solution, we present the change in variables to obtain the solution in Cartesian coordinates.We assume that r ∈ [0, 1], u = 0, w in = w(r), γ = ∂w ∂r , and The invariants I 1 and I 2 that are required in the Papanastasiou function H(I 1 , I 2 ) take the form and the tensor components are given by Taking the change in variables s = t − t , these equations are rewritten as Thus, the equations of continuity and balance of momentum become − ∂p ∂z Integrating Equation ( 26), we obtain Thus, and Equation ( 27) is rewritten as The left hand side of Equation ( 30) is just a function of r, so F must be constant, let us say C, and therefore C = dp/dz.In this way, it follows that where τ rz (r = 0) should be finite and therefore h(z) = 0, leading to The second equation in Equation ( 25) can then be rewritten as The inlet boundary condition allows one to determine the constant C. The inflow velocity w in (r) is given by w in (r) = 1 − r 2 , where w in (0) = 1 and w in (1) = 0, (34) thus, By mass conservation, and integrating by parts leads to Thus, to determine C, we must obtain γ(r) from Equation (33) and verify that The steps to calculate semi-analytical solutions are as follows: Step 1: Set an interval C 0 , C 1 such that F(C 0 ) × F(C 1 ) < 0.
Step 2: Determine the zero for |F(C)| taking |F(C)| < , where is a small value ( is the tolerance for the error).We carefully selected the value of to ensure the attainment of a semi-analytical solution accurate to six significant digits.Using Gauss-Laguerre quadrature in Equation ( 33), obtain γ(r).Using Equation (39), obtain the value of F(C) using Simpson 1/3 quadrature.
The solution in three dimensions is obtained by making the change in coordinates as follows: Thus, in three-dimensional Cartesian coordinates, we have that (41)

Results
The numerical code will now be used to solve confined (see Section 5.1) and free surface flows (see Section 5.2).

Confined Pipe Flows
In this confined pipe flow, the fluid is assumed to have only one relaxation mode.Therefore, it is possible to compare the simulation results with the semi-analytical solution presented before.The parameters used in this simulation are (see Tomé et al. [32] and Quinzani et al. [35]): (δx = δy = 0.01 20 ), M4 = 24 × 24 × 120 (δx = δy = 0.01 24 ) and M5 = 28 × 28 × 140 (δx = δy = 0.01 28 ). Figure 3a,b show the velocity profiles w and u, respectively, with corresponding cross-section velocity distributions in the plane xz (obtained mesh M3).The velocity profiles are fully developed at t = 20 s, suggesting that they have reached a steady-state condition.The velocity profile w shows a parabolic shape in the cross-section represented in Figure 3a, while the influence of the inflow can be observed in the cross-section depicted in Figure 3b, but the solution for u exhibits the expected physical behavior outside this region.Notice that, in all full cells, the initial velocity vector is defined as v = (u, v, w) = (0, 0, 0).Similar to the velocity vector, initial conditions for pressure and tensors are set to zero in full cells.As shown in Figure 4a,b, the pressure and τ xz tensor profiles are in agreement with the physically expected profiles, i.e., linear profiles across the longitudinal direction (flow direction) and transverse direction (perpendicular to the flow), respectively.In addition, the values obtained for τ xz are comparable with the behavior of the analytical solution (see Figure 5b).In the cross-section represented in Figure 4b, there is an influence of the inflow (tensors are defined as zero in the inflow, and the Finger tensor is defined as the identity matrix), but, outside this region, the expected linear profile is obtained as previously stated (for further details, refer to Figure 5b).Figure 5a-c show, respectively, the solution for the velocity w and stress components τ xz and τ zz using meshes M1-M5.The profiles were obtained at y = 0 and were taken in the center of the pipe, z = z max /2 for t = 10 (or t = 25 s).The results were compared with the semi-analytical solution, with good agreement between the numerical results (M1-M5) and the semi-analytical results.The instant t = 10 (t = 25 s) was chosen because the velocity and stresses already show a steady state behavior (the velocity residual is small (see Equation ( 42)).The residual (for the velocities) is calculated as where t is the simulation time, ∆t is the time step and N c is the number of computational cells.The slight variances between the analytical and numerical solutions can be attributed to approximations made in the numerical simulations, which differ from the precise analytical solution.Although the tube's length appears adequate for the complete development of velocity and shear stress profiles, this completeness is not reflected in the τ zz tensor.Consequently, these disparities remain minor.It is worth noting that, across all meshes, the average relative error remains below 5%. Figure 6 shows the calculation of the residuals R es in meshes M1-M5 up to time t = 10 (or, equivalently, t = 25 s).It can be observed that the residuals in the five meshes M1-M5 decrease and show convergence towards a steady-state solution, thus proving the robustness of the numerical method.As expected, we also observe a smaller residual for the most refined meshes.

Free Surface Flows
In this subsection, we test the numerical method's robustness by simulating the extrudate swell phenomenon of Boger fluids.Please note that our aim is not to conduct an in-depth study of this type of flow in Boger fluids; instead, we focus on assessing the reliability of the numerical approach.
The phenomenon known as extrudate swell is very present in various industrial processes.In this phenomenon, the fluid flows over a pipe/die and swells outside the free surface region (the cross-sectional area of the extrudate-the material being extruded-is larger after exiting the die compared to the die orifice).This behavior is mainly due to the elastic recovery of the polymer chains after being subjected to high pressures and shear forces during the extrusion process.Figure 7a shows the domain used in the simulation, and Figure 7b illustrates the phenomenon of extrudate swell (with the contour lines representing the trajectory of the fluid's free surface).It should be remarked that simulating extrudate swell can be challenging due to several numerical difficulties, which are now outlined: the extrusion process involves highly non-uniform and complex flow patterns, especially near the die exit.These flows experience rapid changes in pressure and velocity, making it difficult to model accurately; simulating extrusion processes requires discretizing the computational domain into smaller elements or cells, and the geometry of the die can be quite intricate.Moreover, the simulation must maintain numerical stability, which can be problematic in high-pressure and high-shear regions; simulating extrudate swell is computationally intensive, especially for large-scale industrial extrusion processes and using integral models; extrudate swell is a time-dependent phenomenon as the material continuously deforms and recovers during extrusion.Capturing this transient behavior accurately in numerical simulations requires precise time-stepping algorithms and may increase computational complexity.To address these challenges, researchers often resort to simplifications and assumptions to reduce computational complexity.However, in the context of this work, we take a different approach by considering the complete system of equations and accounting for the full 3D geometry.
To test the robustness of the new numerical procedure, two different Weissenberg numbers were considered (case C1 − Wi = 0.43 and case C2 − Wi = 0.64 ), both using nonshear-thinning highly elastic polymer solutions (Boger fluids − see Table 1).Boger fluids are a type of dilute polymer solution known for their remarkable elasticity, particularly at low apparent shear rates, and this unique characteristic gives rise to a significant extrudate swell during the extrusion process [14,[36][37][38].This makes Boger fluids ideal to test the numerical implementation.Numerical simulation of extruded swelling in two dimensions using the data used here (see Table 1) is presented in Mitsoulis [14].
Table 1.Parameters of the fluid used in the extrudate swell problem (see Mitsoulis [14]).The following parameters were used in the simulations (see Mitsoulis [14] and Tomé et al. [32]): Figure 8 shows the flow development of a Boger fluid for two different Weissenberg number values (cases C1 and C2).We conducted flow measurements at six different time points to analyze the behavior of the fluid in the system.The first four time instants were identical for both C1 and C2 cases, while the last two time points differed.Specifically, for the lower inlet velocity case, we considered time points t = 5 and 6 s, and, for the other case, the time points were t = 4 and 4.4 s.During the initial stages of both cases, the fluid exhibited a smooth flow with a parabolic profile as it exited the tube.However, as the process continued, swell occurred, causing the cross-sectional area of the extrudate to increase after leaving the die.This swelling behavior significantly affects the flow dynamics and needs to be carefully considered in the analysis.The simulation results show the development of the fluid front as it reaches the wall.Notably, there are distinct differences in swelling between two specific time instances: t = 6 s (C1) and t = 4.4 s (C2).
To characterize the extrudate swell phenomenon, an important parameter is the dimensionless swelling rate χ = χ max /(2r), where χ max is the maximum swelling value and r is the pipe radius.For case C1, the maximum swelling value is found to be χ = 1.78, while, for case C2, the swelling rate increases to χ = 1.95.This was expected since the Weissenberg number represents the ratio of the characteristic time scale of the elastic forces acting on a fluid to the characteristic time scale of viscous forces, and, when a polymer melt is subjected to shear flow (for example, in an extruder), the long polymer chains experience deformation due to the flow-induced stretching and alignment.A higher Weissenberg number indicates a more elastic behavior of the polymer melt, leading to more significant elastic recovery and increased extrudate swell.We may therefore conclude that the numerical method employed in the simulations demonstrates its capability to capture the transient physics of the extrudate swell problem in detail, even for a small difference in the Weissenberg number (cases C1 and C2).It accurately predicts the swelling behavior and allows for a better understanding of the process dynamics.By accounting for the material properties and flow conditions, the simulation provides valuable insights into the extrusion process and contributes to the optimization of extrusion operations.

Conclusions
In this work, a novel numerical method was developed to address three-dimensional unsteady free surface flows incorporating integral viscoelastic constitutive equations, specifically the K-BKZ-PSM (Kaye-Bernstein, Kearsley, Zapas-Papanastasiou, Scriven, Macosko) model.To implement this new approach, we integrated it into the FREEFLOW-3D code [27], enhancing its capabilities for handling viscoelastic fluid behavior.
To validate the numerical methodology, we conducted simulations of K-BKZ-PSM fluids in a pipe.The results were compared with a newly derived semi-analytical solution, and we found that the simulations performed on five different meshes yielded excellent agreement with the analytical solution.Furthermore, we applied our methodology to tackle flows with free surfaces.One notable example was the simulation of the classic extrudate swell problem, which involved a highly elastic polymeric solution known as the Boger fluid.
The significance of this work lies in the scarcity of literature concerning the simulation of unsteady three-dimensional flows of K-BKZ-PSM fluids (and integral viscoelastic models in general) using finite differences, especially when considering problems with moving free surfaces.Therefore, we hope that our contributions will inspire and encourage other researchers to further develop and explore the numerical methods we have presented here.
In conclusion, our newly developed numerical method has proven its effectiveness in handling complex fluid flow scenarios, including free surface flows such as extrudate swell simulations.The successful validation against analytical solutions reinforces the reliability of our approach and opens up opportunities for broader applications in the field of viscoelastic fluid dynamics.

Figure 1 .
(a) Typical three-dimensional staggered cell and (b) illustration of cell type classification used.

Figure 2 .
Figure 2. Representation of the pipe and a section in the rz plane.

Figure 3 .Figure 4 .
Figure 3. Velocity components u and w of v(u, v, w) along the plane xz (y = 0) at t = 20 s.(a) Visualization of the velocity profile w.(b) Visualization of the velocity profile u.

Figure 5 .
Comparison between the analytical and numerical solutions for (a) w velocity component and tensors components (b) τ xz and (c) τ zz .

Figure 7 .
Schematic of a free surface simulation in the FREEFLOW-3D software.(a) Schematic representation of the domain; (b) illustration of the extrudate swell.The fluid exits the tube and starts to swell.

Figure 8 .
Figure 8. Flow development of a Boger fluid for two different inlet velocities (cases C1 and C2).