A numerical model based on the curvilinear coordinate system for the MAC method simplified

TL;DR

This paper introduces a simplified MAC numerical method using curvilinear coordinates and spline-based edge geometry to efficiently simulate incompressible fluid flows without free surfaces, validated through standard benchmark problems.

Contribution

The paper presents a novel simplified MAC method with curvilinear coordinates and spline-based geometry for improved fluid flow simulation accuracy and efficiency.

Findings

Accurate simulation of parallel plates, cavity, and atherosclerosis flows.

Results agree well with existing literature.

Demonstrates the method's applicability to complex geometries.

Abstract

In this paper we developed a numerical methodology to study some incompressible fluid flows without free surface, using the curvilinear coordinate system and whose edge geometry is constructed via parametrized spline. First, we discussed the representation of the Navier-Stokes and continuity equations on the curvilinear coordinate system, along with the auxiliary conditions. Then, we presented the numerical method -- a simplified version of MAC (\textit{Marker and Cell}) method -- along with the discretization of the governing equations, which is carried out using the finite differences method and the implementation of the FOU (\textit{First Order Upwind}) scheme. Finally, we applied the numerical methodology to the parallel plates problem, lid-driven cavity problem and atherosclerosis problem, and then we compare the results obtained with those presented in the literature. Keywords: finite differences, simplified MAC, curvilinear coordinates, parallel plates, did-driven cavity, atherosclerosis.

Tables

Table 1. Table 1: Meshes for simulations to first problem.

Mesh	Number of lines in	Number of lines in
	the direction $\xi$	the direction $\eta$
$P1$	9	5
$P2$	17	9
$P3$	33	17
$P4$	65	33
$P5$	129	65

Table 2. Table 2: Convergence speed of the numerical code via free stream V l n ​ u ​ m subscript subscript 𝑉 𝑙 𝑛 𝑢 𝑚 {V_{l}}_{num}

Meshes	$V_l^{}{}_{num}{}^{}$	$h_x$	$h_y$	$h_x\times h_y$	$Error=\left|1.5-V_l^{}{}_{num}{}^{}\right|$
$P_1$	1.1999	1.0	0.25	0.25	0.3001
$P_2$	1.3380	0.5	0.125	0.0625	0.1620
$P_3$	1.4158	0.25	0.0625	0.015625	0.0842
$P_4$	1.4498	0.125	0.03125	0.00390625	0.0502
$P_5$	1.5068	0.0625	0.015625	0.0009765625	0.0068

Table 3. Table 3: Localization of the main vortex center

Reference	$Re=100$	$Re=400$
This work	$(0.6109,0.7335)$	$(0.5699,0.6033)$
Bono et al. [2]	$(0.6157,0.7373)$	$(0.5613,0.6123)$
Ghia et al. [8]	$(0.6172,0.7344)$	$(0.5547,0.6055)$
Gupta e Kalita [10]	$(0.6125,0.7375)$	$(0.5500,0.6125)$
Hou et al. [11]	$(0.6196,0.7373)$	$(0.5608,0.6078)$
Marchi et al. [20]	$(0.6162,0.7373)$	$(0.5537,0.6054)$

Equations

$$\xi =\xi (x,y,t);\eta =\eta (x,y,t);\tau =\tau (t)$$

$$\frac{\partial \xi }{\partial x}=J\frac{\partial y}{\partial \eta };\frac{\partial \xi }{\partial y}=-J\frac{\partial x}{\partial \eta };\frac{\partial \eta }{\partial x}=-J\frac{\partial y}{\partial \xi };\frac{\partial \eta }{\partial y}=J\frac{\partial x}{\partial \xi }$$

$$J={\left(\frac{\partial x}{\partial \xi }\frac{\partial y}{\partial \eta }-\frac{\partial x}{\partial \eta }\frac{\partial y}{\partial \xi }\right)}^{-1}$$

$$\underset{\text{mbox}temporalterm}{\underbrace{\frac{\partial }{\partial \tau }\left(\frac{\rho u}{J}\right)}}$$

$$\underset{\text{mbox}temporalterm}{\underbrace{\frac{\partial }{\partial \tau }\left(\frac{\rho v}{J}\right)}}$$

$$\frac{\partial U}{\partial \xi }+\frac{\partial V}{\partial \eta }=0$$

$$\frac{\partial }{\partial \tau }\left(\frac{u}{J}\right)=-\mathcal{C}(u)+\frac{1}{\rho }{\mathcal{P}}^u+\nu \mathcal{V}(u);\frac{\partial }{\partial \tau }\left(\frac{v}{J}\right)=-\mathcal{C}(v)+\frac{1}{\rho }{\mathcal{P}}^v+\nu \mathcal{V}(v)$$

$${\frac{\partial }{\partial \tau }\left(\frac{u}{J}\right)|}_e^k=\frac{1}{J_e}\left(\frac{{u|}_e^{k+1}-{u|}_e^k}{\Delta \tau }\right)=-\mathcal{C}(u){|}_e^k+\frac{1}{\rho }{\mathcal{P}}^u{|}_e^{k+1}+\nu \mathcal{V}(u){|}_e^k$$

$${\frac{\partial }{\partial \tau }\left(\frac{v}{J}\right)|}_n^k=\frac{1}{J_n}\left(\frac{{v|}_n^{k+1}-{v|}_n^k}{\Delta \tau }\right)=-\mathcal{C}(v){|}_n^k+\frac{1}{\rho }{\mathcal{P}}^v{|}_n^{k+1}+\nu \mathcal{V}(v){|}_n^k$$

$$u{|}_e^{k+1}=F{|}_e^k+\frac{J_e\Delta \tau }{\rho }{\mathcal{P}}^u{|}_e^{k+1}$$

$$v{|}_n^{k+1}=G{|}_n^k+\frac{J_n\Delta \tau }{\rho }{\mathcal{P}}^v{|}_n^{k+1}$$

$${\mathcal{C}(u)|}_e^k\approx {\overline{U}}_E^k{u}_E^k-{\overline{U}}_P^k{u}_P^k+{\overline{V}}_{ne}^k{u}_{ne}^k-{\overline{V}}_{se}^k{u}_{se}^k$$

$${\mathcal{C}(v)|}_n^k\approx {\overline{U}}_{ne}^k{v}_{ne}^k-{\overline{U}}_{nw}^k{v}_{nw}^k+{\overline{V}}_N^k{v}_N^k-{\overline{V}}_P^k{v}_P^k$$

$${\left(\frac{1}{\rho }{\mathcal{P}}^u\right)|}_e^{k+1}\approx \frac{1}{\rho }\left[\left(p_{ne}^{k+1}-p_{se}^{k+1}\right){\frac{\partial y}{\partial \xi }|}_e-\left(p_E^{k+1}-p_P^{k+1}\right){\frac{\partial y}{\partial \eta }|}_e\right]$$

$${\left(\frac{1}{\rho }{\mathcal{P}}^v\right)|}_n^{k+1}\approx \frac{1}{\rho }\left[\left(p_{ne}^{k+1}-p_{nw}^{k+1}\right){\frac{\partial x}{\partial \eta }|}_n-\left(p_N^{k+1}-p_P^{k+1}\right){\frac{\partial x}{\partial \xi }|}_n\right]$$

$${\left(\nu \mathcal{V}(u)\right)|}_e^k={\nu }_e{\left[\frac{\partial }{\partial \xi }\left(J\left(\alpha \frac{\partial u}{\partial \xi }-\beta \frac{\partial u}{\partial \eta }\right)\right)\right]|}_e^k+{\nu }_e{\left[\frac{\partial }{\partial \eta }\left(J\left(\gamma \frac{\partial u}{\partial \eta }-\beta \frac{\partial u}{\partial \xi }\right)\right)\right]|}_e^k$$

$${\left(\nu \mathcal{V}(v)\right)|}_n^k={\nu }_n{\left[\frac{\partial }{\partial \xi }\left(J\left(\alpha \frac{\partial v}{\partial \xi }-\beta \frac{\partial v}{\partial \eta }\right)\right)\right]|}_n^k+{\nu }_n{\left[\frac{\partial }{\partial \eta }\left(J\left(\gamma \frac{\partial v}{\partial \eta }-\beta \frac{\partial v}{\partial \xi }\right)\right)\right]|}_n^k$$

$${\frac{\partial }{\partial \xi }\left(J\left(\alpha \frac{\partial u}{\partial \xi }-\beta \frac{\partial u}{\partial \eta }\right)\right)|}_e^k$$

$$-{(J\alpha )|}_P\left(u_e^k-u_w^k\right)$$

$$\frac{\partial }{\partial \eta }\left(J\left(\gamma \frac{\partial u}{\partial \eta }-\beta \frac{\partial u}{\partial \xi }\right)\right){|}_e^k$$

$$-(J\gamma ){|}_{se}^k\left(u_e^k-u_{sse}^k\right)$$

$$\frac{\partial }{\partial \xi }\left(J\left(\alpha \frac{\partial v}{\partial \xi }-\beta \frac{\partial v}{\partial \eta }\right)\right){|}_n^k$$

$$-(J\alpha ){|}_{nw}^k\left(v_n^k-v_{nww}^k\right)$$

$$\frac{\partial }{\partial \eta }\left(J\left(\gamma \frac{\partial v}{\partial \eta }-\beta \frac{\partial v}{\partial \xi }\right)\right){|}_n^k$$

$$-(J\gamma ){|}_P^k\left(v_n^k-v_s^k\right)$$

$$U{|}_e^{k+1}$$

$$U{|}_w^{k+1}$$

$$V{|}_n^{k+1}$$

$$V{|}_s^{k+1}$$

$${\frac{\partial U}{\partial \xi }|}_P^{k+1}+{\frac{\partial V}{\partial \eta }|}_P^{k+1}=0\Rightarrow U{|}_e^{k+1}-U{|}_w^{k+1}+V{|}_n^{k+1}-V{|}_s^{k+1}=0$$

$$J_e\left\{-{\frac{\partial p}{\partial \xi }|}_e^{k+1}\alpha {|}_e+{\frac{\partial p}{\partial \eta }|}_e^{k+1}\beta {|}_e\right\}+J_w\left\{{\frac{\partial p}{\partial \xi }|}_w^{k+1}\alpha {|}_w-{\frac{\partial p}{\partial \eta }|}_w^{k+1}\beta {|}_w\right\}$$

$$+J_n\left\{{\frac{\partial p}{\partial \xi }|}_n^{k+1}\beta {|}_n-{\frac{\partial p}{\partial \eta }|}_n^{k+1}\gamma {|}_n\right\}+J_s\left\{-{\frac{\partial p}{\partial \xi }|}_s^{k+1}\beta {|}_s+{\frac{\partial p}{\partial \eta }|}_s^{k+1}\gamma {|}_s\right\}$$

$$=\frac{\rho }{\Delta \tau }\{-F{|}_e^k{\frac{\partial y}{\partial \eta }|}_e+G{|}_e^k{\frac{\partial x}{\partial \eta }|}_e+F{|}_w^k{\frac{\partial y}{\partial \eta }|}_w-G{|}_w^k{\frac{\partial x}{\partial \eta }|}_w$$

$$+F{|}_n^k{\frac{\partial y}{\partial \xi }|}_n-G{|}_n^k{\frac{\partial x}{\partial \xi }|}_n-F{|}_s^k{\frac{\partial y}{\partial \xi }|}_s+G{|}_s^k{\frac{\partial x}{\partial \xi }|}_s\}$$

$$\text{mbox}(a)ve{l}_t=\frac{v{|}_{nww}+v{|}_n}{2}=0$$

$$ve{l}_n=u{|}_w=0$$

$$\text{mbox}(b)ve{l}_t=\frac{v{|}_n+v{|}_{nee}}{2}=0$$

$$ve{l}_n=u{|}_e=0$$

Full text

keyword: finite difference method, simplified MAC method, curvilinear coordinate system, parallel plates problem, lid-driven cavity problem, atherosclerosis problem

MSC[2010]: 35Q30, 76M20, 35A30

1 Introduction

In order to find numerical solutions to a differential equations it is indispensable to approximate the equations by a mathematical method that can be programmed in the computers. There are some techniques used for this purpose for example - finite differences [9], finite element [13] , finite volume [16], SPH (Smoothed Particles Hydrodynamics) [22], among others. We use finite differences in this work.

In general, the geometry which represent the physical domain of the problems are irregular forms. That is why we chose to study the curvilinear coordinate system. This system allows to generate mesh that coincides with the boundary of the domain in the applied problems [23, 26, 27]. In addition, we can to represent the differential equations, to be solved, in this system [33, 19].

When the partial differential equations are considered in the formulation of the mathematical model of a specific problem it is important to know the terms that compose these equations. The terms can be temporals, convectives, diffusives and others. The proper way to discretize such terms, respecting inherent physical of the problem, is necessary for the success of obtaining the numerical solution. In particular, the development of approximate methods of convective terms have been studied by many research in recent years [5]. The accuracy of the results obtained from the numerical solutions is directly influenced by the choice of convection scheme. The upwind scheme can be used in this purpose. This approach is done according with the sign of local convection speed. The upwind schemes are classified as first order or high order. In the first category we can mention the FOU scheme (First Order Upwinding) [4], which is stable unconditionally and produces a diffusive character which usually smooths the solution. Among the high-order we can mention the SOU (Second Order Upwinding) [25], QUICK (Quadratic Upstream Interpolation for Convective Kinematics) [15], which contributing to increase the accuracy of the numerical method, but introduce non-physical oscillations that can compromise the convergence. Finally, the third-order accurate upwind compact finite difference schemes, which allowed to obtain accurate numerical results for the benchmark flow problems [30].

In this work the objective is propose a versatile numerical method to study some incompressible fluids flow without surface free, in the curvilinear coordinate system, using the upwind scheme FOU to approximate the convective terms, with the board interpolated by parametrized spline, which allows better simulate a greater amount of complex problems. This numerical method is applied to the study of three problems of incompressible fluid flows: parallel plates problem, lid-driven cavity problem, and atherosclerosis problem. In addition, the work shows that our numerical methodology accurately reproduces the results observed in the literature of these three problems studied.

2 Curvilinear governing equations

In many cases, when we want to get analytical solutions to the problems under study, many simplifications are necessary, and this can depart us of the reality. When our goal is to study a certain physical phenomenon more realistically, we need initially to model the physics problem. As usually, the equations obtained in this modeling process has no analytical solution, numerical methods are then used to obtain the solution of the problem studied. Since it isn’t possible to obtain numerical solutions on the continuous region, as this involves infinite points, we first discretize the domain, that is, divide it in points and only in these is that we’ll find the solution of the problem. This set of points is called mesh. It is important that they are properly distributed so that the numerical solution represents satisfactorily the gradients of interest in the problem studied.

Generally, problems of everyday life are not evaluated in rectangular coordinates, because are represented by irregular geometries. Then, the curvilinear coordinate system it is more appropriate. Its main function is the representation of complex geometries, in these cases the cartesian coordinate system leads to a poor fit of the border, since the physical domain doesn’t match the domain mesh [33].

By means of a transformation (which may be numeric) between the cartesian coordinate system $(x,y)$ and the curvilinear coordinate system $(\xi,\eta)$, it is possible to map the domain written in the $(x,y)$ system to other written on a regular geometry $(\xi,\eta)$. The $(x,y)$ system is termed physical domain, while the $(\xi,\eta)$ is called computational domain.

For example, in the two-dimensional case, the computational domain is taken in rectangular form, regardless of the physical geometry, Figs. 1(a) and 1(b). Thus, the physical domain is transformed so the computational domain is always represented by a rectangle. For convenience it is assumed elementary volumes with unitary dimensions, i.e., $\Delta\xi=\Delta\eta=1$. So, as much as the coordinate lines take arbitrary spacing on the physical plane, in the computational dimensions are fixed.

The coordinates of an arbitrary point in the curvilinear coordinate system $(\xi,\eta)$ are related with the cartesian coordinate system $(x,y)$ for transformation equations as follows

[TABLE]

where $\tau=t$, because we are not admitting mesh movement. Therefore, the metric transformation are given by

[TABLE]

where

[TABLE]

is the jacobian of the transformation [33, 19].

The metrics enabled mapping of the physical domain for the computational, implying in the perform of the required geometric compensations. In this work we build the mesh as detailed in [3, 28] and the board was obtained via the parametrized Spline method. Moreover, mathematical compensations is also carried out, by chain rule, in the governing equations of the studied physical problem. Therefore, the conformity between the computational mesh and the governing equations allows adequately handle the computational numerical simulation.

So we can solve the problems of interest. First, it is need to label the points in the mesh such a way that calculations are performed correctly by the numerical method. Considering Fig. 3(b), we have that the mesh is composed by cells where the topological relationships are arranged and labeled as the cardinal points. The labels $P$, $E$, $W$, $N$, $S$, $NE$, $SE$, $NW$, $SW$ mean center, east, west, north, south, northeast, southeast, northwest, southwest, respectively. The abbreviations in tiny, positioned on the sides, are cardinal changes from the center of the cell labeled by $P$.

Assuming by hypothesis a laminar flow, newtonian, isothermal and incompressible, in two-dimensional, without the existence of the source terms, then the Navier-Stokes equations described in the curvilinear coordinate system are written by

[TABLE]

[TABLE]

where the viscosity $\mu$ is constant. For details of how to describe this equations in these coordinate system see [18, 32]. The terms $U$ and $V$ in (2) and (2) are named contravariant components of the velocity vector, these are normal to $\xi$ and $\eta$ lines respectively, and defined as $\displaystyle U=\frac{1}{J}\left(u\frac{\partial{\xi}}{\partial{x}}+v\frac{\partial{\xi}}{\partial{y}}\right)$ and $\displaystyle V=\frac{1}{J}\left(u\frac{\partial{\eta}}{\partial{x}}+v\frac{\partial{\eta}}{\partial{y}}\right)$. As the incompressibility hypothesis was assumed, the continuity equation is such that

[TABLE]

therefore the governing equations of our computational model, in curvilinear coordinates, are given by (2), (2) and (5).

The choice of auxiliary conditions is critical to the formulation of any problem described by differential equations [9, 6]. As the Navier-Stokes equations may be used to describe the flow in many situations, it is important to properly set the initial and boundary conditions of the problem under study. For the initial conditions is taken a velocity field that satisfies continuity equation. The boundary conditions used in the resolution of problems considered can be classified in no-slip and impermeability condition (CNEI), free-slip condition (CLES), prescribed injection condition (CIPR) and continuous ejection condition (CECO). The details about the numerical considerations made on the governing equations and auxiliary conditions are detailed as follows.

3 Numerical modeling

For the numerical solution of the fluid equations has been chosen the MAC (Marker and cell) method formulated by Harlow and Welch [1]. This method has the advantage of permitting the simulation of different types of flow, in the cartesian coordinate system, with and without free surface [21, 34] and multiphase flows [29].

In [24] the authors only describe the procedure of how to solve the mechanical equations of fluid via the curvilinear coordinates system. However, what innovates in this study is that we describe the numerical procedure considering: (1) the governing equations in dimensional form, (2) to an irregular geometry whose board is obtained by parametric Spline, and (3) the system of linear equations for the pressure resolution is solved by Gauss-Seidel methodology. In this work has been considered a MAC methodology to confined flow. With this action has been disregarded any moving marker particles deduction associated with the method, which results in simplification of the numerical calculations.

Furthermore it is considered displaced mesh, where the unknowns are stored in different positions. The displaced storage for the components of the velocity vector and pressure has a positive impact on numerical calculation, due to the fact that reduce numerical instability [6, 1]. As shown in Figs. 3(a) and 3(b), the pressure ($p$) is located in the center of the cell and the components of velocity vector ($u$ and $v$) in the centers of the faces. Note that $p$ refers to the pressure and $P$ the center of the cell. Rewriting the equations (2) and (2) as follows

[TABLE]

where $\displaystyle\frac{\partial{}}{\partial{\tau}}\left(\frac{u}{J}\right)$ and $\displaystyle\frac{\partial{}}{\partial{\tau}}\left(\frac{v}{J}\right)$ are the temporal terms, $\mathscr{C}(u)$ and $\mathscr{C}(v)$ the convective terms, $\mathscr{P}^{u}$ and $\mathscr{P}^{v}$ the pressure terms and $\mathscr{V}(u)$ and $\mathscr{V}(v)$ the diffusive terms and $\nu=\mu/\rho$ the viscosity term.

For the temporal term were done first-order approximation. In the case of convective terms apply approaches upwind type. The pressure terms are approximated by central differences. Finally, the diffusive terms are approximate by central differences, too. The discretization of the equations on the face $e$, for every cell within the domain, in the time level $k$, is written as

[TABLE]

analogously

[TABLE]

Setting $\displaystyle F|_{e}^{k}=J_{e}\Delta\tau\left[-\mathscr{C}(u)|_{e}^{k}+\nu\mathscr{V}(u)|_{e}^{k}\right]+u|_{e}^{k}$ the cartesian component $u$ of the velocity vector is expressed as

[TABLE]

and denoting $\displaystyle G|_{n}^{k}=J_{n}\Delta\tau\left[-\mathscr{C}(v)|_{n}^{k}+\nu\mathscr{V}(v)|_{n}^{k}\right]+v|_{n}^{k}$ we have

[TABLE]

One of the great difficulties in numerically solving the Navier-Stokes is about the discretization of the convective terms, because these terms are not linear. We use the FOU scheme for these terms. Of course this is not the most appropriate scheme to very low/high Reynolds number, but in this study we want to show that the numerical model is appropriate.

The discretization of the $\mathscr{C}(u)$ term on the face $e$ of the cell centered in the cardinal point $P$, in a time level $k$, using central differences, can be given as

[TABLE]

with convection velocities calculated by arithmetic mean, and applying FOU scheme to approximate $u$ component in the corresponding face. On the other hand, the discretization of the $\mathscr{C}(v)$ on the face $n$, in time level $k$, is given by

[TABLE]

with the convection velocities calculated by arithmetic mean, and approaches with the upwind scheme, too. For more details see [6]. The discretization of the pressure terms on the faces $e$ and $n$, in time level $k+1$, are as follows

[TABLE]

and

[TABLE]

Finally, the discretization of the diffusive terms in faces $e$ and $n$, in time level $k$, can be given as

[TABLE]

[TABLE]

with

[TABLE]

Since the contravariant components $U$ and $V$ can be rewritten as $\displaystyle U=u\frac{\partial{y}}{\partial{\eta}}-v\frac{\partial{x}}{\partial{\eta}}$ and $\displaystyle V=-u\frac{\partial{y}}{\partial{\xi}}+v\frac{\partial{x}}{\partial{\xi}}$, from the expressions (6) and (7) – and similar for the other faces – we get the compact form for contravariant components

[TABLE]

Approaching by central difference scheme in P cardinal point, in time level k + 1, the continuity equation gives us the following expression

[TABLE]

Replacing the equations (8) to (11) in (12) and grouping like terms in the same side of equality, we find

[TABLE]

The equation (3) is the pressure evolution equation, that satisfies the continuity equation (5).

The initial conditions considered should satisfy the continuity equation. We performed simulations always starting from a velocity field and pressure in a state of quiescence. We understand the quiescence state as the one with null velocity and pressure fields. The pressure boundary condition is taken as $\displaystyle\frac{\partial{p}}{\partial{n}}=0$. On the boundary conditions, whatever the two-dimensional problem under study, there are four kinds of configurations between the contour and interior cells in the computational domain. Denoting $vel_{t}$, $vel_{n}$ as the tangential and normal velocities at the border of the cells, $vel_{\bullet}$ (prescribed velocity), $vel_{I}$ (prescribed injection velocity), $vel_{E}$ (prescribed ejection velocity), then the four possible configurations are those shown in Fig. 4(d).

Case 1:

The non-slip and impermeability condition (CNEI) is defined by the following expressions $vel_{t}=0\mbox{ and }vel_{n}=0$, so

[TABLE]

Case 2:

The free-slip condition (CLES) is defined by the following expressions $vel_{t}=vel_{\bullet}\mbox{ and }vel_{n}=0$,

[TABLE]

Case 3:

The prescribed injection condition (CIPR) has the form $vel_{t}=0\mbox{ and }\linebreak vel_{n}=vel_{I}$

[TABLE]

Case 4:

The continuous ejection condition (CECO) is indicated by $\displaystyle\frac{\partial{}}{\partial{n}}vel_{t}=0\mbox{ and }\frac{\partial{}}{\partial{n}}vel_{n}=0$

[TABLE]

4 Numerical results

In this section we present the numerical results obtained from the proposed method in the previous sections. The first problem relates to the study of the flow between two parallel plates, the second deals with the flow in a square cavity with upper wall moving and the third refers to atherosclerosis. The mesh, for each problem studied, it was generated exactly as described in [3, 28]. The edges were built by Spline, and the mesh via numerical solution of the grid generation equations. For mesh generation has been stipulated a maximum of 1000 iterations to solve the linear system that create the grid, considering an error of less than $10^{-4}$ in this system. The numerical method of Gauss-Seidel was used to solve the linear system.

With the first problem we show that our code is able to obtain the numerical solution accurately when compared with analytical solution. The second case in our study aims at demonstrate that even for Reynolds number slightly high ($Re=1000$), where the term convective is dominant, it can properly address the problem. Besides this, with previous cases, show that the proposed numerical model simulates cases whose geometry is cartesian. Finally, the third case, detailing the simulation when the computational mesh is perfect adjusted the geometry of the problem. In this case, even with the dominant convective term ($Re=900$), the results are in accordance with literature too.

The color maps used in this study range from dark blue (lower velocity) to dark red (higher velocity), as shown in the Fig. 4 below. The results obtained for these problems are presented in the following.

4.1 Parallel plates problem

In this problem we consider a geometry in rectangular shape with height $H$ and length $L=8H$, according Fig. 5.

The fluid ($Re=100$) is injected into the geometry by edge $AD$, where the boundary condition is CIPR and $vel_{I}=1.0$ $m/s$. The plates are indicated by edges $AB$ and $CD$, where the CNEI conditions is applied. The fluid output occurs by the edge $BC$, where the boundary condition is CECO. For the study of the problem we performed five simulations, taking $H=1$ $m$ and varying the amount of lines considered in the directions $\xi$ and $\eta$, as shown in the table 1.

The initial condition was taken from the state of quiescence. The simulation was performed until the steady state was reached, and this occurred at $\tau=30$, approximately. In this study, for convergence of the pressure and velocity equations we consider the $\Delta\tau=10^{-2}$ to meshes $P1$ to $P3$. But to $P4$ and $P5$ the $\Delta\tau$ taken was equal to $5.10^{-3}$. This difference values in $\Delta\tau$ occurred because of the mesh refinement.

In Figs. 7(a) and 7(b) the velocity profiles in the output section, in the edge $BC$, are presented 111In this case we are only analyzing component $u$ of the velocity because this is perpendicular to the edge where there is fluid output. for the meshes considered in this work and in [2], respectively. Note that in [2] $D=1$, then $y/D=y$. Likewise, as $v_{e}=1$ follows $v_{l}/v_{e}=v_{l}=u$. Therefore, we can compare the graphs shown in the figures because both describe the same relationship.

Looking at the Fig. 7(a), can be observed that further refinement of the mesh gives better results. In the mesh $P5$, for example, the maximum value found for the speed was equal to $1.4814$ $m/s$. When comparing the profiles obtained by the proposed method in this work and literature [2], note that both are close to the analytical solution. However, in the literature the maximum value was 1.4756 m/s with the finite elements method. It is known in the scientific community which computational cost of finite differences (our proposal) used is lower than finite elements.

In Table 2 is shown that the mesh refinement implies lower error between the theoretical velocity and the calculated velocity in our code. So there is clearly a process of convergence to the theoretical velocity, i.e., the mesh doesn’t interfere in the solution obtained. Finally, the velocity field obtained, when the steady state has been reached, can be seen in the Fig. 7 below.

4.2 Lid-driven cavity problem

In this problem we consider a square cavity with edge measuring $1$ $m$, whose geometry is shown in Fig. 8. It is completely fluid filled, the upper wall is moving and we apply the boundary condition CLES with slip velocity equal to $1m/s$. The others walls are subject to boundary conditions CNEI.

The aim of this problem is the numerical verification of algorithm. The convective terms are more evident and the numerical model is more required numerically. The results obtained for this problem are compared with those presented by [2, 8, 9, 10, 11, 20].

This study was carried out from a mesh with 129 lines $\xi$ and $\eta$. The initial condition was taken from the state of quiescence. The simulation was performed until the steady state was reached, and this occurred at $\tau=50$, approximately. For convergence of the pressure and velocity equations, we consider the $\Delta\tau=10^{-3}$ because of the mesh refinement.

The variations of velocity fields obtained, after steady state has been reached, can be seen in the Figs. 10(a) and 10(b). In each of the figures we can see the formation of a primary vortex near the center of the domain, whose location varies with Reynolds number considered. Moreover, we can observe the formation of two other smaller vortices in the lower regions near to left and right boundaries of the cavity. In the first case, where $Re=100$, there is only an indication of the formation of secondary vortices. For $Re=400$ we observe a vortex formed in the lower right corner, and only one beginning of the vortex in the left corner. It is due to increased Reynolds number.

From the localization of the vortices, the interest was to compare the main vortex coordinates obtained with our methodology in relation to other studies. The table 3 present these coordinates for cases where $Re=100$ and $Re=400$. Note that our results are in agreement with the literature.

One third case ($Re=1000$), the secondary vortices are very defined in the lower right and left corners. There is evidence of the formation of another vortex in the upper left corner. Analogously, it is due to increased Reynolds number too.

Graphically, the results displayed in Fig. 11(a), compared with that obtained in Fig. 11(b), shows that our simulation is in accordance with literature [9]. The coordinates of the primary vortex obtained with our numerical code were (0.5528, 0.5698).

4.3 Atherosclerosis problem

Atherosclerosis is a disease associated with accumulation of lipids, complex carbohydrates, blood components, cells and other elements in large and medium-sized arteries, it is the main cause of heart disease [31]. In general, the development of this problem starts from accumulation of cholesterol LDL in the artery walls type, which may be higher or lower depending on the availability of this substance in the blood [17]. The accumulation of these compounds in the walls of an artery causes the same hardening through the formation of atherosclerotic plaques, which can lead to stenosis, or narrowing of the blood vessel, reducing blood flow in the artery [7].

Motivated by this problem, our objective was to reproduce a single case of blood flow in the region of a large caliber artery containing a stenosis. The constriction in the upper and lower border, both with same size, is located in the same point as presented in Fig. 11.

We consider a fluid with $Re=900$. This consideration is acceptable to the scientific community, and we can approach the blood as a Newtonian incompressible viscous fluid [14]. The fluid is injected into the geometry by edge $AD$, subject to the boundary condition CIPR with $vel_{I}=0.1467$. The output occurs through edge $BC$, where applied condition is CECO, and in the other walls we consider the boundary condition CNEI. As in previous cases, the initial condition was taken from the state of quiescence. For the study we consider the mesh shown in Fig. 12, constructed from the dimensions set out above, containing $\xi\times\eta=129\times 20$ lines, and $\Delta\tau=5.10^{-3}$.

The velocity field obtained after steady state is reached is shown in Fig. 13. The fluid enters in the geometry, suffers the performance of CNEI boundary condition, and the fully developed profile is established. When the fluid begins its entry path by stenosis undergoes a continuous drop in pressure since the Reynolds is moderated, thereby increasing the velocity. The maximum velocity value occurs in the narrowing region. After passing through narrowing the fluid velocity decreases and flows throughout the rest of the geometry, but in a lower level than that of the input. This pattern in the flow is also observed by [7], which leads to implications for human health.

In Fig. 14 are plotted the values of $u$ and $v$ components, respectively. From these graphs we can see that after the stenosis there is occurrence of vortex formation. Note that although the geometry has symmetry, the velocity components are asymmetric in the domain. The vortex near the top wall, in size and position with respect to the abscissa, is different from the vortex that is side of the bottom wall. This asymmetry in the velocity field is due to the moderate value of the Reynolds number used in the simulation. This asymmetric behavior was also obtained by [14]. For Reynolds values less than $900$ the flow tends to be symmetrical, because the vortices tend to disappear. But for values greater than $900$ the convective terms are strongly dominant, the flow tends to be asymmetric and vortices increase in size, intensity and position with respect to abscissa. The simulation here displayed when $Re=900$ is a turbulent presage.

The study of the disordered formation of these vortices beyond the stenosis is critical [14], once they can hinder the continuity of blood flow in arteries and thus, can aggravate health problems that may lead the patient to death.

5 Conclusion

The objective of this study was to present an effective numerical model to solve the equations of fluid dynamics in the case: laminar, newtonian, incompressible, isothermal, two-dimensional in curvilinear coordinates system. In addition, show that even applying a first order scheme (FOU) in the convective terms, it was possible to obtain satisfactory results with the numerical model.

In the cases of flows between parallel plates and in the lid-driven cavity, the numeric code was able to simulate the velocity field accurately, that is, the numerical model converged to values of interest does not depend on the mesh (see table 2) and also for various values of Reynolds (see Figs. 7(b), 10(b), 11(b)). Furthermore, we show that in the case of geometries whose cartesian coordinate system is completely inappropriate, the curvilinear coordinates system appears as a prominent alternative. This was the case of atherosclerosis problem. In this problem we show that the mesh is perfectly adequate to the domain, and that the proposed numerical model allowed to obtain numerical solutions in agreement with the literature.

For future work we intend to extend the methodology adopted in this study. First we will implement a higher-order convective scheme which solves problems in a more comprehensive range of Reynolds number, and then extend the code for three-dimensional cases.