**1. Introduction**

The polymers' processing techniques are predominantly non-isothermal, such as injection molding [1–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

**Citation:** Fernandes, C. A Fully Implicit Log-Conformation Tensor Coupled Algorithm for the Solution of Incompressible Non-Isothermal Viscoelastic Flows. *Polymers* **2022**, *14*, 4099. https://doi.org/10.3390/ polym14194099

Academic Editors: Patrick Ilg and Maria Bercea

Received: 5 September 2022 Accepted: 26 September 2022 Published: 30 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

dissipated and partly stored in the fluid. Therefore, the energy equation should predict which part of the mechanical power is dissipated and which part is accumulated as elastic energy [4,8]. For that purpose, an additional term is needed in the energy equation [14]. Peters and Baaijens [13] developed an internal energy equation for multiple rate-type fluids based on a constant weighting factor that characterizes the ratio of entropy to energy elasticity [15]. Several numerical studies have also used this concept [16–18] and we will employ this in the current work.

Numerical simulations can describe these complex flow mechanisms and help to gain a better understanding of, and improvements in, the processes where they occur. For that purpose, Computational Fluid Dynamics (CFD) are used to guide the theoretical researchers and practitioner engineers, through the use of both open-source [3–5,7,19] and proprietary software [20]. In the last decade, a significant effort has been made in research on the nonisothermal flows of viscoelastic fluids. A survey of the scientific literature finds different works that describe the non-isothermal viscoelastic fluid flows based on iterative numerical algorithms. For example, Shahbani-Zahiri et al. [21] studied the recirculation and thermal regions of viscoelastic flow in the symmetric planar problem for different expansion angles. Kunisch and Marduel [22] employed the finite element (FE) method to study the optimal control of non-isothermal viscoelastic fluids to minimize vortices and control the heat flux. Spanjaards et al. [23] performed a 3D transient non-isothermal simulation to predict the extrudate shape of viscoelastic fluids emerging from an asymmetric keyhole-shaped die. However, the current state-of-the-art codes depend on iterative algorithms, such as the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) procedure [24], which are known to delay the convergence of the problem of interest when compared to monolithic or coupled algorithms [25–30]. The iterative algorithms, also known as segregated algorithms, are characterized to provide a separate solution of the linear momentum, mass, viscoelastic polymer stress tensor and energy conservation equations, which are then iterated until convergence. Recently, with the increase in computational resources and due to scalability problems in the segregated algorithms, the monolithic approach has been used, with the advantage of decreasing the computational wall time of the simulation, particularly for finer meshes, as shown by Fernandes et al. [28]. Thus, a methodology based on the monolithic approach for the simulation of non-isothermal viscoelastic flows would be of major importance.

In addition, the benchmark problems of both planar and axisymmetric contraction flows are also extensively studied to evaluate the stability of newly developed numerical algorithms [31–33]. These benchmark problems are especially important because, near the contraction, complex flow profiles are generated, and thus large stress gradients are developed, which can cause numerical difficulties, leading to the overall failure of the algorithms. Bearing this in mind, we will revisit the axisymmetric sudden contraction benchmark flow to validate the newly developed, fully implicit, coupled solver for nonisothermal viscoelastic fluid flows.

The rheological data available in the literature for the validation of non-isothermal viscoelastic fluids are scarse due to the complex fluid behaviour, which generally requires several modes to capture the full range of operating conditions. In this work, we will employ the highly elastic, polyisobutylene-based polymer solution (PIB-Boger fluid), which is typically described by the quasi-linear Oldroyd-B viscoelastic fluid model. Another important issue to consider when solving viscoelastic fluid flows is the elastic effects, specifically, the flow at high Weissenberg numbers [34]. It is well known that numerical simulations tend to become unstable at increased Weissenberg numbers, the so-called High Weissenberg Number Problem (HWNP). The seminal work of Fattal and Kupferman [35] proposed a reformulation of the viscoelastic stress-tensor-based formulation to solve the HWNP, where the logarithm of the conformation tensor is used as the main variable in the constitutive transport equation. Different methods have been also used to solve the HWNP, such as the square root of the conformation tensor [36]. A detailed discussion of this topic

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.

#### **2. 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.

#### *2.1. 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:

$$\frac{\partial \mu\_i}{\partial \mathbf{x}\_i} = \mathbf{0},\tag{1}$$

$$
\rho \left( \frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial x\_j} \right) + \frac{\partial p}{\partial x\_i} - \frac{\partial \tau\_{ij}}{\partial x\_j} = 0,\tag{2}
$$

where Einstein's summation convention applies, *u<sup>i</sup>* are the velocity components along the Cartesian co-ordinates *x<sup>i</sup>* , *ρ* 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:

$$(\tau\_N)\_{ij} = \mathfrak{D}\eta\_N(T)\mathcal{D}\_{ij} = \eta\_N(T)\left(\frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i}\right),\tag{3}$$

$$\begin{split} & \frac{1}{\eta\_{\rm E}(T)} \Big( \mathcal{I}\_{\rm ij} + \mathcal{h}((\tau\_{\rm E})\_{\rm ij}) \big) (\tau\_{\rm E})\_{\rm ij} + \frac{\lambda(T)}{\eta\_{\rm E}(T)} \left( \frac{\partial(\tau\_{\rm E})\_{\rm ij}}{\partial t} + u\_{k} \frac{\partial(\tau\_{\rm E})\_{\rm ij}}{\partial \mathbf{x}\_{k}} \right) - \left( \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial u\_{j}}{\partial \mathbf{x}\_{i}} \right) - \\ & \frac{\lambda(T)}{\eta\_{\rm E}(T)} \Big( (\tau\_{\rm E})\_{\rm ik} \frac{\partial u\_{j}}{\partial \mathbf{x}\_{k}} + (\tau\_{\rm E})\_{\rm jk} \frac{\partial u\_{i}}{\partial \mathbf{x}\_{k}} \Big) = 0, \end{split} \tag{4}$$

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], h((*τE*)*ij*) = 0. For the Giesekus model [41] h((*τE*)*ij*) = *κλ ηE* (*τE*)*ij*, 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 h((*τE*)*ij*) = *eλ ηE* tr((*τE*)*ij*)I*ij*, where tr((*τE*)*ij*) = (*τE*)*ii* is the trace of the polymeric stress tensor and *e* 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*/(*η<sup>N</sup>* + *ηE*) = *ηN*/*η*0, known as the retardation ratio, where *η*<sup>0</sup> 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:

$$
\rho \mathbb{C}\_p \left( \frac{\partial T}{\partial t} + u\_l \frac{\partial T}{\partial \mathbf{x}\_l} \right) - k \frac{\partial^2 T}{\partial \mathbf{x}\_l^2} = (\tau\_N)\_{\bar{\imath}\bar{\jmath}} \mathbf{D}\_{\bar{\jmath}\bar{\imath}} + a(\tau\_E)\_{\bar{\imath}\bar{\jmath}} \mathbf{D}\_{\bar{\jmath}\bar{\imath}} + (1 - a) \frac{(\tau\_E)\_{\bar{\imath}i}}{2\bar{\lambda}(T)},\tag{5}
$$

where *k* is the thermal conductivity of the fluid, without dependence on temperature *T* and polymer orientation, *C<sup>p</sup>* 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, *λ*¯(*T*) = *λ*(*T*) 1 + *λ*(*T*) *e* (*τE*)*ii ηE*(*T*) −<sup>1</sup> for the PTT model and *λ*¯(*T*) = *λ*(*T*) 1 + *λ*(*T*) *κ* (*τE*)*ij ηE*(*T*) −<sup>1</sup> 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

$$
\lambda(T) = a\_T \lambda(T\_0),
\tag{6}
$$

$$
\eta\_N(T) = a\_T \eta\_N(T\_0) \,. \tag{7}
$$

$$
\eta\_E(T) = a\_T \eta\_E(T\_0). \tag{8}
$$

where *T*<sup>0</sup> is a reference temperature and *a<sup>T</sup>* is a shift factor obeying the Williams–Landel– Ferry (WLF) relation:

$$a\_T = \exp\left(\frac{-\mathbb{C}\_1(T - T\_0)}{\mathbb{C}\_2 + T - T\_0}\right),\tag{9}$$

in which *C*<sup>1</sup> and *C*<sup>2</sup> are the WLF parameters and *T*<sup>0</sup> 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 *Tg*, enabling the thermorheological coupling, and (15, 50) for temperatures near *T<sup>g</sup>* [18].

#### *2.2. 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

$$(\tau\_E)\_{ij} = \frac{\eta\_E(T)}{\lambda(T)} (\sigma\_{ij} - \mathbf{I}\_{ij}). \tag{10}$$

Subsequently, the conformation tensor *σij* is replaced by its matrix logarithmic Ψ*ij* = log(*σij*), and Equations (2), (4) and (5) are substituted by

$$\rho \left( \frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right) + \frac{\partial p}{\partial \mathbf{x}\_i} - 2\eta\_N(T) \frac{\partial \mathcal{D}\_{ij}}{\partial \mathbf{x}\_j} - \frac{\eta\_E(T)}{\lambda(T)} \frac{\partial \left( \mathbf{e}^{\mathbf{v}\_{ij}} - \mathbf{I}\_{ij} \right)}{\partial \mathbf{x}\_j} = \mathbf{0},\tag{11}$$

$$\frac{1}{\lambda(T)}\left[\left(\mathbf{I}\_{lj} + \mathbf{h}(\mathbf{e}^{\mathbf{V}\_{ij}})\right)(\mathbf{e}^{\mathbf{V}\_{ij}} - \mathbf{I}\_{lj})\right] + \frac{\partial \mathbf{e}^{\mathbf{V}\_{ij}}}{\partial t} + \mu\_k \frac{\partial \mathbf{e}^{\mathbf{V}\_{ij}}}{\partial \mathbf{x}\_k} - \mathbf{e}^{\mathbf{V}\_{ik}} \frac{\partial u\_j}{\partial \mathbf{x}\_k} - \mathbf{e}^{\mathbf{V}\_{jk}} \frac{\partial u\_l}{\partial \mathbf{x}\_k} = \mathbf{0},\tag{12}$$

$$\rho \mathbf{C}\_p \left( \frac{\partial T}{\partial t} + u\_i \frac{\partial T}{\partial \mathbf{x}\_i} \right) - k \frac{\partial^2 T}{\partial \mathbf{x}\_i^2} = (\tau\_\mathcal{N})\_{\bar{\mathcal{U}}} \mathbf{D}\_{\bar{\mathcal{U}}} + \frac{\eta\_E(T)}{\lambda(T)} \left[ a (\mathbf{e}^{\mathbf{Y}\_{\bar{\mathcal{U}}}} - \mathbf{I}\_{\bar{\mathcal{U}}}) \mathbf{D}\_{\bar{\mathcal{U}}} + (1 - a) \frac{(\mathbf{e}^{\mathbf{Y}} - \mathbf{I})\_{\bar{\mathcal{U}}}}{2\bar{\lambda}(T)} \right], \tag{13}$$

with *λ*¯(*T*) = *λ*(*T*) 1 + *e* (e <sup>Ψ</sup> <sup>−</sup> <sup>I</sup>)*ii*−<sup>1</sup> for the PTT model, *λ*¯(*T*) = *λ*(*T*) 1 + *κ* (e <sup>Ψ</sup> <sup>−</sup> <sup>I</sup>)*ij*−<sup>1</sup> for the Giesekus model and *λ*¯(*T*) = *λ*(*T*) for the Oldroyd-B model. Here, e <sup>Ψ</sup>*ij* is the matrix *d*

exponential function, defined as e <sup>Ψ</sup>*ij* = ∑ *m*=1 e *<sup>ξ</sup><sup>m</sup>* **P***m*, with *d* as the dimension of the physical

space (*d* = 2 for the calculations performed in this work), *ξ<sup>m</sup>* the eigenvalues of Ψ*ij* and **P***<sup>m</sup>* is the projection matrix onto the corresponding eigenspace. Therefore, if **e**ˆ*<sup>m</sup>* is the eigenvector corresponding to *ξm*, then **P***<sup>m</sup>* = **e**ˆ*<sup>m</sup>* ⊗ **e**ˆ*<sup>m</sup>* [38]. In addition, note that *λ*(*T*0) and *ηP*(*T*0) are known values for the reference temperature *T*0, and, therefore, the quotient *ηE*(*T*) *λ*(*T*) in Equations (11) and (13) is a constant, considering that *λ* and *η<sup>E</sup>* scale in the same way [45].

#### **3. 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.

#### *3.1. 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*) *VP*, where the subscript *P* refers to values of the variables at cell with centroid *P*, as shown in Figure 1, to yield

$$\begin{split} &\rho \left( \int\_{V\_P} \frac{\partial u\_i}{\partial t} \, \mathrm{d}V\_P + \int\_{V\_P} u\_j \frac{\partial u\_i}{\partial \mathbf{x}\_j} \, \mathrm{d}V\_P \right) + \int\_{V\_P} \frac{\partial p}{\partial \mathbf{x}\_i} \, \mathrm{d}V\_P - \int\_{V\_P} (\eta\_N(T) + \eta^\star) \frac{\partial^2 u\_i}{\partial \mathbf{x}\_j^2} \, \mathrm{d}V\_P - \\ & \quad - \int\_{V\_P} \frac{\eta\_E(T)}{\lambda(T)} \frac{\partial (\mathbf{e}^\mathbf{v} - \mathbf{I})\_{ij}}{\partial \mathbf{x}\_j} \, \mathrm{d}V\_P = - \int\_{V\_P} \frac{\partial}{\partial \mathbf{x}\_j} \left( \eta^\star \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right) \, \mathrm{d}V\_{P\prime} \end{split} \tag{14}$$

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 *<sup>∂</sup> ∂x<sup>j</sup> ∂u<sup>j</sup> ∂x<sup>i</sup>* = 0.

**Figure 1.** Schematic representation of the control volume *V<sup>P</sup>* with centroid *P* (owner), with distance vector to the origin **r***P*, and neighboring control volumes with centroids *F*, *G*, *H* and *I*. The face shared by the control volumes with centroids *P* and *F* is represented by *f* , with area *S<sup>f</sup>* and face unit normal vector **n***<sup>f</sup>* . The distance vector between *P* and *F* is represented by **d***PF*.

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:

$$\begin{split} &\rho \left( \int\_{V\_{P}} \frac{\partial u\_{i}}{\partial t} \, \mathrm{d}V\_{P} + \oint\_{S} n\_{j}(u\_{i}u\_{i}) \, \mathrm{d}S \right) + \oint\_{S} n\_{i}p \, \mathrm{d}S - \oint\_{S} n\_{j} \left[ (\eta\_{N}(T) + \eta^{\star}) \frac{\partial u\_{i}}{\partial x\_{j}} \right] \mathrm{d}S - \\ & \quad - \int\_{V\_{P}} \frac{\eta\_{E}(T)}{\lambda(T)} \frac{\partial (\mathrm{e}^{\mathbf{Y}} - \mathrm{I})\_{ij}}{\partial x\_{j}} \, \mathrm{d}V\_{P} = - \oint\_{S} n\_{j} \eta^{\star} \frac{\partial u\_{i}}{\partial x\_{j}} \, \mathrm{d}S, \end{split} \tag{15}$$

where *S* is the boundary of control volume *V<sup>P</sup>* and *n<sup>i</sup>* 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:

$$\begin{split} & V\_{P}\rho\_{P} \frac{(u\_{i})\_{P} - (u\_{i})\_{P}^{0}}{\Delta t} + \sum\_{f=\text{ub}(P)} (s\_{f}\rho u\_{f}\mu\_{i})\_{f} + \sum\_{f=\text{ub}(P)} (s\_{i}p)\_{f} - \\ & - \sum\_{f=\text{ub}(P)} \left[ s\_{f}(\eta\_{N}(T) + \eta^{\star}) \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} \right]\_{f} - \int\_{V\_{P}} \frac{\eta\_{E}(T)}{\lambda(T)} \frac{\partial (\mathbf{e}^{\mathbf{x}} - \mathbf{I})\_{ij}}{\partial \mathbf{x}\_{j}} \cdot \mathbf{d}V\_{P} \\ & = - \sum\_{f=\text{ub}(P)} \left( s\_{i}\eta^{\star} \frac{\partial u\_{i}}{\partial \mathbf{x}\_{j}} \right)\_{f} \end{split} \tag{16}$$

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<sup>i</sup>* 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

 *a uu P u<sup>P</sup>* + *a uv P v<sup>P</sup>* + *a up P p<sup>P</sup>* + *a u*Ψ*xx P* (Ψ*xx*)*<sup>P</sup>* + *a u*Ψ*xy P* (Ψ*xy*)*P*+ + ∑ *F*=*NB*(*P*) *a uu F u<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a uv F v<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a up F pF*+ + ∑ *F*=*NB*(*P*) *a u*Ψ*xx F* (Ψ*xx*)*<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a u*Ψ*xy F* (Ψ*xy*)*<sup>F</sup>* = *b u P* , *a vu P u<sup>P</sup>* + *a vv P v<sup>P</sup>* + *a vp P p<sup>P</sup>* + *a v*Ψ*xy P* (Ψ*xy*)*<sup>P</sup>* + *a v*Ψ*yy P* (Ψ*yy*)*P*+ + ∑ *F*=*NB*(*P*) *a vu F u<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a vv F v<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a vp F pF*+ ∑ *F*=*NB*(*P*) *a v*Ψ*xy F* (Ψ*xy*)*<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a v*Ψ*yy F* (Ψ*yy*)*<sup>F</sup>* = *b v P* , (17)

where *a uiφ P* and *a uiφ F* are the owner and neighbor coefficients in the discretized linear momentum equation representing the velocity component *u<sup>i</sup>* and the variable *φ* interactions, respectively; *b ui 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, *∂*(e <sup>Ψ</sup>−I)*ij ∂x<sup>j</sup>* , to the linear momentum equations, we will employ an implicit discretization by considering the following Taylor approximation for the functional (e <sup>Ψ</sup> <sup>−</sup> <sup>I</sup>)*ij* [28,39]

$$\begin{split} (\mathbf{e}^{\mathbf{v}} - \mathbf{I})\_{ij} &\approx (\mathbf{e}^{\mathbf{v}} - \mathbf{I})^{0}\_{ij} + \frac{\partial (\mathbf{e}^{\mathbf{v}} - \mathbf{I})\_{ij}}{\partial \mathbf{\varPsi}\_{kl}} \bigg|\_{\mathbf{\varPsi}^{0}\_{kl}} \left( \mathbf{\varPsi}\_{kl} - \mathbf{\varPsi}^{0}\_{kl} \right) \\ &= (\mathbf{e}^{\mathbf{v}} - \mathbf{I})^{0}\_{ij} + \sum\_{i,j}^{d} \mathbf{e}^{\lambda\_{m}/2} \mathbf{e}^{\lambda\_{l}/2} \frac{\sinh((\lambda\_{m} - \lambda\_{n})/2)}{(\lambda\_{m} - \lambda\_{n})/2} \sum\_{k,l} p^{m}\_{\mathbf{i\bar{k}}} p^{n}\_{lj} \left( \mathbf{\varPsi}\_{kl} - \mathbf{\varPsi}^{0}\_{kl} \right), \end{split} \tag{18}$$

where the derivative of the functional (e <sup>Ψ</sup> <sup>−</sup> <sup>I</sup>)*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 Ψ<sup>0</sup> *ij*. Therefore, the divergence of the functional (e <sup>Ψ</sup> <sup>−</sup> <sup>I</sup>)*ij* can be written as:

$$\begin{split} & \int\_{V\_{P}} \left( \frac{\eta\_{E}(T)}{\lambda(T)} \frac{\partial \left( \mathbf{e}^{\mathbf{Y}} - \mathbf{I} \right)\_{ij}}{\partial \mathbf{x}\_{j}} \right) \mathrm{d}V\_{P} \\ \approx & \frac{\eta\_{E}(T)}{\lambda(T)} \int\_{V\_{P}} \left( \frac{\partial \left( \mathbf{e}^{\mathbf{Y}} - \mathbf{I} \right)\_{ij}}{\partial \mathbf{x}\_{j}} \right)^{0} \mathrm{d}V\_{P} + \frac{\eta\_{E}(T)}{\lambda(T)} \int\_{V\_{P}} \partial \left[ \left( \frac{\partial \left( \mathbf{e}^{\mathbf{Y}} - \mathbf{I} \right)\_{ij}}{\partial \mathbf{Y}\_{kl}} \right)^{0} \mathbf{Y}\_{kl} \right] / \partial \mathbf{x}\_{j} \, \mathrm{d}V\_{P} \\ & - \frac{\eta\_{E}(T)}{\lambda(T)} \int\_{V\_{P}} \partial \left[ \left( \frac{\partial \left( \mathbf{e}^{\mathbf{Y}} - \mathbf{I} \right)\_{ij}}{\partial \mathbf{Y}\_{kl}} \right)^{0} \mathbf{Y}\_{kl}^{0} \right] / \partial \mathbf{x}\_{j} \, \mathrm{d}V\_{P}, \end{split} \tag{19}$$

Subsequently, applying the Gauss divergence theorem we can transform the volume integrals with derivatives into surface integrals

$$\int\_{V\_P} \partial \left[ \left( \frac{\partial (\mathbf{e}^\mathbf{v} - \mathbf{I})\_{ij}}{\partial \mathbf{\varPsi}\_{kl}} \right)^0 \mathbf{\varPsi}\_{kl} \right] / \partial \mathbf{x}\_j \, \mathbf{d}V\_P = \oint\_S n\_j \left( \frac{\partial (\mathbf{e}^\mathbf{v} - \mathbf{I})\_{ij}}{\partial \mathbf{\varPsi}\_{kl}} \right)^0 \mathbf{\varPsi}\_{kl} \, \mathbf{d}S,\tag{20a}$$

$$\int\_{V\_P} \partial \left[ \left( \frac{\partial \left( \mathbf{e}^\mathbf{F} - \mathbf{I} \right)\_{ij}}{\partial \mathbf{\varPsi}\_{kl}} \right)^0 \mathbf{\varPsi}\_{kl}^0 \right] / \mathbf{\varPsi}\_j \, \mathrm{d}V\_P = \oint\_S n\_j \left( \frac{\partial \left( \mathbf{e}^\mathbf{F} - \mathbf{I} \right)\_{ij}}{\partial \mathbf{\varPsi}\_{kl}} \right)^0 \mathbf{\varPsi}\_{kl}^0 \, \mathrm{d}S,\tag{20b}$$

and obtain the discretized form for the divergence of the functional (e <sup>Ψ</sup> <sup>−</sup> <sup>I</sup>)*ij* (i.e., the log-conformation tensor term) in the linear momentum equations as follows:

$$\begin{split} &\int\_{V\_{\mathbb{P}}} \left(\frac{\eta\_{E}(T)}{\lambda(T)} \frac{\partial \left(\mathbf{e}^{\mathbf{Y}} - \mathbf{I}\right)\_{ij}}{\partial \mathbf{x}\_{j}}\right) \mathrm{d}V\_{\mathbb{P}} \\ \approx &\frac{\eta\_{E}(T)}{\lambda(T)} \left(\frac{\partial \left(\mathbf{e}^{\mathbf{Y}} - \mathbf{I}\right)\_{ij}}{\partial \mathbf{x}\_{j}}\right)^{0} \mathrm{V}\_{\mathbb{P}} + \frac{\eta\_{E}(T)}{\lambda(T)} \sum\_{f=\mathrm{nb}(P)} \left(\mathbf{s}\_{j} \left(\frac{\partial \left(\mathbf{e}^{\mathbf{Y}} - \mathbf{I}\right)\_{ij}}{\partial \mathbf{Y}\_{kl}}\right)^{0} \mathbf{Y}\_{kl}\right)\_{f} \\ & - \frac{\eta\_{E}(T)}{\lambda(T)} \sum\_{f=\mathrm{nb}(P)} \left(\mathbf{s}\_{j} \left(\frac{\partial \left(\mathbf{e}^{\mathbf{Y}} - \mathbf{I}\right)\_{ij}}{\partial \mathbf{Y}\_{kl}}\right)^{0} \mathbf{Y}\_{kl}^{0}\right)\_{f} .\end{split} \tag{21}$$

Lastly, the contributions of the divergence of the log-conformation (DLC) tensor term for the linear momentum equations are given by

$$a\_{F, \text{DLC}}^{\
u\_{\dot{j}} \psi\_{kl}} = -\frac{\eta\_E(T)}{\lambda(T)} \left[ s\_{\dot{j}} \left( \frac{\partial \left( \mathbf{e}^{\mathbf{Y}} - \mathbf{I} \right)\_{\dot{i}\dot{j}}}{\partial \mathbf{\varPsi}\_{kl}} \right)^0 \right]\_f (1 - w\_f), \tag{22a}$$

$$\mathbf{a}\_{P,DLC}^{\boldsymbol{\mu}\_{\text{j}}\boldsymbol{\Psi}\_{kl}} = -\frac{\eta\_{E}(T)}{\lambda(T)} \sum\_{f=\text{nb}(P)} \left[ s\_{j} \left( \frac{\partial (\mathbf{e}^{\mathbf{y}} - \mathbf{I})\_{ij}}{\partial \mathbf{Y}\_{kl}} \right)^{0} \right]\_{f} w\_{f'} \tag{22b}$$

$$\boldsymbol{b}\_{\mathrm{P,DLC}}^{\boldsymbol{u}\_{\mathrm{j}}} = -\frac{\eta\_{\mathrm{E}}(T)}{\lambda(T)} \left[ \left( \frac{\partial (\mathbf{e}^{\mathbf{Y}} - \mathbf{I})\_{\mathrm{ij}}}{\partial \mathbf{x}\_{\mathrm{j}}} \right)^{0} \boldsymbol{V}\_{\mathrm{P}} - \sum\_{f=\mathrm{nb}(P)} \left( \mathbf{s}\_{j} \left( \frac{\partial (\mathbf{e}^{\mathbf{Y}} - \mathbf{I})\_{\mathrm{ij}}}{\partial \mathbf{Y}\_{\mathrm{kl}}} \right)^{0} \mathbf{Y}\_{\mathrm{kl}}^{0} \right)\_{f} \right], \tag{22c}$$

where *w<sup>f</sup>* 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.

#### *3.2. 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<sup>P</sup>* as follows:

$$\int\_{V\_P} \frac{\partial \mu\_l}{\partial \mathbf{x}\_l} \, \mathbf{d}V\_P = 0. \tag{23}$$

Subsequently, by employing the divergence theorem, the volume integral is transformed into a surface integral as follows:

$$\oint\_{\mathcal{S}} n\_i u\_i \,\mathrm{d}\mathcal{S} = 0.\tag{24}$$

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:

$$\sum\_{f=nb(P)} (s\_i u\_i)\_f = 0.\tag{25}$$

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

$$\sum\_{f=nb(P)} \left[ s\_i \left( -\frac{\partial p}{\partial \mathbf{x}\_i} \overline{D}\_i \right) \right]\_f + \sum\_{f=nb(P)} (s\_i \overline{u}\_i)\_f = \sum\_{f=nb(P)} \left[ s\_i \left( -\frac{\overline{\partial p}}{\partial \mathbf{x}\_i} \overline{D}\_i \right) \right]\_f. \tag{26}$$

Finally, we can write the mass conservation algebraic equation as:

$$a\_P^{pp} p\_P + a\_P^{pu} u\_P + a\_P^{pv} v\_P + \sum\_{F=NB(P)} a\_F^{pp} p\_F + \sum\_{F=NB(P)} a\_F^{pu} u\_F + \sum\_{F=NB(P)} a\_F^{pv} v\_F = b\_{P'}^p \tag{27}$$

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:

$$
\left[s\_i \left(-\frac{\partial p}{\partial \mathbf{x}\_i} \overline{D}\_i\right)\right]\_f = -\frac{\left[s\_i(s\_i \overline{D}\_i)\right]\_f}{(d\_{PF})\_i(s\_i)\_f}(p\_F - p\_P),\tag{28}
$$

where the coefficients of the implicit pressure gradient term for the mass conservation equation are given by

$$a\_F^{pp} = -\frac{\left[s\_i(s\_i\overline{D}\_i)\right]\_f}{(d\_{PF})\_i(s\_i)\_f} \tag{29a}$$

$$a\_P^{pp} = -\sum\_{F=NB(P)} a\_F^{pp} \,. \tag{29b}$$

The implicit coefficients for the second term in Equation (26) (corresponding to the velocity contribution) are given by

$$a\_F^{pu} = s\_f^x (1 - w\_f)\_\prime \qquad a\_F^{pv} = s\_f^y (1 - w\_f)\_\prime \tag{30a}$$

$$a\_P^{pu} = \sum\_{f=nb(P)} s\_f^x w\_{f\prime} \qquad a\_P^{pv} = \sum\_{f=nb(P)} s\_f^y w\_f \,. \tag{30b}$$

Lastly, the coefficients of the explicit pressure gradient term contribution for the mass conservation equation are given by

$$b\_P^p = \sum\_{f=nb(P)} \left[ s\_i \left( -\frac{\overline{\partial p}}{\partial \mathbf{x}\_i} \overline{D}\_i \right) \right]\_f \tag{31}$$

#### *3.3. 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])

$$\begin{split} \frac{\partial \Psi\_{ij}}{\partial t} &+ \mu\_k \frac{\partial \Psi\_{ij}}{\partial x\_k} + [\Psi\_{ij}, \Omega\_{ij}] - 2f(\text{ad}(\Psi\_{ij})) \mathbf{D}\_{ij} \\ &= -\frac{1}{\lambda(T)} \left( \mathbf{I}\_{ij} + \kappa \left( \mathbf{e}^{\mathbf{Y}\_{ij}} - \mathbf{I}\_{ij} \right) \right) \left( \mathbf{e}^{\mathbf{Y}\_{ij}} - \mathbf{I}\_{ij} \right) \mathbf{e}^{-\mathbf{Y}\_{ij}} \end{split} \tag{32}$$

where Ω*ij* = *∂u<sup>i</sup> ∂x<sup>j</sup>* − *∂u<sup>j</sup> ∂x<sup>i</sup>* /2 is the vorticity tensor and [Ψ*ij*, Ω*ij*] = Ψ*ij*Ω*ij* − Ω*ij*Ψ*ij* is the commutator term. Following Knechtges [38], *f*(ad(Ψ*ij*)) is defined by the application of the function *f*(*x*) = *<sup>x</sup>*/2 tanh(*x*/2) to the adjoint operator ad(Ψ*ij*), such as

$$f(\operatorname{ad}(\Psi\_{ij}))y = \sum\_{m,n=1}^{d} f(\lambda\_m - \lambda\_n) \mathbf{P}\_m y \mathbf{P}\_n. \tag{33}$$

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<sup>P</sup>* to yield

$$\begin{split} \int\_{V\_{P}} \frac{\partial \mathbf{V}\_{ij}}{\partial t} \, \mathrm{d}V\_{P} + \int\_{V\_{P}} u\_{k} \frac{\partial \mathbf{V}\_{ij}}{\partial \mathbf{x}\_{k}} \, \mathrm{d}V\_{P} + \int\_{V\_{P}} [\mathbf{\varPsi}\_{ij}, \mathrm{\varOmega}\_{ij}] \, \mathrm{d}V\_{P} - \int\_{V\_{P}} 2f(\mathrm{ad}(\mathbf{V}\_{ij})) \mathbf{D}\_{ij} \, \mathrm{d}V\_{P} &= \\ - \int\_{V\_{P}} \frac{1}{\lambda(T)} \Big( \mathrm{I}\_{ij} + \kappa \Big( \mathbf{e}^{\mathbf{V}\_{ij}} - \mathrm{I}\_{ij} \Big) \Big) \Big( \mathbf{e}^{\mathbf{V}\_{ij}} - \mathrm{I}\_{ij} \Big) \mathbf{e}^{-\mathbf{V}\_{ij}} \, \mathrm{d}V\_{P}. \end{split} \tag{34}$$

This leads to the algebraic equation with the following form:

 *a* Ψ*xx*Ψ*xx P* (Ψ*xx*)*<sup>P</sup>* + *a* Ψ*xxu P u<sup>P</sup>* + *a* Ψ*xxv P v<sup>P</sup>* + ∑ *F*=*NB*(*P*) *a* Ψ*xx*Ψ*xx F* (Ψ*xx*)*F*+ + ∑ *F*=*NB*(*P*) *a* Ψ*xxu F u<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a* Ψ*xxv F v<sup>F</sup>* = *b* Ψ*xx P* , *a* Ψ*xy*Ψ*xy P* (Ψ*xy*)*<sup>P</sup>* + *a* Ψ*xyu P u<sup>P</sup>* + *a* Ψ*xyv P v<sup>P</sup>* + ∑ *F*=*NB*(*P*) *a* Ψ*xy*Ψ*xy F* (Ψ*xy*)*F*+ + ∑ *F*=*NB*(*P*) *a* Ψ*xyu F u<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a* Ψ*xyv F v<sup>F</sup>* = *b* Ψ*xy P* , *a* Ψ*yy*Ψ*yy P* (Ψ*yy*)*<sup>P</sup>* + *a* Ψ*yyu P u<sup>P</sup>* + *a* Ψ*yyv P v<sup>P</sup>* + ∑ *F*=*NB*(*P*) *a* Ψ*yy*Ψ*yy F* (Ψ*yy*)*F*+ + ∑ *F*=*NB*(*P*) *a* Ψ*yyu F u<sup>F</sup>* + ∑ *F*=*NB*(*P*) *a* Ψ*yyv F v<sup>F</sup>* = *b* Ψ*yy P* , (35)

where *a* Ψ*ijφ P* and *a* Ψ*ijφ F* 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 *<sup>∂</sup>*Ψ*ij ∂t* and advective *u<sup>k</sup> ∂*(Ψ*ij*) *∂x<sup>k</sup>* 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:

$$\mathbb{E}\left[\boldsymbol{\Psi}\_{ij}, \boldsymbol{\Omega}\_{ij}\right] = \frac{1}{2} \left(\boldsymbol{\Psi}\_{ik}\frac{\partial u\_k}{\partial \boldsymbol{\chi}\_j} - \boldsymbol{\Psi}\_{ik}\frac{\partial u\_j}{\partial \boldsymbol{\chi}\_k} - \frac{\partial u\_i}{\partial \boldsymbol{\chi}\_k}\boldsymbol{\Psi}\_{kj} + \frac{\partial u\_k}{\partial \boldsymbol{\chi}\_i}\boldsymbol{\Psi}\_{kj}\right), \tag{36}$$

and, subsequently, we can use Taylor approximations such as

$$\begin{split} & -\frac{1}{2} \left( \boldsymbol{\Psi}\_{ik} \frac{\partial \boldsymbol{u}\_{j}}{\partial \mathbf{x}\_{k}} + \frac{\partial \boldsymbol{u}\_{i}}{\partial \mathbf{x}\_{k}} \boldsymbol{\Psi}\_{kj} \right) \\ & \approx +\frac{1}{2} \left[ \boldsymbol{\Psi}\_{ik}^{0} \left( \frac{\partial \boldsymbol{u}\_{j}}{\partial \mathbf{x}\_{k}} \right)^{0} + \left( \frac{\partial \boldsymbol{u}\_{i}}{\partial \mathbf{x}\_{k}} \right)^{0} \boldsymbol{\Psi}\_{kj}^{0} \right] - \frac{1}{2} \left[ \boldsymbol{\Psi}\_{ik} \left( \frac{\partial \boldsymbol{u}\_{j}}{\partial \mathbf{x}\_{k}} \right)^{0} + \left( \frac{\partial \boldsymbol{u}\_{i}}{\partial \mathbf{x}\_{k}} \right)^{0} \boldsymbol{\Psi}\_{kj} \right] \\ & -\frac{1}{2} \left[ \boldsymbol{\Psi}\_{ik}^{0} \frac{\partial \boldsymbol{u}\_{j}}{\partial \mathbf{x}\_{k}} + \frac{\partial \boldsymbol{u}\_{i}}{\partial \mathbf{x}\_{k}} \boldsymbol{\Psi}\_{kj}^{0} \right], \end{split} \tag{37}$$

and

$$\begin{split} &+\frac{1}{2}\left(\boldsymbol{\Psi}\_{ik}\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{j}}+\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{i}}\boldsymbol{\Psi}\_{kj}\right) \\ &\approx -\frac{1}{2}\left[\boldsymbol{\Psi}\_{ik}^{0}\left(\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{j}}\right)^{0}+\left(\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{i}}\right)^{0}\boldsymbol{\Psi}\_{kj}^{0}\right]+\frac{1}{2}\left[\boldsymbol{\Psi}\_{ik}\left(\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{j}}\right)^{0}+\left(\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{i}}\right)^{0}\boldsymbol{\Psi}\_{kj}\right] \\ &+\frac{1}{2}\left[\boldsymbol{\Psi}\_{ik}^{0}\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{j}}+\frac{\partial u\_{k}}{\partial \boldsymbol{x}\_{i}}\boldsymbol{\Psi}\_{kj}^{0}\right]. \end{split} \tag{38}$$

Starting with the terms in Equation (37), the negative commutator terms, the first contribution is explicit, being given by

$$(b\_P^{\Psi\_{ij}})\_{\text{neg}} = -\frac{1}{2} V\_P \left[ \Psi\_{ik}^0 \left(\frac{\partial u\_j}{\partial \mathbf{x}\_k}\right)^0 + \left(\frac{\partial u\_i}{\partial \mathbf{x}\_k}\right)^0 \Psi\_{kj}^0 \right]. \tag{39}$$

The second contribution is implicit in <sup>Ψ</sup>*ij* and explicit in *<sup>∂</sup>u<sup>i</sup> ∂x<sup>j</sup>* , being given by

$$(a\_P^{\Psi\_{ij}\Psi\_{ik}})\_{\text{neg},1} = -\frac{1}{2} V\_P \left(\frac{\partial u\_j}{\partial \mathbf{x}\_k}\right)^0,\tag{40a}$$

$$(a\_P^{\Psi\_{ij}\Psi\_{kj}})\_{\text{neg},2} = -\frac{1}{2}V\_P \left(\frac{\partial u\_i}{\partial \mathbf{x}\_k}\right)^0. \tag{40b}$$

The third contribution is implicit in *<sup>∂</sup>u<sup>i</sup> ∂x<sup>j</sup>* 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:

$$\begin{split} & \int\_{V\_{\mathbb{P}}} \left( \Psi\_{ik}^{0} \frac{\partial u\_{j}}{\partial \mathbf{x}\_{k}} + \frac{\partial u\_{i}}{\partial \mathbf{x}\_{k}} \Psi\_{kj}^{0} \right) \mathrm{d}V\_{\mathbb{P}} \\ &= \oint\_{\mathcal{S}} n\_{k} \left( \Psi\_{ik}^{0} u\_{j} + u\_{i} \Psi\_{kj}^{0} \right) \mathrm{d}\mathcal{S} - \int\_{V\_{\mathbb{P}}} \left[ \left( \frac{\partial \Psi\_{ik}}{\partial \mathbf{x}\_{k}} \right)^{0} u\_{j} + u\_{i} \left( \frac{\partial \Psi\_{kj}}{\partial \mathbf{x}\_{k}} \right)^{0} \right] \mathrm{d}V\_{\mathbb{P}}. \end{split} \tag{41}$$

Applying the Gauss divergence theorem, the discretized form of the terms on the right-hand side of Equation (41) is

$$\sum\_{f'=nb(P)} \left[ s\_k \left( \Psi\_{ik}^0 \boldsymbol{\mu}\_j + \boldsymbol{\mu}\_i \mathbf{Y}\_{kj}^0 \right) \right]\_f - V\_P \left[ \left( \frac{\partial \Psi\_{ik}}{\partial \mathbf{x}\_k} \right)^0 \boldsymbol{\mu}\_j + \boldsymbol{\mu}\_i \left( \frac{\partial \Psi\_{kj}}{\partial \mathbf{x}\_k} \right)^0 \right]. \tag{42}$$

Using linear weighted interpolation, we can write the contributions of the third term as:

$$(a\_F^{\Psi\_{ij}\mu\_j})\_{\text{reg.1}} = -\frac{1}{2}(s\_k \Psi\_{ik}^0)\_f (1 - w\_f)\_\text{.}\tag{43a}$$

$$(a\_P^{\Psi\_{ij} u\_j})\_{n \in \mathcal{g}, 1} = -\frac{1}{2} \left[ \sum\_{f = nb(P)} (s\_k \Psi\_{ik}^0)\_f w\_f - V\_P \left(\frac{\partial \Psi\_{ik}}{\partial \mathbf{x}\_k}\right)^0 \right],\tag{43b}$$

$$(a\_F^{\Psi\_{ij} u\_i})\_{\text{reg.2}} = -\frac{1}{2} (\Psi\_{kb}^0 \mathbf{s}\_k)\_f (1 - w\_f)\_\prime \tag{43c}$$

$$\langle a\_P^{\Psi\_{ij} u\_i} \rangle\_{n \text{g.2}} = -\frac{1}{2} \left[ \sum\_{f=nb(P)} (s\_k \Psi\_{kj}^0)\_f w\_f - V\_P \left(\frac{\partial \Psi\_{kj}}{\partial \mathbf{x}\_k}\right)^0 \right]. \tag{43d}$$

Following the same reasoning as given above, the terms in Equation (38), the positive commutator terms, generate the following coefficients:

$$(b\_P^{\Psi\_{ij}})\_{pos} = \frac{1}{2} V\_P \left[ \Psi\_{ik}^0 \left( \frac{\partial u\_k}{\partial \mathbf{x}\_j} \right)^0 + \left( \frac{\partial u\_k}{\partial \mathbf{x}\_i} \right)^0 \Psi\_{kj}^0 \right], \tag{44}$$

$$(a\_P^{\Psi\_{i\bar{j}}\Psi\_{i\bar{k}}})\_{pos,1} = \frac{1}{2} V\_P \left(\frac{\partial u\_k}{\partial x\_{\bar{j}}}\right)^0 \,\,\,\,\,\tag{45a}$$

$$(a\_P^{\Psi\_{ij}\Psi\_{kj}})\_{pos,2} = \frac{1}{2} V\_P \left(\frac{\partial u\_k}{\partial \mathbf{x}\_i}\right)^0 \,\,\,\,\,\tag{45b}$$

$$(a\_{\rm F}^{\Psi\_{ij} u\_k})\_{pos, 1} = \frac{1}{2} (s\_f \Psi\_{ik}^0)\_f (1 - w\_f)\_\prime \tag{46a}$$

$$(a\_P^{\Psi\_{ij} u\_k})\_{pos, 1} = \frac{1}{2} \left[ \sum\_{f = nb(P)} (s\_f \mathbf{\bar{Y}}\_{ik}^0)\_f w\_f - V\_P \left(\frac{\partial \mathbf{\bar{Y}}\_{ik}}{\partial x\_j}\right)^0 \right],\tag{46b}$$

$$(a\_F^{\Psi\_{i\bar{j}}\mu\_k})\_{pos,2} = \frac{1}{2}(s\_i\Psi\_{k\bar{j}}^0)\_f(1-w\_f)\_\prime \tag{46c}$$

$$(a\_P^{\Psi\_{ij} u\_k})\_{pos, 2} = \frac{1}{2} \left[ \sum\_{f = nb(P)} (s\_i \Psi\_{kj}^0)\_f w\_f - V\_P \left(\frac{\partial \Psi\_{kj}}{\partial x\_i}\right)^0 \right]. \tag{46d}$$

Lastly, we implicitly discretize the adjoint operator *f*(ad(Ψ))*ijkl*D*kl* by considering the following Taylor approximation [39]

$$\begin{split} & f(\operatorname{ad}(\Psi))\_{ijkl} \mathbf{D}\_{kl} \\ & \approx f(\operatorname{ad}(\Psi^{0}))\_{ijkl} \mathbf{D}\_{kl}^{0} + \sum\_{kl} \frac{\partial \Big( f(\operatorname{ad}(\Psi^{0}))\_{iqrj} \mathbf{D}\_{qr}^{0} \Big)}{\partial \Psi\_{kl}} (\Psi\_{kl} - \Psi\_{kl}^{0}) \\ & + \sum\_{kl} \frac{\partial \Big( f(\operatorname{ad}(\Psi^{0}))\_{iqrj} \mathbf{D}\_{qr}^{0} \Big)}{\partial \mathbf{D}\_{kl}} (\mathbf{D}\_{kl} - \mathbf{D}\_{kl}^{0}), \end{split} \tag{47}$$

where *<sup>∂</sup>*(*f*(ad(Ψ<sup>0</sup> ))*iqrj*D 0 *qr*) *∂*D*kl* is the derivative of the adjoint operator with respect to D*kl*, given by [38]

$$\begin{aligned} \frac{\partial \left( f(\text{ad}(\mathbf{Y}^0))\_{l\neq r} \mathbf{D}\_{qr} \right)}{\partial \mathbf{D}\_{kl}} \\ = f(\text{ad}(\mathbf{Y}^0))\_{l\neq kl} \\ = \sum\_{m,n=1}^d f(\lambda\_m - \lambda\_n) p\_{\text{ik}}^m p\_{l\neq r}^n \end{aligned} \tag{48}$$

and *<sup>∂</sup>*(*f*(ad(Ψ<sup>0</sup> ))*iqrj*D 0 *qr*) *∂*Ψ*kl* is the derivative of the adjoint operator with respect to Ψ*kl*, given by [49]

$$\begin{split} \frac{\partial \left( f(\text{ad}(\Psi^0))\_{ilq\gamma} \mathbf{D}\_{q\gamma}^0 \right)}{\partial \Psi\_{kl}} \\ = \sum\_{m,n,o=1}^d \frac{f(\lambda\_m - \lambda\_o) - f(\lambda\_n - \lambda\_o)}{\lambda\_m - \lambda\_n} \left( p\_{ik}^m p\_{lq}^n \mathbf{D}\_{q\gamma} p\_{rj}^o + p\_{iq}^o \mathbf{D}\_{q\gamma} p\_{rk}^n p\_{lj}^m \right). \end{split} \tag{49}$$

Note that, if *λ<sup>m</sup>* is equal to *λn*, then the denominator of Equation (49) needs to be replaced by *f* 0 (*λ<sup>m</sup>* − *λn*) [38].

Subsequently, we integrate the adjoint term over a control volume *V<sup>P</sup>* as follows:

$$\begin{split} &\int\_{V\_{\rm{P}}} 2f(\operatorname{ad}(\mathbf{Y}))\_{ijkl}\mathbf{D}\_{kl}\,\mathrm{d}V\_{\rm{P}} \\ &\approx \int\_{V\_{\rm{P}}} 2f(\operatorname{ad}(\mathbf{Y}^{0}))\_{ijkl}\mathbf{D}\_{kl}^{0}\,\mathrm{d}V\_{\rm{P}} + \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{i\dot{q}r}\mathbf{D}\_{qr})}{\partial\mathbf{Y}\_{kl}}\right]^{0}\Psi\_{kl}\,\mathrm{d}V\_{\rm{P}} \\ & \quad - \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{i\dot{q}r}\mathbf{D}\_{qr})}{\partial\mathbf{Y}\_{kl}}\right]^{0}\Psi\_{kl}^{0}\,\mathrm{d}V\_{\rm{P}} + \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{i\dot{q}r}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}\mathbf{D}\_{kl}\,\mathrm{d}V\_{\rm{P}} \\ & \quad - \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{i\dot{q}r}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}\mathbf{D}\_{kl}^{0}\,\mathrm{d}V\_{\rm{P}}. \end{split} \tag{50}$$

and substituting D*kl* = <sup>1</sup> 2 *∂u<sup>k</sup> ∂x<sup>l</sup>* + *∂u<sup>l</sup> ∂x<sup>k</sup>* in Equation (50), we obtain

$$\begin{split} &\int\_{V\_{P}} 2f(\operatorname{ad}(\mathbf{Y}))\_{ijkl} \mathbf{D}\_{kl} \operatorname{d}V\_{P} \\ &\approx \int\_{V\_{P}} 2f(\operatorname{ad}(\mathbf{Y}^{0}))\_{ijkl} \mathbf{D}\_{kl}^{0} \operatorname{d}V\_{P} + \int\_{V\_{P}} 2 \left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr})\mathbf{D}\_{qr}}{\partial \mathbf{Y}\_{kl}}\right]^{0} \mathbf{Y}\_{kl} \operatorname{d}V\_{P} \\ &\quad - \int\_{V\_{P}} 2 \left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr})\mathbf{D}\_{qr}}{\partial \mathbf{Y}\_{kl}}\right]^{0} \mathbf{Y}\_{kl}^{0} \operatorname{d}V\_{P} \\ &\quad + \int\_{V\_{P}} \left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr})\mathbf{D}\_{qr}}{\partial \mathbf{D}\_{kl}}\right]^{0} \left(\frac{\partial u\_{k}}{\partial \mathbf{x}\_{l}} + \frac{\partial u\_{l}}{\partial \mathbf{x}\_{k}}\right) \operatorname{d}V\_{P} \\ &\quad - \int\_{V\_{P}} 2 \left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr})\mathbf{D}\_{qr}}{\partial \mathbf{D}\_{kl}}\right]^{0} \mathbf{D}\_{kl}^{0} \operatorname{d}V\_{P}. \end{split} \tag{51}$$

As the fourth term on the right hand side of Equation (51) contains implicit velocity gradients, we employ integration by parts to linearize them, obtaining

$$\begin{split} &\int\_{V\_{\rm{P}}} 2f(\operatorname{ad}(\mathbf{F}))\_{ijkl}\mathbf{D}\_{kl}\operatorname{d}V\_{P} \\ &\approx \int\_{V\_{\rm{P}}} 2f(\operatorname{ad}(\mathbf{F}^{0}))\_{ijkl}\mathbf{D}\_{kl}^{0}\operatorname{d}V\_{P} + \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{F}))\_{ijqr}\mathbf{D}\_{qr})}{\partial\mathbf{V}\_{kl}}\right]^{0}\mathbf{Y}\_{kl}\operatorname{d}V\_{P} \\ &\quad - \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr}\mathbf{D}\_{qr})}{\partial\mathbf{Y}\_{kl}}\right]^{0}\mathbf{Y}\_{kl}^{0}\operatorname{d}V\_{P} \\ &\quad + \oint\_{S} \left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}(u\_{k}n\_{l} + u\_{l}n\_{k})\operatorname{dS} \\ &\quad - \int\_{V\_{\rm{P}}} \left[\frac{\partial\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}}{\partial\mathbf{x}\_{l}}u\_{k} + \frac{\partial\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr}\mathbf{D}\_{qr})}{\partial\mathbf{x}\_{k}}\right]^{0}}{\partial\mathbf{x}\_{k}}u\_{l}\right] \operatorname{d}V\_{P} \\ &\quad - \int\_{V\_{\rm{P}}} 2\left[\frac{\partial(f(\operatorname{ad}(\mathbf{Y}))\_{ijqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}\operatorname{\mathbf{D}^{0}\_{kl}}\operatorname{d}{$$

Thus, the discretized form of Equation (52) can then be written as:

$$\begin{split} &\int\_{V\_{P}} 2f(\operatorname{ad}(\mathbf{Y}))\_{ijkl}\mathbf{D}\_{kl}\operatorname{d}V\_{P} \\ \approx & 2f(\operatorname{ad}(\mathbf{Y}^{0}))\_{ijkl}\mathbf{D}\_{kl}^{0}V\_{P} + 2\left[\frac{\partial f(\operatorname{ad}(\mathbf{Y}))\_{ijr}\mathbf{D}\_{qr}}{\partial\mathbf{Y}\_{kl}}\right]^{0}\Psi\_{kl}V\_{P} \\ & -2\left[\frac{\partial (f(\operatorname{ad}(\mathbf{Y}))\_{iqr}\mathbf{D}\_{qr})}{\partial\mathbf{Y}\_{kl}}\right]^{0}\Psi\_{kl}^{0}V\_{P} + \sum\_{f=\operatorname{ad}(P)}\left[\left[\frac{\partial (f(\operatorname{ad}(\mathbf{Y}))\_{iqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}(\boldsymbol{u}\_{k\boldsymbol{l}}\boldsymbol{s} + \boldsymbol{u}\_{l}\boldsymbol{s}\_{k})\right]\_{f} \\ & -\left[\frac{\partial \left[\frac{\partial (f(\operatorname{ad}(\mathbf{Y}))\_{iqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}}{\partial\mathbf{x}\_{l}}\boldsymbol{u}\_{k} + \frac{\partial \left[\frac{\partial (f(\operatorname{ad}(\mathbf{Y}))\_{iqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}}{\partial\mathbf{x}\_{k}}\boldsymbol{u}\_{l}\right]V\_{P} \\ & -2\left[\frac{\partial (f(\operatorname{ad}(\mathbf{Y}))\_{iqr}\mathbf{D}\_{qr})}{\partial\mathbf{D}\_{kl}}\right]^{0}\mathbf{D}\_{kl}^{0}V\_{P}. \end{split} \tag{53}$$

The contributions of the adjoint term for the log-conformation tensor constitutive equations, using linear weighted interpolation, read as

$$\mathbb{E}\left(a\_P^{\mathbf{w}\_{ij}\mathbf{w}\_{kl}}\right)\_{adj} = 2\left[\frac{\partial \left(f(\mathbf{ad}(\mathbf{\varPsi}))\_{i\not\!p\varPi}\mathbf{D}\_{q\varGamma}\right)}{\partial \mathbf{\varPsi}\_{kl}}\right]^0 V\_{P\varGamma} \tag{54a}$$

$$(a\_F^{\mathbf{w}\_{ij}\mu\_k})\_{adj} = \left[ \left[ \frac{\partial \left( f(\text{ad}(\mathbf{\varPsi}))\_{ilpj} \mathbf{D}\_{qr} \right)}{\partial \mathbf{D}\_{kl}} \right]^0 s\_l \right]\_f (1 - w\_f)\_i \tag{54b}$$

$$\delta(a\_{P}^{\mathbf{V}\_{ij};\mu\_{k}})\_{\mathrm{adj}} = \sum\_{f=\mathrm{nb}(P)} \left[ \left[ \frac{\partial \left( f(\mathrm{ad}(\mathbf{V}))\_{\mathrm{i}\eta f} \mathbf{D}\_{\mathrm{q}T} \right)}{\partial \mathbf{D}\_{\mathrm{k}l}} \right]^{0} s\_{l} \right]\_{f} w\_{f} - \frac{\partial \left[ \frac{\partial \left( f(\mathrm{ad}(\mathbf{V}))\_{\mathrm{i}\eta f} \mathbf{D}\_{\mathrm{l}l} \right)}{\partial \mathbf{D}\_{\mathrm{l}l}} \right]^{0}}{\partial \mathbf{x}\_{l}} V\_{\mathrm{P}} \tag{54c}$$

$$(a\_F^{\mathbf{w}\_{ij}\mu\_l})\_{adj} = \left[ \left[ \frac{\partial (f(\text{ad}(\mathbf{Y}))\_{i\hat{q}r\hat{f}} \mathbf{D}\_{qr})}{\partial \mathbf{D}\_{kl}} \right]\_f^0 s\_k \right]\_f (1 - w\_f)\_i \tag{54d}$$

$$(a\_P^{\mathbf{F}\_{ijkl}})\_{adj} = \sum\_{f=\text{nb}(P)} \left[ \left[ \frac{\partial \left( f(\text{ad}(\mathbf{Y}))\_{lqij} \text{D}\_{qr} \right)}{\partial \text{D}\_{kl}} \right]\_f^0 s\_k \right]\_f \text{w}\_f - \frac{\partial \left[ \frac{\partial \left( f(\text{ad}(\mathbf{Y}))\_{lqij} \text{D}\_{qr} \right)}{\partial \text{D}\_{kl}} \right]^0}{\partial \text{x}\_k} \text{V}\_{P\prime} \quad \text{(54e)}$$

$$(b\_P^{\mathbf{F}\_{ij}})\_{adj} = 2f(\text{ad}(\mathbf{Y}^0))\_{ijkl} D\_{kl}^0 V\_P - 2 \left[ \frac{\partial \left( f(\text{ad}(\mathbf{Y}))\_{lqij} \text{D}\_{qr} \right)}{\partial \text{F}\_{kl}} \right]^0 \mathbf{Y}\_{kl}^0 V\_P$$

$$\left[ \text{d}(f(\text{ad}(\mathbf{Y}))\_{l:\cdots D\_r}) \right]^0$$

$$-2\left[\frac{\partial \left(f(\text{ad}(\Psi))\_{iqrj}\mathbf{D}\_{qr}\right)}{\partial \mathbf{D}\_{kl}}\right]^0 \mathbf{D}\_{kl}^0 V\_P. \tag{54f}$$

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.

#### *3.4. 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 *VP*, to yield

$$\begin{split} &\rho \mathbf{C}\_{p} \Big( \int\_{V\_{P}} \frac{\partial T}{\partial t} \, \mathbf{d}V\_{P} + \int\_{V\_{P}} u\_{i} \frac{\partial T}{\partial \mathbf{x}\_{i}} \, \mathbf{d}V\_{P} \Big) - \int\_{V\_{P}} k \frac{\partial^{2}T}{\partial \mathbf{x}\_{i}^{2}} \, \mathbf{d}V\_{P} \\ &= \int\_{V\_{P}} \Big[ (\tau\_{\mathbf{N}})\_{\dot{\boldsymbol{\eta}}} \mathbf{D}\_{\dot{\boldsymbol{\eta}}} + \frac{\eta\_{E}(T)}{\lambda(T)} \Big( a(\mathbf{e}^{\mathbf{y}} - \mathbf{I})\_{\dot{\boldsymbol{\eta}}} \mathbf{D}\_{\dot{\boldsymbol{\eta}}} + (1 - a) \frac{(\mathbf{e}^{\mathbf{y}} - \mathbf{I})\_{\dot{\boldsymbol{\eta}}}}{2\bar{\lambda}(T)} \Big) \Big] \, \mathbf{d}V\_{P}. \end{split} \tag{55}$$

Using the Gauss divergence theorem, the volume integrals of the advection and diffusion terms in Equation (55) are transformed into surface integrals as:

$$\begin{split} &\rho \mathbf{C}\_{p} \left( \int\_{V\_{P}} \frac{\partial T}{\partial t} \, \mathbf{d}V\_{P} + \oint\_{S} n\_{i} \boldsymbol{\mu}\_{i} T \, \mathbf{d}S \right) - \oint\_{S} n\_{i} \left( k \frac{\partial T}{\partial \mathbf{x}\_{i}} \right) \, \mathbf{d}S \\ &= \int\_{V\_{P}} \left[ (\mathbf{r}\_{N})\_{ij} \mathbf{D}\_{ji} + \frac{\eta\_{E}(T)}{\lambda(T)} \left( \boldsymbol{\mu}(\mathbf{e}^{\mathbf{v}} - \mathbf{I})\_{ij} \mathbf{D}\_{ji} + (1 - \boldsymbol{\alpha}) \frac{(\mathbf{e}^{\mathbf{v}} - \mathbf{I})\_{ii}}{2\bar{\lambda}(T)} \right) \right] \, \mathbf{d}V\_{P} . \end{split} \tag{56}$$

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

$$\begin{split} &V\_{P}\rho\_{P}(\mathbf{C}\_{P})\_{P}\frac{T\_{P}-T\_{P}^{0}}{\Delta t} + \sum\_{f=\text{nb}(P)} (s\_{i}\rho\mathbf{C}\_{p}u\_{i}T)\_{f} - \sum\_{f=\text{nb}(P)} \left(s\_{i}k\frac{\partial T}{\partial x\_{i}}\right)\_{f} \\ &= V\_{P}\left[ (\tau\_{N})\_{\stackrel{\text{ii}}{\text{j}}}\mathbf{D}\_{\stackrel{\text{iii}}{\text{j}}} + \frac{\eta\_{E}(T)}{\lambda(T)}\left(\alpha(\mathbf{e}^{\mathbf{Y}}-\mathbf{I})\_{\stackrel{\text{i}}{\text{j}}}\mathbf{D}\_{\stackrel{\text{j}}{\text{i}}} + (1-\alpha)\frac{(\mathbf{e}^{\mathbf{Y}}-\mathbf{I})\_{\stackrel{\text{ii}}{\text{j}}}}{2\overline{\lambda}(T)}\right)\right]. \end{split} \tag{57}$$

This leads to the algebraic equation for the energy balance with the following form:

$$a\_P^{TT}T\_P + a\_P^{Tu}u\_P + a\_P^{Tv}v\_P + \sum\_{F=NB(P)} a\_F^{TT}T\_F + \sum\_{F=NB(P)} a\_F^{Tu}u\_F + \sum\_{F=NB(P)} a\_F^{Tv}v\_F = b\_{P\prime}^{T} \tag{58}$$

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:

$$a\_{P, \text{clg}}^{TT} = \frac{V\_P \,\rho\_P \, (\text{C}\_p)\_P}{\Delta t} \,\text{},\tag{59a}$$

$$b\_{P, \text{cchg}}^T = \frac{V\_P \,\rho\_P \, (\text{C}\_p)\_P T\_P^0}{\Delta t}. \tag{59b}$$

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* (*m*˙ *<sup>f</sup>* = (*siρui*)*<sup>f</sup>* ) 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:

$$a\_{\rm F,adv}^{\rm T\mu} = a\_{\rm F,adv}^{\rm T\upsilon} = (\mathbb{C}\_p)\_f \max(\dot{m}\_f, \, 0), \tag{60a}$$

$$a\_{P,adv}^{Tu} = -\sum\_{F=NB(P)} a\_{F,adv}^{Tu} \qquad a\_{P,adv}^{Tv} = -\sum\_{F=NB(P)} a\_{F,adv}^{Tv} \tag{60b}$$

where the superscript *Tφ* means the influence of the *φ* variable in the *T* energy equation. The term max(*m*˙ *f* , 0) represents the maximum of *m*˙ *<sup>f</sup>* and 0, where the mass flux is positive if it goes from owner to neighbor cells, i.e., leaves the control-volume *VP*.

The implicit diffusion (*idi f f*) contribution, third term of Equation (57), is discretized, taking a linear profile (see page 86 of [48]) as

$$\left(s\_i k \frac{\partial T}{\partial \mathbf{x}\_i}\right)\_f = k\_f \frac{(s\_i \mathbf{s}\_i)\_f}{(d\_{PF})\_i (s\_i)\_f} (T\_\mathcal{F} - T\_\mathcal{P})\_\prime \tag{61}$$

where (*dPF*)*<sup>i</sup>* 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

$$a\_{\rm F,diff}^{TT} = -k\_f \frac{(s\_i s\_i)\_f}{(d\_{PF})\_i (s\_i)\_f} \,\prime \tag{62a}$$

$$a\_{P,diff}^{TT} = -\sum\_{F=NB(P)} a\_{F,diff}^{TT}.\tag{62b}$$

Finally, the coefficients of the explicit term contribution (right-hand side of Equation (57)) for the conservation of energy equation are given by

$$b\_{\rm P,source}^{\rm T} = V\_P \left[ (\tau\_N)\_{ij} \mathbf{D}\_{ji} + \frac{\eta\_E(T)}{\lambda(T)} \left( a(\mathbf{e}^{\mathbf{y}} - \mathbf{I})\_{ij} \mathbf{D}\_{ji} + (1 - a) \frac{(\mathbf{e}^{\mathbf{y}} - \mathbf{I})\_{ii}}{2\bar{\lambda}(T)} \right) \right]. \tag{63}$$

#### *3.5. 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:

 *a uu P a uv P a up P a u*Ψ*xx P a u*Ψ*xy P a u*Ψ*yy P a uT P a vu P a vv P a vp P a v*Ψ*xx P a v*Ψ*xy P a v*Ψ*yy P a vT P a pu P a pv P a pp P a p*Ψ*xx P a p*Ψ*xy P a p*Ψ*yy P a pT P a* Ψ*xxu P a* Ψ*xx v P a* Ψ*xx p P a* Ψ*xx*Ψ*xx P a* Ψ*xx*Ψ*xy P a* Ψ*xx*Ψ*yy P a* Ψ*xxT P a* Ψ*xyu P a* Ψ*xy v P a* Ψ*xy p P a* Ψ*xy*Ψ*xx P a* Ψ*xy*Ψ*xy P a* Ψ*xy*Ψ*yy P a* Ψ*xyT P a* Ψ*yyu P a* Ψ*yy v P a* Ψ*yy p P a* Ψ*yy*Ψ*xx P a* Ψ*yy*Ψ*xy P a* Ψ*yy*Ψ*yy P a* Ψ*yyT P a Tu P a Tv P a T p P a T*Ψ*xx P a T*Ψ*xy P a T*Ψ*yy P a TT P uP vP pP* (Ψ*xx*)*<sup>P</sup>* (Ψ*xy*)*<sup>P</sup>* (Ψ*yy*)*<sup>P</sup> TP* + + ∑ *F*=*nb*(*P*) *a uu F a uv F a up F a u*Ψ*xx F a u*Ψ*xy F a u*Ψ*yy F a uT F a vu F a vv F a vp F a v*Ψ*xx F a v*Ψ*xy F a v*Ψ*yy F a vT F a pu F a pv F a pp F a p*Ψ*xx F a p*Ψ*xy F a p*Ψ*yy F a pT F a* Ψ*xxu F a* Ψ*xx v F a* Ψ*xx p F a* Ψ*xx*Ψ*xx F a* Ψ*xx*Ψ*xy F a* Ψ*xx*Ψ*yy F a* Ψ*xxT F a* Ψ*xyu F a* Ψ*xy v F a* Ψ*xy p F a* Ψ*xy*Ψ*xx F a* Ψ*xy*Ψ*xy F a* Ψ*xy*Ψ*yy F a* Ψ*xyT F a* Ψ*yyu F a* Ψ*yy v F a* Ψ*yy p F a* Ψ*yy*Ψ*xx F a* Ψ*yy*Ψ*xy F a* Ψ*yy*Ψ*yy F a* Ψ*yyT F a Tu F a Tv F a T p F a T*Ψ*xx F a T*Ψ*xy F a T*Ψ*yy F a TT F uF vF pF* (Ψ*xx*)*<sup>F</sup>* (Ψ*xy*)*<sup>F</sup>* (Ψ*yy*)*<sup>F</sup> TF* = *b u P b v P b p P b* Ψ*xx P b* Ψ*xy P b* Ψ*yy P b T P* (64)

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<sup>i</sup>* , *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* , Ψ*<sup>n</sup> ij*, *T n* ).


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−<sup>6</sup> ), 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.

#### **4. 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].

#### *4.1. Geometry, Meshes, and Initial and Boundary Conditions*

An axisymmetric 4:1 sudden contraction with ratio of the radii *R*1/*R*<sup>2</sup> = 4 was chosen as test geometry (Figure 2a), because of the availability of numerical data in the literature [18]. The upstream length was *l*<sup>1</sup> = 80*R*<sup>2</sup> and the downstream length was *l*<sup>2</sup> = 50*R*2. The downstream channel height was chosen as *R*<sup>2</sup> = 0.0020604 m.

**Figure 2.** (**a**) Schematic of the 4:1 planar sudden contraction cross-section used to simulate the non-isothermal flow of a viscoelastic matrix fluid described by the Oldroyd-B constitutive model. The upstream length was *l*<sup>1</sup> and the downstream length was *l*2. The upstream channel height was *R*1 and the downstream channel height was *R*2. The temperature at the upstream wall is *Tw*,1, while, for the downstream wall, the temperature was *Tw*,2. (**b**) The geometrical domain was divided into five blocks for mesh discretization. (**c**) A detailed view of the contraction area with the minimum normalized cell size at the corner (∆*z*)min = (∆*r*)min = 0.01*R*2.

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<sup>z</sup>* = ∆*zi*+1/∆*z<sup>i</sup>* ,with ∆*z<sup>i</sup>* being the length of the cell *i* in the *z*-direction). In this way, since *f<sup>z</sup>* > 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.

**Table 1.** Characteristics of the meshes employed to simulate the non-isothermal viscoelastic flow of a viscoelastic matrix fluid described by the Oldroyd-B constitutive model in the axisymmetric 4:1 planar sudden contraction geometry.


Nz and Nr are number of cells along *z* and *r* directions, respectively, inside each block. *f<sup>z</sup>* and *f<sup>r</sup>* 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:


#### *4.2. 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:

$$Re = \frac{\rho R\_2 \overline{U}\_{z,2}}{\eta\_0} \tag{65}$$

$$
\dot{W}i = \frac{\lambda \overline{\mathcal{U}}\_{z,2}}{R\_2} \tag{66}
$$

$$Pe = \frac{\lambda R\_2 \overline{\mathcal{U}}\_{z,2} \mathcal{C}\_P}{k} \tag{67}$$

where *Uz*,2 is the mean velocity in the axial direction in the downstream channel, and *<sup>λ</sup>* is the relaxation time. In this case study, *Re* <sup>=</sup> 3.9 <sup>×</sup> <sup>10</sup>−<sup>5</sup> , corresponding to creeping flow conditions. The retardation ratio was equal to *β* = *ηS*/*η*<sup>0</sup> = 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<sup>P</sup>* = 1500 J/kg K and *k* = 0.17 W/mK. The WLF parameters were *C*<sup>1</sup> = 4.54 and *C*<sup>2</sup> = 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/*Uz*,2) of 10−<sup>4</sup> allowed for converged solutions to be obtained for all the runs performed. The maximum local Courant number corresponding to the normalized time-step 10−<sup>4</sup> obtained for the axisymmetric 4:1 sudden contraction is 0.07.

#### *4.3. Effects of the Energy Partitioning Parameter α*

In Figure 3, the temperature profile on the line *r* = 0.97*R*<sup>2</sup> (Figure 3a) and the temperature contour plots (Figure 3b) are shown as a function of the axial position *z*/*R*<sup>1</sup> 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 *Tw*,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.97*R*2. In both cases, we see the formation of a larger temperature-rising region for *z*/*R*<sup>1</sup> > 1, which is also extended through the downstream channel radial direction.

**Figure 3.** (**a**) Temperature profiles on the line *r* = 0.97*R*<sup>2</sup> [18] and (**b**) temperature contours as a function of the axial position *z*/*R*<sup>1</sup> for *α* = 0 and *α* = 1 at *Wi* = 5 and ∆*T* = 0 K using M2.

In Figure 4, the axial normal stress profile on the line *r* = 0.97*R*<sup>2</sup> (Figure 4a) and the axial normal stress contour plots (Figure 4b) as a function of the axial position *z*/*R*<sup>1</sup> 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*<sup>1</sup> < 0.2, which then increases in width, but with smaller normal stress values, at 0.2 < *z*/*R*<sup>1</sup> < 1.

**Figure 4.** (**a**) Axial normal stress profiles τ*P*,*zz*/(*η*0*Uz*,2/*R*2) on the line *r* = 0.97*R*<sup>2</sup> [18] and (**b**) axial normal stress contours as a function of the axial position *z*/*R*<sup>1</sup> for *α* = 0 and *α* = 1 at *Wi* = 5 and ∆*T* = 0 K using M2.

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*<sup>1</sup> = 0. Finally, at the contraction vertical wall, the shear stress is positive, allowing for the extension of the corner vortex size.

**Figure 5.** Contours of (**a**) first normal stress difference (τ*P*,*zz* − τ*P*,*rr*)/(*η*0*Uz*,2/*R*2) and (**b**) shear stress τ*P*,*rz*/(*η*0*Uz*,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.97*R*<sup>2</sup> (Figure 6a) and the axial velocity contour plots (Figure 6b) as a function of the axial position *z*/*R*<sup>1</sup> 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.

**Figure 6.** (**a**) Axial velocity profiles *Uz*/*Uz*,2 on the line *r* = 0.97*R*<sup>2</sup> [18] and (**b**) axial velocity contours as a function of the axial position *z*/*R*<sup>1</sup> for *α* = 0 and *α* = 1 at *Wi* = 5 and ∆*T* = 0 K using M2.

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.


M1: 4293 CV; M2: 17172 CV.

#### *4.4. 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*<sup>1</sup> = 0.25) temperature changes by more than 60 K and the centerline (*r*/*R*<sup>1</sup> = 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.

**Figure 7.** (**a**) Temperature *T* and (**b**) axial velocity *Uz*/*Uz*,2 as a function of the radial distance *r*/*R*<sup>1</sup> at the outlet, with *Wi* = 5, *α* = 0 and temperature jumps ∆*T* = −30 K, ∆*T* = 0 K and ∆*T* = +30 K using M2 [18].

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.

**Figure 8.** (**a**) Temperature *T* and (**b**) axial velocity *Uz*/*Uz*,2 contours at *Wi* = 5, *α* = 0 and temperature jumps ∆*T* = −30 K and ∆*T* = +30 K using M2.

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.

**Figure 9.** Contours of (**a**) first normal stress difference (τ*P*,*zz* − τ*P*,*rr*)/(*η*0*Uz*,2/*R*2) and (**b**) shear stress τ*P*,*rz*/(*η*0*Uz*,2/*R*2) for *α* = 0, *Wi* = 5, ∆*T* = −30 K (top) and ∆*T* = +30 K (bottom) using M2.

#### **5. 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.

**Funding:** This work is funded by FEDER funds through the COMPETE 2020 Programme and National Funds through FCT (Portuguese Foundation for Science and Technology) under the projects UID-B/05256/2020 and UID-P/05256/2020.

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to acknowledge the Minho University cluster under the project NORTE-07-0162-FEDER-000086 (URL: http://search6.di.uminho.pt (accessed on 1 February 2022)) for providing HPC resources that have contributed to the research results reported within this paper. The author thank M. Spahn from RWTH Aachen University for insightful comments regarding this work, specifically during the implementation of the isothermal log-conformation discretization schemes.

**Conflicts of Interest:** The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

