1. Introduction

2. Governing Equations

• continuity equation:

  $$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0.$$

  (8)
• linear momentum equations:

  $$\begin{array}{cc}\hfill & \frac{\partial u}{\partial t}+\frac{\partial \left(uu\right)}{\partial x}+\frac{\partial \left(vu\right)}{\partial y}+\frac{\partial \left(wu\right)}{\partial z}=-\frac{\partial p}{\partial x}+\frac{1}{Re}\left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial z^2}\right]+\frac{\partial {\mathsf{\Phi }}^{xx}}{\partial x}+\frac{\partial {\mathsf{\Phi }}^{xy}}{\partial y}+\frac{\partial {\mathsf{\Phi }}^{xz}}{\partial z}+\frac{1}{F{r}^2}{g}_x,\hfill \\ \hfill & \frac{\partial v}{\partial t}+\frac{\partial \left(uv\right)}{\partial x}+\frac{\partial \left(vv\right)}{\partial y}+\frac{\partial \left(wv\right)}{\partial z}=-\frac{\partial p}{\partial y}+\frac{1}{Re}\left[\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}+\frac{{\partial }^{2}v}{\partial z^2}\right]+\frac{\partial {\mathsf{\Phi }}^{xy}}{\partial x}+\frac{\partial {\mathsf{\Phi }}^{yy}}{\partial y}+\frac{\partial {\mathsf{\Phi }}^{yz}}{\partial z}+\frac{1}{F{r}^2}{g}_y,\hfill \\ \hfill & \frac{\partial w}{\partial t}+\frac{\partial \left(uw\right)}{\partial x}+\frac{\partial \left(vw\right)}{\partial y}+\frac{\partial \left(ww\right)}{\partial z}=-\frac{\partial p}{\partial z}+\frac{1}{Re}\left[\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}+\frac{{\partial }^{2}w}{\partial z^2}\right]+\frac{\partial {\mathsf{\Phi }}^{xz}}{\partial x}+\frac{\partial {\mathsf{\Phi }}^{yz}}{\partial y}+\frac{\partial {\mathsf{\Phi }}^{zz}}{\partial z}+\frac{1}{F{r}^2}{g}_z,\hfill \end{array}$$

  (9)

  where $g_x,g_y,g_z$ are the Cartesian components of the gravity vector.
• stress tensor $\mathsf{\Phi }$:

  $$\begin{array}{cc}\hfill {\mathsf{\Phi }}^{xx}& ={\tau }^{xx}-\frac{2}{Re}\frac{\partial u}{\partial x},\hfill \\ \hfill {\mathsf{\Phi }}^{xy}& ={\tau }^{xy}-\frac{1}{Re}\left[\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right],\hfill \\ \hfill {\mathsf{\Phi }}^{xz}& ={\tau }^{xz}-\frac{1}{Re}\left[\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right],\hfill \\ \hfill {\mathsf{\Phi }}^{yy}& ={\tau }^{yy}-\frac{2}{Re}\frac{\partial v}{\partial y},\hfill \\ \hfill {\mathsf{\Phi }}^{yz}& ={\tau }^{yz}-\frac{1}{Re}\left[\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right],\hfill \\ \hfill {\mathsf{\Phi }}^{zz}& ={\tau }^{zz}-\frac{2}{Re}\frac{\partial w}{\partial z}.\hfill \end{array}$$

  (10)
• stress tensor τ:

  $$\begin{array}{cc}\hfill {\tau }^{xx}\left(t\right)& ={\int }_{-\infty }^t\sum _{k=1}^{m_1}\frac{a_k}{Wi{\lambda }_k}{e}^{-\frac{(t-t^{\prime })}{Wi{\lambda }_k}}\frac{\alpha }{\alpha -3+\beta I_1+(1-\beta )I_2}{B}_{t^{\prime }}^{xx}\left(t\right)d{t}^{\prime },\hfill \\ \hfill {\tau }^{xy}\left(t\right)& ={\int }_{-\infty }^t\sum _{k=1}^{m_1}\frac{a_k}{Wi{\lambda }_k}{e}^{-\frac{(t-t^{\prime })}{Wi{\lambda }_k}}\frac{\alpha }{\alpha -3+\beta I_1+(1-\beta )I_2}{B}_{t^{\prime }}^{xy}\left(t\right)d{t}^{\prime },\hfill \\ \hfill {\tau }^{xz}\left(t\right)& ={\int }_{-\infty }^t\sum _{k=1}^{m_1}\frac{a_k}{Wi{\lambda }_k}{e}^{-\frac{(t-t^{\prime })}{Wi{\lambda }_k}}\frac{\alpha }{\alpha -3+\beta I_1+(1-\beta )I_2}{B}_{t^{\prime }}^{xz}\left(t\right)d{t}^{\prime },\hfill \\ \hfill {\tau }^{yy}\left(t\right)& ={\int }_{-\infty }^t\sum _{k=1}^{m_1}\frac{a_k}{Wi{\lambda }_k}{e}^{-\frac{(t-t^{\prime })}{Wi{\lambda }_k}}\frac{\alpha }{\alpha -3+\beta I_1+(1-\beta )I_2}{B}_{t^{\prime }}^{yy}\left(t\right)d{t}^{\prime },\hfill \\ \hfill {\tau }^{yz}\left(t\right)& ={\int }_{-\infty }^t\sum _{k=1}^{m_1}\frac{a_k}{Wi{\lambda }_k}{e}^{-\frac{(t-t^{\prime })}{Wi{\lambda }_k}}\frac{\alpha }{\alpha -3+\beta I_1+(1-\beta )I_2}{B}_{t^{\prime }}^{yz}\left(t\right)d{t}^{\prime },\hfill \\ \hfill {\tau }^{zz}\left(t\right)& ={\int }_{-\infty }^t\sum _{k=1}^{m_1}\frac{a_k}{Wi{\lambda }_k}{e}^{-\frac{(t-t^{\prime })}{Wi{\lambda }_k}}\frac{\alpha }{\alpha -3+\beta I_1+(1-\beta )I_2}{B}_{t^{\prime }}^{zz}\left(t\right)d{t}^{\prime }.\hfill \end{array}$$

  (11)
• Finger tensor $\mathbf{B}$:

  $$\begin{array}{cc}\hfill \frac{\partial B^{xx}}{\partial t}+& \frac{\partial \left(u{B}^{xx}\right)}{\partial x}+\frac{\partial \left(v{B}^{xx}\right)}{\partial y}+\frac{\partial \left(w{B}^{xx}\right)}{\partial z}=2\left[\frac{\partial u}{\partial x}{B}^{xx}+\frac{\partial u}{\partial y}{B}^{xy}+\frac{\partial u}{\partial z}{B}^{xz}\right],\hfill \\ \hfill \frac{\partial B^{xy}}{\partial t}+& \frac{\partial \left(u{B}^{xy}\right)}{\partial x}+\frac{\partial \left(v{B}^{xy}\right)}{\partial y}+\frac{\partial \left(w{B}^{xy}\right)}{\partial z}=\frac{\partial v}{\partial x}{B}^{xx}+\left[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right]B^{xy}+\frac{\partial v}{\partial z}{B}^{xz}+\frac{\partial u}{\partial y}{B}^{yy}+\frac{\partial u}{\partial z}{B}^{yz},\hfill \\ \hfill \frac{\partial B^{xz}}{\partial t}+& \frac{\partial \left(u{B}^{xz}\right)}{\partial x}+\frac{\partial \left(v{B}^{xz}\right)}{\partial y}+\frac{\partial \left(w{B}^{xz}\right)}{\partial z}=\frac{\partial w}{\partial x}{B}^{xx}+\frac{\partial w}{\partial y}{B}^{xy}+\left[\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}\right]B^{xz}+\frac{\partial u}{\partial y}{B}^{yz}+\frac{\partial u}{\partial z}{B}^{zz},\hfill \\ \hfill \frac{\partial B^{yy}}{\partial t}+& \frac{\partial \left(u{B}^{yy}\right)}{\partial x}+\frac{\partial \left(v{B}^{yy}\right)}{\partial y}+\frac{\partial \left(w{B}^{yy}\right)}{\partial z}=2\left[\frac{\partial v}{\partial x}{B}^{xy}+\frac{\partial v}{\partial y}{B}^{yy}+\frac{\partial v}{\partial z}{B}^{yz}\right],\hfill \\ \hfill \frac{\partial B^{yz}}{\partial t}+& \frac{\partial \left(u{B}^{yz}\right)}{\partial x}+\frac{\partial \left(v{B}^{yz}\right)}{\partial y}+\frac{\partial \left(w{B}^{yz}\right)}{\partial z}=\frac{\partial w}{\partial x}{B}^{xy}+\frac{\partial v}{\partial x}{B}^{xz}+\frac{\partial w}{\partial y}{B}^{yy}+\left[\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right]B^{yz}+\frac{\partial v}{\partial z}{B}^{zz},\hfill \\ \hfill \frac{\partial B^{zz}}{\partial t}+& \frac{\partial \left(u{B}^{zz}\right)}{\partial x}+\frac{\partial \left(v{B}^{zz}\right)}{\partial y}+\frac{\partial \left(w{B}^{zz}\right)}{\partial z}=2\left[\frac{\partial w}{\partial x}{B}^{xz}+\frac{\partial w}{\partial y}{B}^{yz}+\frac{\partial w}{\partial z}{B}^{zz}\right].\hfill \end{array}$$

  (12)

3. Numerical Method

4. Semi-analytical solution

5. Results

5.1. Confined pipe flows

In this confined pipe flow (see Figure 3), the fluid is assumed to have only one relaxation mode. Therefore, it is possible to compare the simulation results with the semi-analytical solution presented before. The parameters used in this simulation are:

• Diameter $L=0.01\mathrm{m},U=0.025{\mathrm{m}.\mathrm{s}}^{-1},\rho =801.5{\mathrm{Kg}.\mathrm{m}}^{-3}$;
• ${\lambda }_{ref}=0.1396\mathrm{s},a_1=1.6648\mathrm{Pa},{\eta }_0=0.2324\mathrm{Pa}.\mathrm{s}$;
• Number of deformation fields $N=50$;
• $Re=\frac{\rho LU}{{\eta }_0}=0.8621,Wi={\lambda }_{ref}\frac{U}{L}=0.349$;
• Geometry: $0.01\mathrm{m}$ x $0,01\mathrm{m}$ x $0.05\mathrm{m}$;
• Meshes (number of cells in the x, y and z directions): $M1=12$x12x60 ($\delta x=\delta y=\frac{0.01}{12}$), $M2=16$x16x80 ($\delta x=\delta y=\frac{0.01}{16}$), $M3=20$x20x100 ($\delta x=\delta y=\frac{0.01}{20}$), $M4=24$x24x120 ($\delta x=\delta y=\frac{0.01}{24}$) and $M5=28$x28x140 ($\delta x=\delta y=\frac{0.01}{28}$).

Figure 4a and Figure 4b shows the velocity profiles w and u, respectively, with corresponding cross-section velocity distributions in the plane $xz$ (obtained mesh $M3$). The velocity profiles are fully developed at $\overline{t}=20\mathrm{s}$ (see also Figure 6a and Figure 7), suggesting that they have reached a steady-state condition. The velocity profile w shows a parabolic shape in the cross-section represented in Figure 4a, while the influence of the inflow can be observed in the cross-section depicted in Figure 4b, but the solution for v exhibits the expected physical behavior outside this region. Notice that, in all Full cells, the initial velocity vector is defined as $\mathbf{v}=(u,v,w)=(0,0,0)$.

Similarly to the velocity vector, initial conditions for pressure and tensors are set to zero in Full cells. As shown in Figure 5a and Figure 5b, the pressure and ${\tau }^{xz}$ tensor profiles are in agreement with the physically expected profiles, i.e., linear profiles across the longitudinal direction (flow direction) and transverse direction (perpendicular to the flow), respectively. In addition, the values obtained for ${\tau }^{xz}$ are comparable with the behavior of the analytical solution (see Figure 6b). In the cross-section represented in Figure 5b, there is an influence of the inflow (tensors are defined as zero in the inflow, and the finger tensor is defined as the identity matrix), but outside this region, the expected linear profile is obtained as previously stated (for further details, refer to Figure 6b).

Figure 5. Pressure and ${\tau }^{xz}$ distribution along the $xz$ plane ($y=0$).

Figure 5. Pressure and ${\tau }^{xz}$ distribution along the $xz$ plane ($y=0$).

Figure 6. Comparison between the analytical and numerical solutions for (a) w velocity component and tensors components (b) ${\tau }^{xz}$ and (c) ${\tau }^{zz}$.

Figure 6. Comparison between the analytical and numerical solutions for (a) w velocity component and tensors components (b) ${\tau }^{xz}$ and (c) ${\tau }^{zz}$.

Figure 6a-Figure 6c show, respectively, the solution for the velocity w and stress components ${\tau }^{xz}$ and ${\tau }^{zz}$, using meshes $M1$ - $M5$. The profiles were obtained at $y=0$ and were taken in the center of the pipe, $z=z_{max}/2$ in the time $t=10$ (or $\overline{t}=25\mathrm{s}$). The results were compared with the semi-analytical solution, with good agreement between the numerical results ($M1-M5$) and the semi-analytical results. The instant $t=10$ ($\overline{t}=25\mathrm{s}$) was chosen because the velocity and stresses already show a steady state behavior (the velocity residual is small (see Equation (42)). The residual (for the velocities) is calculated as

$$R_{es}={\displaystyle \frac{\sqrt{{\sum }_{i=1}^{N_c}\{{(u_i^t-u_i^{t-\delta t})}^2+{(v_i^t-v_i^{t-\delta t})}^2+{(w_i^t-w_i^{t-\delta t})}^2\}}}{N_c}},$$

(42)

where t is the simulation time, $\delta t$ is the time step, and $N_c$ is the number of computational cells.

Figure 7 shows the calculation of the residuals $R_{es}$ in meshes $M1-M5$ up to time $t=10$ (or equivalently, $\overline{t}=25\mathrm{s}$). It can be observed that the residuals in the five meshes $M1-M5$ decrease and show convergence towards a steady-state solution, thus proving the robustness of the numerical method. As expected, we also observe a smaller residual for the most refined meshes.

Figure 7. Total residual in meshes $M1$, $M2$, $M3$, $M4$ e $M5$ up to $t=10$ (or $\overline{t}=25\mathrm{s}$).

Figure 7. Total residual in meshes $M1$, $M2$, $M3$, $M4$ e $M5$ up to $t=10$ (or $\overline{t}=25\mathrm{s}$).