A Fully Implicit Log-Conformation Tensor Coupled Algorithm for the Solution of Incompressible Non-Isothermal Viscoelastic Flows

In this work, a block-coupled algorithm is presented, which can compute laminar, incompressible, non-isothermal, viscoelastic flow problems based on the log-conformation tensor approach. The inter-equation coupling of the incompressible Cauchy linear momentum and mass conservation equations is obtained in a procedure based on the Rhie–Chow interpolation. The divergence of the log-conformation tensor term in the linear momentum equations is implicitly discretized in this work. In addition, the velocity field is considered implicitly in the log-conformation tensor constitutive equations by expanding the advection, rotation and the rate of deformation terms with a Taylor series expansion truncated at the second-order error term. Finally, the advection and diffusion terms in the energy equation are also implicitly discretized. The mass, linear momentum, log-conformation tensor constitutive model and energy-discretized linear equations are joined into a block-matrix following a monolithic framework. Validation of the newly developed algorithm is performed for the non-isothermal viscoelastic matrix-based Oldroyd-B fluid flow in the axisymmetric 4:1 planar sudden contraction benchmark problem.


Introduction
The polymers' processing techniques are predominantly non-isothermal, such as injection molding [1][2][3], heat exchange problems [4,5], or in plastication, including heating and cooling sequences [6,7]. The thermal conductivity and heat transfer are usually low in this processes; however, due to the heating or cooling of the machine's operations, large temperature gradients arise in the fluid [4,8]. In addition, the viscoelastic behavior of polymers acts on the temperature field as well as on the fluid deformation [4,8]. Therefore, flow properties are strongly dependent on both rheology and temperature; and, thus, it is essential to understand and make predictions regarding non-isothermal viscoelastic fluid flows.
The temperature dependence of linear viscoelastic properties (such as the relaxation time λ) can be included in constitutive equations using the time-temperature superposition principle [9]. In this way, the material properties can be defined through a function of the temperature, the so-called shift factor [10]. Two empirical correlations of the shift factor are widely employed: the William-Landel-Ferry (WLF) [11] and Arrhenius [12] models. Thus, the temperature is considered an independent variable in the constitutive equations employed to compute the components of the polymeric stress tensor (see the work of Peters and Baaijens [13] for a detailed discussion on this topic). In addition, when solving non-isothermal viscoelastic flows, the internal energy of the fluid is not only a function of the temperature [13]. The conversion mechanisms of internal energy need to be taken into account for non-isothermal viscoelastic flows; specifically, the thermal energy is partly can be found in Afonso et al. [37]. In this work, we will employ the log-conformation tensor approach to handle the HWNP.
In this manuscript, a new numerical code is developed in the context of the Finite Volume Method (FVM), following a monolithic framework to compute the non-isothermal flow of viscoelastic fluids. To the author's knowledge, the other CFD codes, which provide a fully implicit block-coupled solution for a discretized, log-conformation, viscoelastic, fluidflow system, were developed by Knechtges [38] using the Finite Element Method (FEM) and Spahn [39] using FVM; however, these studies did not consider the non-isothermal effects. In this work, the solution to the enlarged system of equations, composed of continuity, linear momentum, log-conformation tensor constitutive equation and energy, is obtained using a Bi-Conjugate Gradient Stabilized solver. The validation of the fully implicit, blockcoupled, non-isothermal, viscoelastic, log-conformation tensor-based solver is performed for the Oldroyd-B fluid flow in the axisymmetric 4:1 planar sudden-contraction benchmark problem. For assessment purposes, the results obtained with the newly-developed code are compared with numerical results found in the scientific literature. We study flows at a high Weissenberg number and we investigate the limits of pure energy elasticity and pure entropy elasticity. Lastly, we also analyzed the effect of the jump in wall temperature near the contraction for positive and negative increments.
The remaining sections of the manuscript are organized as follows. In Section 2, the governing equations for the stress tensor and log-conformation tensor-based formulations of non-isothermal viscoelastic flows are presented. Subsequently, in Section 3, the numerical procedure of the block-coupled algorithm will be described in detail, including the finite-volume discretization process for all the implicit terms considered in the governing equations. In Section 4, the validation of the newly-developed numerical algorithm is performed, and in Section 5 the main conclusions of the work are addressed.

Governing Equations
In this section, the equations that involve non-isothermal viscoelastic fluid flow problems are presented for both stress-tensor-and log-conformation tensor-based formulations.

Stress-Tensor-Based Formulation
The governing equations for laminar, incompressible, non-isothermal viscoelastic flows are the conservation of mass and linear momentum, together with a constitutive equation modeling the polymer stress behavior and the energy equation to account for thermodynamical effects.
The conservation of mass and linear momentum equations read as follows: where Einstein's summation convention applies, u i are the velocity components along the Cartesian co-ordinates x i , ρ is the fluid density, t is the time, p is the pressure and τ ij are the components of the total extra-stress tensor (i, j = 1, 2 for 2D flows), which is split into Newtonian (solvent), (τ N ) ij , and elastic (polymeric), (τ E ) ij , contributions, such that τ ij = (τ N ) ij + (τ E ) ij . The calculation of the stress terms is completed using the following relations: 1 where η N (T) and η E (T) are the temperature-dependent solvent and polymeric viscosities, respectively, D ij is the strain rate tensor, which describes the rate of stretching and shearing, λ(T) is the temperature-dependent polymer relaxation time, I ij is the identity tensor and h((τ E ) ij ) is a tensor that can be given by different expressions, related to the constitutive equation chosen to model the viscoelastic fluid. For the Oldroyd-B model [40], where κ is a positive constant, the so-called mobility factor, which is related to the elongational behavior of the fluids. For the Phan-Thien-Tanner (PTT) model [42,43], the tensor h is of the form ii is the trace of the polymeric stress tensor and is a material parameter called the extensibility factor, related to the fluid behavior in extensional flow. In addition, the Giesekus and PTT models present one more non-linearity, which is given by the product h((τ E ) ij )(τ E ) ij . This term is responsible for the shear-thinning, the non-zero second normal stress coefficient and the stress overshoot in transient shear flows of viscoelastic fluids. In this work, we will provide a preliminary assessment of the merits of the fully implicit, block-coupled, non-isothermal, viscoelastic, log-conformation tensor-based algorithm for calculations using the Oldroyd-B fluid model, which is commonly used to validate newly-developed viscoelastic codes due to the stress singular behavior near sharp corners or at stagnation points. For these models, a characteristic (solvent) viscosity ratio can be defined by β = η N /(η N + η E ) = η N /η 0 , known as the retardation ratio, where η 0 is the total viscosity in the limit of vanishing shear rate.
Following the work of Peters and Baaijens [13] the energy balance equation for the case of viscoelastic flows is as follows: where k is the thermal conductivity of the fluid, without dependence on temperature T and polymer orientation, C p is the specific heat capacity, also without temperature and polymer orientation dependence [44], and α is the energy partitioning coefficient. When α = 1, all mechanical energy is dissipated as heat (pure entropy elasticity), and if α = 0, all mechanical energy is stored as elastic energy (pure elastic material) [13,18]. Habla et al. [18] concluded that the effect of the parameter α is negligible because, with a fully developed shear flow, only stress work occurs, and the internal energy storage is absent. In addition, for the Giesekus model. For the Oldroyd-B model calculations considered in the validation section of this work (Section 4),λ(T) = λ(T). The temperature dependencies of the relaxation time, λ(T), solvent and polymeric viscosities, η N (T) and η E (T), respectively, are given by where T 0 is a reference temperature and a T is a shift factor obeying the Williams-Landel-Ferry (WLF) relation: in which C 1 and C 2 are the WLF parameters and T 0 is the reference temperature. Frequently used sets of WLF parameters (C 1 , C 2 ) are (5, 150) for temperatures relatively far from the glass transition temperature T g , enabling the thermorheological coupling, and (15,50) for temperatures near T g [18].

Log-Conformation Tensor-Based Formulation
In this section, we write the viscoelastic stress tensor-based formulation in terms of the log-conformation tensor variable, which was proposed by Fattal and Kupferman [35] to address the HWNP. For that purpose, the polymeric stress tensor, (τ E ) ij , is related to the conformation tensor, σ ij , by the following equation Subsequently, the conformation tensor σ ij is replaced by its matrix logarithmic Ψ ij = log(σ ij ), and Equations (2), (4) and (5) are substituted by 1 λ(T) e ξ m P m , with d as the dimension of the physical space (d = 2 for the calculations performed in this work), ξ m the eigenvalues of Ψ ij and P m is the projection matrix onto the corresponding eigenspace. Therefore, ifê m is the eigenvector corresponding to ξ m , then P m =ê m ⊗ê m [38]. In addition, note that λ(T 0 ) and η P (T 0 ) are known values for the reference temperature T 0 , and, therefore, the quotient (11) and (13) is a constant, considering that λ and η E scale in the same way [45].

Numerical Method
In this section, we will describe a finite volume numerical method to set up a block-coupled solver procedure to simultaneously solve the continuity (Equation (1)), linear momentum (Equation (11)), log-conformation tensor (Equation (12)) and energy (Equation (13)) equations.
Within the framework of the block-coupled solver developed in this work, the advection, pressure gradient, diffusion and log-conformation tensor terms in the conservation of linear momentum equations are implicitly discretized (see Section 3.1). Subsequently, the velocity field term in the conservation of mass equation is also treated in an implicit manner (see Section 3.2). In addition, and as an extension to our previous work [28], where we have discretized implicitly the advection term in the stress constitutive equation, here the rotation and the rate of deformation terms are implicitly discretized (see Section 3.3). Lastly, the advection and diffusion terms in the energy equation are also implicitly discretized (see Section 3.4). The rate of change term in all the equations is implicitly discretized using the backward implicit Euler scheme.

Discretization of the Equations for Conservation of Linear Momentum
In the framework of the FVM, the discretization process starts by integrating the conservation of linear momentum equations (Equation (11)) over a general control volume (also called representative volume or cell) V P , where the subscript P refers to values of the variables at cell with centroid P, as shown in Figure 1, to yield where the additional terms involving η are related to the improved both-side diffusion technique [46], which can solve the checkerboard pattern due to numerical instabilities caused by a velocity-stress decoupling. Note that we also used the identity ∂ ∂x j ∂u j ∂x i = 0. The following step of the discretization process is to apply the Gauss divergence theorem to transform the volume integrals of the advection, pressure and diffusion terms in Equation (14) into surface integrals as follows: where S is the boundary of control volume V P and n i are the components of the outward pointing unit vector normal to S.
Subsequently, a second-order integration scheme is employed to approximate the surface integrals and the following linear momentum semi-discretized equations are obtained: where ∆t is the time-step, the superscript 0 represents the previous time step value, nb refers to values at the faces f , obtained by interpolation between P and its neighbors, and s i are the components of the area normal vector to face f . Finally, the linear momentum semi-discretized equations are transformed into an algebraic linear system of equations by expressing the variation in the dependent variable and its derivatives in terms of the control volume P and its neighbors' values at the respective centroids, such as where a u i φ P and a u i φ F are the owner and neighbor coefficients in the discretized linear momentum equation representing the velocity component u i and the variable φ interactions, respectively; b u i P is the source term, where NB(P) refers to the neighbors of the controlvolume with centroid P.
For the sake of conciseness, the contributions of the rate of change, advection, pressure gradient, implicit diffusion and explicit diffusion terms shown in Equation (16) can be found in Fernandes et al. [28]. Regarding the contribution of the log-conformation tensor term, , to the linear momentum equations, we will employ an implicit discretization by considering the following Taylor approximation for the functional (e Ψ − I) ij [28,39] where the derivative of the functional (e Ψ − I) ij in order to the log-conformation tensor variable is substituted by the expression found in Knechtges [38] for the finite element method, and p m ik and p n lj are the coefficients of the projector matrix belonging to the m-th and n-th eigenvalues (λ m , λ n ) of Ψ 0 ij . Therefore, the divergence of the functional (e Ψ − I) ij can be written as: Subsequently, applying the Gauss divergence theorem we can transform the volume integrals with derivatives into surface integralŝ and obtain the discretized form for the divergence of the functional (e Ψ − I) ij (i.e., the log-conformation tensor term) in the linear momentum equations as follows: Lastly, the contributions of the divergence of the log-conformation (DLC) tensor term for the linear momentum equations are given by where w f are the interpolation weights.
It should be noted that the total contribution for the owner and neighbor coefficients related to the linear momentum equations' interactions is given by the sum of the rate of change, advection, pressure gradient, implicit diffusion and divergence of log-conformation tensor terms. In addition, the total contribution for the explicit coefficient related to the linear momentum equations is given by the sum of the rate of change, explicit diffusion and divergence of log-conformation tensor terms.

Discretization of the Equation for Conservation of Mass
Following the same steps as presented in Section 3.1, we begin the discretization of the continuity equation, Equation (1), with the integration over the control volume V P as follows:ˆV Subsequently, by employing the divergence theorem, the volume integral is transformed into a surface integral as follows: Then, by applying a second-order integration scheme to approximate the surface integral, we can write the semi-discretized form of the continuity equation as: In a collocated framework, the velocity at the face is obtained by reconstructing a pseudomomentum equation at the face from the linear momentum equations of the straddling cells P and F, known as the Rhie-Chow interpolation [47]. For the sake of conciseness, the derivation of the discretized form for the equation of conservation of mass will not be given in detail, but it can be found in our previous work [28], which reads as Finally, we can write the mass conservation algebraic equation as: where a pφ P and a pφ F are the owner and neighbor coefficients in the discretized mass conservation equation representing the pressure field and the variable φ interactions, respectively; b p P is the source term.
The implicit pressure gradient term is discretized (see page 86 of [48]) as follows: where the coefficients of the implicit pressure gradient term for the mass conservation equation are given by The implicit coefficients for the second term in Equation (26) (corresponding to the velocity contribution) are given by Lastly, the coefficients of the explicit pressure gradient term contribution for the mass conservation equation are given by

Discretization of the Log-Conformation Tensor Constitutive Equations
The constitutive equations, Equation (4), can be written, without loss of generality, for a Giesekus fluid model by (see Theorem A.42 in [38]) where commutator term. Following Knechtges [38], f (ad(Ψ ij )) is defined by the application of the function f (x) = x/2 tanh(x/2) to the adjoint operator ad(Ψ ij ), such as where y is a symmetric matrix satisfying ad(Ψ ij )y = [Ψ ij , y], and y is continuously differentiable (C 1 ). Notice that, Equation (32) simplifies to the constitutive equation for an Oldroyd-B fluid model when κ = 0. The discretization of the log-conformation tensor constitutive equations, Equation (32), starts with the integration over the control volume V P to yield This leads to the algebraic equation with the following form: where a are the owner and neighbor coefficients in the discretized logconformation tensor constitutive equations representing the tensor component Ψ ij and the variable φ interactions, respectively, and b Ψij P is the source term. Again, for the sake of conciseness, the discretization of the rate of change terms in Equation (32) are not detailed here because similar discretizations were performed for the polymeric stress-tensor (τ E ) ij (see previous work [28]).
Regarding the commutator term [Ψ ij , Ω ij ], we can write the following expanded form: and, subsequently, we can use Taylor approximations such as and Starting with the terms in Equation (37), the negative commutator terms, the first contribution is explicit, being given by The second contribution is implicit in Ψ ij and explicit in ∂u i ∂x j , being given by (a The third contribution is implicit in ∂u i ∂x j and explicit in Ψ ij ; therefore, we need an implicit discretization of the gradient operator for velocity, which requires the integration by parts of this term, such as: Applying the Gauss divergence theorem, the discretized form of the terms on the right-hand side of Equation (41) is Using linear weighted interpolation, we can write the contributions of the third term as: Following the same reasoning as given above, the terms in Equation (38), the positive commutator terms, generate the following coefficients: Lastly, we implicitly discretize the adjoint operator f (ad(Ψ)) ijkl D kl by considering the following Taylor approximation [39] f (ad(Ψ)) ijkl D kl where is the derivative of the adjoint operator with respect to D kl , given by [38] and is the derivative of the adjoint operator with respect to Ψ kl , given by [49] ∂ f (ad(Ψ 0 )) iqrj D 0 Note that, if λ m is equal to λ n , then the denominator of Equation (49) needs to be replaced by f (λ m − λ n ) [38]. Subsequently, we integrate the adjoint term over a control volume V P as follows: and substituting D kl = 1 in Equation (50), we obtain As the fourth term on the right hand side of Equation (51) contains implicit velocity gradients, we employ integration by parts to linearize them, obtaininĝ Thus, the discretized form of Equation (52) can then be written as: The contributions of the adjoint term for the log-conformation tensor constitutive equations, using linear weighted interpolation, read as (a Again, it should be noted that the total contribution of the owner and neighbor coefficients related to the log-conformation tensor components' interactions is given by the sum of the log-conformation tensor, rate of change, advection, commutator and adjoint terms. In addition, the total contribution for the explicit coefficient related to the logconformation tensor constitutive equations is given by the sum of the rate of change, commutator and adjoint terms.

Discretization of the Equation for Conservation of Energy
The discretization starts by integrating the equation for the conservation of energy (Equation (13)) over a general control volume V P , to yield Using the Gauss divergence theorem, the volume integrals of the advection and diffusion terms in Equation (55) are transformed into surface integrals as: The semi-discretized equation for the conservation of energy is obtained by evaluating the surface integrals using a second-order integration scheme and approximating the rate of change term with a backward implicit Euler scheme, such as This leads to the algebraic equation for the energy balance with the following form: a TT P T P + a Tu P u P + a Tv where a Tφ P and a Tφ F are the owner and neighbor coefficients in the discretized conservation of energy equation, representing the temperature T and the variable φ interactions, respectively, and b T P is the source term. The rate of change (rchg) term in Equation (57) (first term), contributes to both the diagonal of the system of equations and to the explicit term as: Then, and in the framework of the FVM, the advection term in Equation (57) (second term) is linearized by computing the mass flow rate at control-volume face f (ṁ f = (s i ρu i ) f ) using the previous iteration values. Here, we used the UDS differentiating scheme to approximate the advection term. However, many high-order schemes could be used, such as the MINMOD or CUBISTA schemes. For the sake of readability, the discretization procedure will be presented for the UDS scheme, but it is important to stress that the methodology is independent of the adopted discretization scheme. In this way, the coefficients of the advection (adv) term contribution for the conservation of energy equation are given by: where the superscript Tφ means the influence of the φ variable in the T energy equation. The term max(ṁ f , 0) represents the maximum ofṁ f and 0, where the mass flux is positive if it goes from owner to neighbor cells, i.e., leaves the control-volume V P . The implicit diffusion (idi f f ) contribution, third term of Equation (57), is discretized, taking a linear profile (see page 86 of [48]) as where (d PF ) i is the vector joining the centroids P and F (see Figure 1). The coefficients of the implicit diffusion term for the conservation of energy equation are given by Finally, the coefficients of the explicit term contribution (right-hand side of Equation (57)) for the conservation of energy equation are given by

Block-Coupled Algorithm
Combining the discretized conservation of linear momentum (Equation (17)), conservation of mass (Equation (27)), log-conformation tensor (Equation (35)) and conservation of energy equations (Equation (58)), the following linear system of equations, written in matrix form, is obtained for each control volume: The linear systems (Equation (64)) obtained for each control volume of the computational domain are merged in a full system of equations, which can be written in the form AΦ = b where all variables Φ = (u i , p, Ψ ij , T) are solved simultaneously. In this procedure, all variables in the different equations are treated implicitly, which is expected to be advantageous to the stability of the overall calculation process. The fully implicit coupled algorithm can be summarized into the following steps: 1. Initialize the fluid variables with the latest known values (u n i , p n , Ψ n ij , T n ). (17), (27), (35), (58)) and solve for u i , p, Ψ ij and T (Equation (64)).

Iterate until convergence.
For the solution of the global system of discretized algebraic equations, it is fundamental that an efficient linear solver is used to obtain the best overall convergence. In this work, the iterative solver Bi-Conjugate Gradient Stabilized (BiCGStab) [50], combined with an LU preconditioner, was used to retrieve the solution of the global system of discretized algebraic equations (see detailed discussion in Pimenta and Alves [30]). The initial residual for each iteration is evaluated based on the current values of the field, before solving the block-coupled system. After each block solver linear iteration, the residual is re-evaluated (final residual). When the maximum number of linear iterations (in this work defined as 100) or the final residual falls below the solver absolute tolerance (set as 10 −6 ), the block-coupled system current iteration stops and advances in time. All computations are performed on a computer with a 2.00-GHz 64 cores AMD EPYC 7662 CPU processor and 128 GB of RAM.

Results and Discussion
The validation of the newly-developed, fully implicit, block-coupled, non-isothermal, viscoelastic, log-conformation tensor-based algorithm was performed for the laminar, incompressible, non-isothermal viscoelastic flow of an Oldroyd-B fluid in an axisymmetric 4:1 sudden contraction geometry. For assessment purposes, the results computed with the newly-developed code were compared with the available data from the scientific literature [18].

Geometry, Meshes, and Initial and Boundary Conditions
An axisymmetric 4:1 sudden contraction with ratio of the radii R 1 /R 2 = 4 was chosen as test geometry (Figure 2a), because of the availability of numerical data in the literature [18]. The upstream length was l 1 = 80R 2 and the downstream length was l 2 = 50R 2 . The downstream channel height was chosen as R 2 = 0.0020604 m. The flow had a rotational symmetry that was normal to the rz−plane and, to save computational resources and reduce the CPU times, only half of the domain was considered. The characteristics of the meshes used to discretize the geometrical domain are presented in Table 1. The meshes employed in the current work (M1 and M2) resulted from a mesh convergence analysis performed by Habla et al. [18], and had a similar refinement level to the two most refined meshes (M3 and M4) therein. The expansion or contraction geometrical factors were defined for each direction as the ratio of two consecutive cells lengths ( f z = ∆z i+1 /∆z i ,with ∆z i being the length of the cell i in the z-direction). In this way, since f z > 1 in Block V (see Table 1), in the z direction, the cells expanded from left to right. Figure 2c shows the details of the level of mesh refinement at the contraction region for M2. The higher refinement that occurs near the walls and in the contraction region allowed for the highest gradients of the computed flow variables in these locations to be captured. Nz and Nr are number of cells along z and r directions, respectively, inside each block. f z and f r are the expansion/contraction ratios inside each block. NC is the number of cells in the mesh. ∆z min and ∆r min are the minimum cell size in each direction.
The following boundary and initial conditions were used for all the runs that were performed: • For velocity, no-slip at the walls, symmetry at the centerline, parabolic velocity profile at the inlet (with average velocity U z,1 = 0.00129 m/s), and a zero-gradient condition at the outlet, i.e., assuming a fully developed flow; • For pressure, the inlet and wall boundary conditions were set as zero-gradient and the centerline as symmetry boundary condition. At the outlet Dirichlet boundary condition was used, with a fixed value p = 0. Notice that, although the zero-pressure gradient specified at the inlet did not match the fully developed Poiseuille flow with the average velocity U z,1 , this inconsistency did not affect the results, because the length of the upstream channel was sufficiently large to achieve fully developed flow conditions; • For the log-conformation tensor components, zero values were assumed at the inlet, a symmetry boundary condition was used at the centerline, a linear extrapolation of the tensor components to the boundary was used at the walls, and a zero-gradient condition was imposed at the outlet; • For the temperature, a Dirichlet condition was imposed at the inlet (T inl = 462 K), a symmetry boundary condition was used at the centerline, at the upstream wall, (z < l 1 ), T w,1 = 462 K, while, for the downstream wall, (z ≥ l 1 ) the temperature T w,2 was chosen such as to give temperature jumps of ∆T = T w,2 − T w,1 = −30 K, 0 K, 30 K. A zero-gradient condition was imposed at the outlet; • All fields were set to zero at the initial time.

Numerical Parameters
The dimensionless numbers governing the flow are the Reynolds number Re, the Weissenberg number Wi, the Peclet number Pe and the retardation ratio β, defined as: where U z,2 is the mean velocity in the axial direction in the downstream channel, and λ is the relaxation time. In this case study, Re = 3.9 × 10 −5 , corresponding to creeping flow conditions. The retardation ratio was equal to β = η S /η 0 = 19/20, thus assuming the solvent contribution to be negligibly small, which approximately recovered an Upper-Convected-Maxwell model. The Peclet number was kept constant at Pe = 345 by setting c P = 1500 J/kg K and k = 0.17 W/mK. The WLF parameters were C 1 = 4.54 and C 2 = 150.36. The split coefficient varied between pure energy elasticity and entropy elasticity, α = 0 or 1, respectively. The use of a normalized time-step ∆t/(R 2 /U z,2 ) of 10 −4 allowed for converged solutions to be obtained for all the runs performed. The maximum local Courant number corresponding to the normalized time-step 10 −4 obtained for the axisymmetric 4:1 sudden contraction is 0.07.

Effects of the Energy Partitioning Parameter α
In Figure 3, the temperature profile on the line r = 0.97R 2 (Figure 3a) and the temperature contour plots (Figure 3b) are shown as a function of the axial position z/R 1 for α = 0 and α = 1 at Wi = 5 and ∆T = 0 K. As illustrated in Figure 3a, the temperature profile computed by the newly-developed, fully implicit, block-coupled, non-isothermal, viscoelastic, log-conformation tensor-based algorithm is in agreement with the results presented by Habla et al. [18]. Due to the increase in the deformation rate near the contraction, the dissipation also increases, which leads to a temperature rise shortly before the contraction. At the contraction, the fluid moves to the wall with the imposed fixed wall temperature T w,2 = 462 K and, therefore, due to heat conduction towards the wall, a fast decrease in temperature is observed. Note that this decrease is remarkably higher for entropy elasticity (α = 1) due to its higher temperature, which promotes a larger heat conduction rate. Subsequently, just after the re-entrant corner, due to the larger normal stresses developed here (see Figure 4), which lead to an increase of dissipation, the temperature rises again at the steepest rate. Further along the downstream channel, the temperatures also increase, but now at a smaller rate and, ultimately, the difference in the temperatures between energy elasticity and entropy elasticity cases is small, because the energy is now more released as pure energy elasticity (α = 0) [18]. The temperature contour plots shown in Figure 3b are similar for both the energy and entropy elasticities, as expected from the marginal differences shown in Figure 3a for the temperature profile at r = 0.97R 2 . In both cases, we see the formation of a larger temperature-rising region for z/R 1 > 1, which is also extended through the downstream channel radial direction. In Figure 4, the axial normal stress profile on the line r = 0.97R 2 ( Figure 4a) and the axial normal stress contour plots ( Figure 4b) as a function of the axial position z/R 1 for α = 0 and α = 1 at Wi = 5 and ∆T = 0 K, are shown. As illustrated in Figure 4a, the axial normal stress profile computed by the newly-developed, fully implicit, block-coupled, non-isothermal, viscoelastic, log-conformation tensor-based algorithm and the isothermal version is in agreement with the results presented by Habla et al. [18]. In the non-isothermal cases, the axial normal stress is smaller than the one obtained for the isothermal calculation. This behaviour can be attributed to the fact that increasing the temperature leads to a reduction in the viscosity value (see Equation (9)), which decreases the stress values, and also leads to a reduction in the relaxation time, resulting in a smaller local Weissenberg number and, therefore, fewer elastic effects (i.e., decrease in stresses). In addition, just after the re-entrant corner, we see an abrupt increase of the normal stresses due to the increase in the fluid deformation in this region, followed by a fast relaxation, before it starts to increase further down the downstream channel, but now at a smaller rate. The axial normal stress contour plots shown in Figure 4b are again similar for both the energy and entropy elasticities. In both cases, we see the formation of a thinner layer of normal stresses rising region at 0 < z/R 1 < 0.2, which then increases in width, but with smaller normal stress values, at 0.2 < z/R 1 < 1. Additionally, Figure 5 shows the contours of the first normal stress difference and shear stress predicted on M2 at α = 0, Wi = 5 and ∆T = 0 K, for a zoomed region around the re-entrant corner. The maximum first normal stress difference is located on the downstream channel wall near the re-entrant corner (see Figure 5a). Moreover, at the contraction vertical wall, the first normal stress difference is negative, which is responsible for fluid re-circulation and the formation of the corner vortex. The maximum magnitude of the shear stress component is also located on the downstream channel wall near the re-entrant corner (see Figure 5b); however, in that case, the shear stress value is negative, pushing the fluid towards the symmetry line at r/R 1 = 0. Finally, at the contraction vertical wall, the shear stress is positive, allowing for the extension of the corner vortex size. Contours of (a) first normal stress difference (τ P,zz − τ P,rr )/(η 0 U z,2 /R 2 ) and (b) shear stress τ P,rz /(η 0 U z,2 /R 2 ) for α = 0, Wi = 5 and ∆T = 0 K using M2.
In Figure 6, the axial velocity profile on the line r = 0.97R 2 (Figure 6a) and the axial velocity contour plots (Figure 6b) as a function of the axial position z/R 1 for α = 0 and α = 1 at Wi = 5 and ∆T = 0 K are shown. As illustrated in Figure 6a, the axial velocity profiles computed by the newly-developed, fully implicit, block-coupled, non-isothermal, viscoelastic, log-conformation tensor-based algorithm and with the isothermal version are in agreement with the results presented by Habla et al. [18]. In addition, we can see that the influence of temperature on the local velocity profile for both cases of entropy elasticity and energy elasticity is negligible, being similar to the axial velocity of the isothermal case. The axial velocity contour plots shown in Figure 6b are, again, similar for both the energy and entropy elasticities, showing the accommodation of the fluid near the re-entrant corner, i.e., the fluid is accelerated in the center while decelerating in the wall-near regions, and a fully developed velocity profile at the downstream channel. In Table 2, we provide a summary of the number of iterations and execution time of the segregated/iterative and coupled/monolithic approaches for the calculation of the non-isothermal, viscoelastic matrix-based Oldroyd-B fluid flow in the two-dimensional, axisymmetric 4:1 planar sudden-contraction geometry, using the two different meshes, M1 and M2, with α = 0 at Wi = 5 and ∆T = 0 K. The ratio of the number of iterations required by the segregated algorithm to that required by the coupled one increases from 432 to 524 for M1 and M2, respectively, with accompanying computational cost ratios of 17 and 19, respectively, which clearly shows the benefits obtained by using the fully implicit coupled approach. Table 2. Comparison of the number of iterations and CPU time required by the segregated (S) and coupled (C) solvers for the calculation of the non-isothermal viscoelastic matrix-based Oldroyd-B fluid flow in the two-dimensional axisymmetric 4:1 planar sudden-contraction geometry, using the two different meshes M1 and M2, with α = 0 at Wi = 5 and ∆T = 0 K.

Mesh
Number

Effects of Wall Temperature Jumps
The temperature and axial velocity profiles at the outlet of the downstream section are shown in Figure 7 for Wi = 5, α = 0 and temperature jumps ∆T = −30 K, ∆T = 0 K and ∆T = +30 K. The results for the case of energy elasticity at the outlet of the downstream section, for all temperature jumps, did not present any differences [18]. As shown in Figure 7a, the wall (r/R 1 = 0.25) temperature changes by more than 60 K and the centerline (r/R 1 = 0) temperature varied less than 5 K from ∆T = −30 K to ∆T = +30 K, in agreement with the results obtained by Habla et al. [18]. These temperature changes are responsible for the limited effect of external heating or cooling in the thermal control of the flow at the bulk region. In Figure 7b, the axial velocity profile at the centerline increases with the decrease in temperature jump, due to the smaller viscosity in the near-wall region. The temperature and axial velocity contours are shown in Figure 8 for Wi = 5, α = 0 and temperature jumps ∆T = −30 K and ∆T = +30 K. For the cooling case, i.e., ∆T = −30 K (top figures in Figure 8), the reduction in temperature on the downstream channel walls promotes an increase in fluid centerline velocity of approximately 3.1 times more than the one obtained for the case ∆T = 0 K (see Figure 6). For the heating case, i.e., ∆T = +30 K (bottom figures in Figure 8), the increase in temperature on the downstream channel walls promotes an increase in fluid centerline velocity of approximately 2.3 times more than the one obtained for the case ∆T = 0 K (see Figure 6). Notice that the increase in centerline velocity for the heating case is approximately 75% smaller than the one obtained for the cooling case. Additionally, Figure 9 shows the contours of the first normal stress difference and shear stress predicted on M2 at α = 0, Wi = 5, ∆T = −30 K (top) and ∆T = +30 K (bottom), for a zoomed region around the re-entrant corner. The maximum first normal stress difference was found to decrease by 35% and increase by 30% for the cases ∆T = −30 K and ∆T = +30 K, respectively, when compared to the case ∆T = 0 K. The minimum value of the shear stress is found to both increase and decrease by 50% for the cases ∆T = −30 K and ∆T = +30 K, respectively, when compared to the case ∆T = 0 K. . Contours of (a) first normal stress difference (τ P,zz − τ P,rr )/(η 0 U z,2 /R 2 ) and (b) shear stress τ P,rz /(η 0 U z,2 /R 2 ) for α = 0, Wi = 5, ∆T = −30 K (top) and ∆T = +30 K (bottom) using M2.

Conclusions
This paper presented a fully implicit coupled method (also known as a monolithic approach) for the solution of laminar, incompressible, non-isothermal, viscoelastic flows based on the log-conformation tensor framework for high Weissenberg number problems. The fully implicit coupled solver is a pressure-based method in which the pressure equation is derived through the Rhie-Chow interpolation, allowing for coupling between the pressure and velocity fields. In addition, the explicit diffusion term added by the improved both-sides-diffusion (iBSD) technique to the linear momentum equations is discretized with a special second-order derivative of the velocity field, allowing for coupling between the velocity and log-conformation tensor fields. Subsequently, the divergence of the logconformation tensor term in the linear momentum equations is implicitly discretized, and the velocity field is considered implicitly in the log-conformation tensor constitutive equations by expanding the advection, rotation and the rate of deformation terms, all by considering a Taylor series expansion truncated at the second-order error term. Finally, the advection and diffusion terms in the energy equation are also implicitly discretized.
The validation of the newly-developed algorithm was performed for the non-isothermal viscoelastic matrix-based Oldroyd-B fluid flow in the benchmark problem of a two-dimensional axisymmetric 4:1 planar sudden-contraction geometry, and the results obtained by the fully implicit coupled algorithm were shown to be in remarkable agreement with other results found in the scientific literature (less than 5% error differences), which validated the present implementation. The results were obtained at a high Weissenberg number, and allowed to study the influence of the energy splitting factor at the limit of pure energy elasticity and pure entropy elasticity, and the effect of the wall temperature jump near the contraction for positive and negative temperature increments.
In future works, the algorithm will be further assessed by looking at 3D case studies and employing non-linear viscoelastic fluid behaviors, such as the shear-thinning phenomenon predicted by the Giesekus fluid model.