Three-Dimensional Simulations of Viscous Incompressible Fluid Flows on Grids with Unmatched Grid Interfaces

Abstract

The paper describes a numerical method that considers specific computational fluid dynamics (CFD) aspects of viscous incompressible flow simulations in the vicinity of interfaces between unmatched fragments of unstructured grids. The method is based on the general grid interface (GGI) principle, which involves conservative flux interpolation and does not require original grid modification at unmatched interfaces. A method is presented which combines adjacent unmatched fragments of an unstructured grid into a single domain by means of virtual interfaces considering connections between adjacent cells through virtual faces. The performance of the method is illustrated by the finite-volume discretization of the transport equation in the region of matched interfaces and its modification for the case of unmatched interfaces. The efficiency of the proposed method is demonstrated by three-dimensional CFD simulations with grid models composed of unmatched unstructured grid fragments. The simulation results are compared with equivalent simulations on matched grids. The influence of unmatched interfaces on the convergence rate and accuracy of the solution is assessed.

1. Introduction

At present, demand for CFD simulations of industrial objects based on the solution of the three-dimensional Navier–Stokes equations and involving large-scale models has grown significantly. Examples of such applications are simulations in nuclear power engineering, shipbuilding, and aircraft engineering, as well as other simulations that model hydrodynamic flows in complex-geometry structures [1,2,3,4,5,6]. In most cases, such simulations are done on a single grid that can contain more than a billion computational points (cells). The generation of such grid models representing all the details of physical processes taking place in the computational model is a nontrivial and generally time- and resource-intensive task. Dividing a complex initial geometry into simpler fragments makes the task of grid generation much easier, provides higher-quality grids, and efficiently saves computational resources. For example, optional grid refinement in individual fragments delivers a more accurate solution without significantly increasing the cell count in the grid model. Such grid models are usually composed of unmatched grid fragments generally having a different number and position of nodes at their adjacent boundaries. CFD simulations on such grids cannot be performed by standard approaches and require special methods to be developed. A key feature of CFD simulations using such grid models is that the unmatched grid fragments must be combined into a single computational model by means of special unmatched grid interfaces. The main problem of the unmatched grid interfaces is that they are nonconservative, which may lead to an unstable solution or perturbations [7].

Among the grid interface matching approaches proposed in the literature, there are a number of algorithms that geometrically merge fragments of arbitrary unstructured grids into a single grid model [8,9,10]. This is done by the modification of the original grid and the transformation of interfaces into a set of inner faces. These methods keep the grid model calculations conservative, but the process of the geometric merging of arbitrary unstructured grids requires using a large number of computational geometry algorithms (exact positioning of new nodes, node merging, projection of resulting nodes to edges and faces, etc.) and does not exclude the risk of creating domains with “poor” cells, the quality of which falls beyond admissible criteria.

There is a different group of methods, the so-called buffer/zipper layer [11,12], which is based on the creation in the region of unmatched grid interfaces of an intermediate unstructured layer to link the unmatched grid fragments. These methods are also conservative in grid model calculations. The drawback of this group of methods is that they cannot be used with grid models containing arbitrary unstructured grids, which limits the range of their applications in industry-oriented simulations on arbitrary-geometry grids.

In terms of computational efficiency, the most attractive approach is to use conservative methods based on the interpolation of quantities in the region of the unmatched interfaces. Such methods neither modify nor reconstruct the original grid at interfaces, which makes their programming much easier. One such method is the “patched grid” [13], which is widely used in calculations on multidomain grids. The drawback of this method is that it requires a high-quality structured grid.

One of the interpolation-based methods for combining unmatched unstructured grid interfaces is the general grid interface (GGI) method [14]. The GGI algorithm does not require modifying the original grid and is based on the use of weights determined by the face area ratio of virtual faces created by face projection from one of the interfaces to the other. The GGI method formed the basis for the development and successful implementation of a method enabling computations of nonoverlapping moving grids [15]. The authors of the GGI method demonstrated the conservativeness of the GGI interface and its modifications, given that the weights satisfy some criteria, but did not demonstrate that such weights can be found in all cases (for example, for unmatched unstructured interfaces). The authors also did not describe the procedure of interpolation between the faces of unmatched interfaces, which requires additional algorithms.

This paper describes a numerical method, which considers specific aspects of solving the Navier–Stokes equations in viscous incompressible flow simulations in the vicinity of interfaces between unmatched arbitrary unstructured grid fragments. The method is based on the philosophy of the GGI method. As distinct from the latter, our method does not require calculating the weights because it creates an identical set of virtual faces at adjacent unmatched interfaces, which guarantees its conservativeness. We propose a method of interpolation between the faces of unmatched interfaces based on area averaging. This method makes it possible to combine adjacent unmatched grid fragments into a single domain by linking adjacent cells through virtual faces, which saves time for problem-solving.

The numerical methods presented in this paper are implemented within the Russian software package LOGOS [16,17]. The performance of our method is demonstrated by simulations of turbulent flows in a circular pipe with a sudden restriction and in a circular diffuser. Each of the simulations is done using grid models consisting of matched and unmatched unstructured grid fragments. The influence of the unmatched interfaces on the convergence rate and accuracy of the solution is assessed.

2. Mathematical Model and Numerical Method

2.1. Mathematical Model

We consider a mathematical model of an isothermal flow of a viscous incompressible fluid. Under the given assumptions, the flow is described by the Navier–Stokes equations in the following form [18,19]:

$$\{\begin{array}{l}\nabla \cdot u=0,\\ \mathsf{\rho }\left(\frac{\partial u}{\partial t}+(u\cdot \nabla )u\right)-\nabla \cdot \mathsf{\tau }=-\nabla p.\end{array}$$

(1)

where ρ is the simulated fluid density, t is the time, u is the velocity vector of the averaged flow in three-dimensional Cartesian coordinates, τ = τ_μ + τ_t is the sum of the molecular and turbulent components in the viscous part of the stress tensor, p is the pressure, Γ is the diffusion coefficient, and φ is a scalar quantity.

The system of Equation (1) is solved by numerical integration on a finite-volume grid. As an illustration of the final-volume discretization of Equation (1), consider the discretization of the equation of transfer of scalar quantity φ:

$$\mathsf{\rho }\frac{\partial \mathsf{\phi }}{\partial t}+\mathsf{\rho }\nabla (u\mathsf{\phi })-\nabla (\mathsf{\Gamma }\nabla \mathsf{\phi })=0$$

(2)

The first term in (2) is an unsteady member, the second one is a convective member, and the next one is a diffusive member. To simplify our presentation below, we consider the diffusion coefficient Γ to be a constant. In the general case, the complete form for the discrete equivalent of the differential equation of transfer of the scalar quantity (2), considering the assumptions above, with the time member approximated by an implicit Eulerian scheme [20], for the reference volume (cell) P is given by:

$$\mathsf{\rho }\frac{{\mathsf{\phi }}_P-{\mathsf{\phi }}_P^{t-1}}{\Delta t}{V}_P+C_P-D_P=0.$$

(3)

Here, φ_P is the scalar quantity at the center of cell P at the current time step, ${\mathsf{\phi }}_P^{t-1}$ is the scalar quantity at the center of cell P at the previous time step, C_P is the convective member for cell P at the current time step, and D_P is the diffusive member for cell P at the current time step.

In what follows, the method is described for the case of discrete equivalents of the convective C_P and diffusive D_P members for a matched interface and their modification for an unmatched interface.

Consider a computational model composed of two grid fragments representing areas A and B, the adjacent boundaries of which fit each other in both the number and position of their nodes and create a pair of matched interfaces (Figure 1).

Approximated on the finite-volume grid, the convective and diffusive members at the current time step for cell P lying at the matched interface are given by

$$C_P={\displaystyle \sum _{f=face(P_{inner},P_{bound})}\mathsf{\rho }{\mathsf{\phi }}_f(u_f\cdot S_f)}+{\displaystyle \sum _{f=face(P_{interface})}\mathsf{\rho }{\mathsf{\phi }}_f(u_f\cdot S_f)},$$

(4)

$$D_P={\displaystyle \sum _{f=face(P_{inner},P_{bound})}\mathsf{\Gamma }(\nabla {\mathsf{\phi }}_f\cdot S_f)}+{\displaystyle \sum _{f=face(P_{interface})}\mathsf{\Gamma }(\nabla {\mathsf{\phi }}_f\cdot S_f)},$$

(5)

where S_f = S_f n_f is the area normal to face f, S_f is the area of face f, φ_f is the scalar quantity at the center of face f at the current time step, u_f is the velocity at the center of face f at the current time step, $\nabla$φ_f is the gradient of the scalar quantity on face f at the current time step, f = face(P_inner, P_bound) is summing over all the inner and boundary faces of cell P, and f = face(P_interface) is summing over the interface faces of cell P.

The values of unknowns φ_f and u_f for the inner and interface faces can be calculated by any of the known approximation schemes [18]. For the approximation scheme CD, the values of φ_f and u_f are determined on face f assuming the linear variation in φ_f and u_f between the centers of cells P and N (Figure 1), and are defined by

φ_f = λ_f φ_P + (1 − λ_f) φ_N,

(6)

u_f = λ_f u_P + (1 − λ_f) u_N,

(7)

where the value of the geometric interpolation coefficient λ_f is defined as the ratio of distances fN and PN projected on face normal n_f for cell P adjacent to cell N through face f (Figure 1), and is given by

$${\mathsf{\lambda }}_f=\frac{\left|n_f\cdot d^{Nf}\right|}{\left|n_f\cdot d^{Nf}\right|+\left|n_f\cdot d^{Pf}\right|},$$

(8)

where d^Nf and d^Pf are the vectors constructed between the cell centers and face center f (Figure 1).

The calculations of the gradients on the faces of arbitrary unstructured grids are corrected for nonorthogonality [21]. The value of gradient $\nabla {\mathsf{\phi }}_f^n$ on interface face f of cell P (Figure 1) at the n-th iteration of the current time step is defined by

$$\nabla {\mathsf{\phi }}_f^n=({\mathsf{\phi }}_N^n-{\mathsf{\phi }}_P^n)\frac{S_f}{S_f\cdot d_{PN}}+\overline{\nabla {\mathsf{\phi }}_{{}_f}^{n-1}}-(\overline{\nabla {\mathsf{\phi }}_f^{n-1}}\cdot d_{PN})\frac{S_f}{S_f\cdot d_{PN}},$$

(9)

where $\overline{\nabla {\mathsf{\phi }}_f^{n-1}}={\mathsf{\lambda }}_f^{}\nabla {\mathsf{\phi }}_P^{n-1}+(1-{\mathsf{\lambda }}_f^{})\nabla {\mathsf{\phi }}_N^{n-1}$ is the interpolated gradient of the scalar quantity on face f at the previous iteration.

In turn, gradient $\nabla$φ_P at the center of cell P can be calculated by the Green–Gauss algorithm [18] (the same is for cell N):

$$\nabla {\mathsf{\phi }}_P=\frac{1}{V_P}{\displaystyle \sum _{f=face(P)}{\mathsf{\phi }}_f{S}_f},$$

(10)

where f = face(P) is summing over all faces of cell P, including inner, boundary, and interface faces.

The discrete equivalent of the differential equation of transfer of the scalar quantity (3) can finally be represented as a system of algebraic equations, in which the following equation is constructed for each cell P in the domain (Figure 1):

$$a_P{\mathsf{\phi }}_P+a_{PN}{\mathsf{\phi }}_N+{\displaystyle \sum _{i=nb(P)}{a}_i{\mathsf{\phi }}_i=b_P},$$

(11)

where a_P is the diagonal coefficient of cell P, a_PN is the nondiagonal coefficient relating cells P and N through the matched interface, a_i are the nondiagonal coefficients relating cell P to the cells in domain A through common (inner) faces, b_P is the right-hand member, and i = nb(P) is summing over all neighbor cells i having common faces with cell P.

As a result, for the grid model composed of two matched domains A and B (Figure 2), a system of algebraic equations, (12), is developed, which can be solved by one of the iteration methods [22]:

$$\left(\begin{array}{cccccc}{a}_{A_1}& \cdots & 0& a_{A_1,B_1}& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & a_{A_N}& 0& \cdots & a_{A_N,B_N}\\ a_{A_1,B_1}& \cdots & 0& a_{B_1}& \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & a_{A_N,B_N}& 0& \cdots & a_{B_N}\end{array}\right)(\mathsf{\phi })=(b),$$

(12)

where $a_{A_i}$, i $\in$ [1…N] are the diagonal coefficients of the cells in area A, $a_{B_i}$ are the diagonal coefficients of the cells in area B, $a_{A_i{B}_i}$ are the nondiagonal coefficients relating cells A_i to cells B_i through the interface faces, (φ) is the vector of the sought scalar quantity, and (b) is the vector of the right-hand member.

To simplify the presentation, we do not show the nondiagonal coefficients relating the cells through the inner faces in Equation (12). Let us write the diagonal coefficients in Equation (12) for the cells in area A (the same is for the cells in area B) without considering the boundary conditions at the n-th iteration of the current time step:

$$a_{A_i}=\mathsf{\rho }\frac{{\mathsf{\phi }}_{A_i}^n}{\Delta t}{V}_{A_i}+S_{inner}^{A_i}+S_{interface}^{A_i},$$

(13)

where the term $S_{inner}^{A_i}$ is created due to the adjacency of cell A_i with the cells in area A through inner faces f:

$$S_{inner}^{A_i}={\displaystyle \sum _{f=face(A_{inner})}\left(\mathsf{\rho }{\mathsf{\lambda }}_f{u}_f^n\cdot S_f-\mathsf{\Gamma }\frac{S_f}{n_f\cdot d^{A_i{A}_j}}\right)},$$

(14)

where A_i is the cell adjacent to cell A_j through the inner face, the term $S_{interface}^{A_i}$ is created due to the adjacency of cell A_i with cell B_i through interface faces f:

$$S_{interface}^{A_i}={\displaystyle \sum _{f=face(A_{interface})}\left(\mathsf{\rho }{\mathsf{\lambda }}_f{u}_f^n\cdot S_f-\mathsf{\Gamma }\frac{S_f}{n_f\cdot d^{A_i{B}_i}}\right)}.$$

(15)

The nondiagonal coefficients in the system of linear algebraic Equation (12) relating the cells in area A (the same is for the cells in area B) through the inner faces at the n-th iteration of the current time step are given by

$$a_{A_i{A}_j}={\displaystyle \sum _{f=face(A_{inner})}\left(\mathsf{\rho }(1-{\mathsf{\lambda }}_f)u_f^n\cdot S_f+\mathsf{\Gamma }\frac{S_f}{n_f\cdot d^{A_i{A}_j}}\right)}.$$

(16)

The nondiagonal coefficients in the system of linear algebraic Equation (12) relating cell A_i to cell B_i through the interface faces at the n-th iteration of the current time step are given by

$$a_{A_i{B}_i}={\displaystyle \sum _{f=face(A_{interface})}\left(\mathsf{\rho }(1-{\mathsf{\lambda }}_f)u_f^n\cdot S_f+\mathsf{\Gamma }\frac{S_f}{n_f\cdot d^{A_i{B}_i}}\right)}.$$

(17)

The coefficients on the right side of the system of linear algebraic Equation (12) for the cells in area A (the same is for the cells in area B) without considering the boundary conditions at the n-th iteration of the current time step are given by

$$b_{A_i}=\mathsf{\rho }\frac{{\mathsf{\phi }}_{A_i}^{t-1}}{\Delta t}{V}_{A_i}+b_{A_i}^{inner}+b_{A_i}^{interface},$$

(18)

where $b_{A_i}^{inner}$ and $b_{A_i}^{interface}$ are given by

$$b_{A_i}^{inner}={\displaystyle \sum _{f=face(P_{inner})}\mathsf{\Gamma }{S}_f}{\left[\overline{\nabla {\mathsf{\phi }}_f}\cdot n_f-\frac{\overline{\nabla {\mathsf{\phi }}_f}\cdot d^{A_i{A}_j}}{n_f\cdot d^{A_i{A}_j}}\right]}^{n-1},$$

(19)

$$b_{A_i}^{interface}={\displaystyle \sum _{f=face(P_{interface})}\mathsf{\Gamma }{S}_f}{\left[\overline{\nabla {\mathsf{\phi }}_f}\cdot n_f-\frac{\overline{\nabla {\mathsf{\phi }}_f}\cdot d^{A_i{B}_i}}{n_f\cdot d^{A_i{B}_i}}\right]}^{n-1}.$$

(20)

The finite-volume discretization algorithm described above applies to all of the equations in (1). The method can be used in calculations of viscous incompressible flows on grid models of any dimensionality. However, the potential of its application is limited to calculations on models composed of a single grid or a grid composed of matched grid fragments. At present, this limitation makes the method much more difficult to use in most industry-oriented calculations, when grids are built of unmatched fragments to resolve details of specific physical processes. Such calculations require special interfaces, in which quantities φ_f, u_f in (4) and (5), and quantities $\nabla$φ_f in (9) on the faces belonging to unmatched interfaces cannot be calculated by standard approximation schemes, and the neighborhood of cell P with several cells on the opposite interface must be taken into account.

Below we describe a modification of the numerical method, which considers specific CFD aspects of viscous incompressible flow simulations in the vicinity of interfaces between unmatched grid fragments.

2.2. Unmatched Interface Linking Method

We consider a computational model composed of two grid fragments representing areas A and B, the adjacent boundaries of which create a pair of unmatched interfaces (Figure 3). The unmatched interfaces are understood to be adjacent boundaries of unmatched grid fragments, in which the number and position of nodes are different in the general case.

To merge the unmatched interfaces, the original faces are replaced by virtual ones. The geometric parameters of the virtual faces are determined by the successive projection of the original faces in area A to all of the original faces in area B. To identify the coordinates of the nodes created at the face intersections, one can use any of the polygon intersection algorithms [23]. The resulting set of virtual faces constitutes a virtual interface, which links the cells of the unmatched adjacent grid areas (Figure 3).

In the general case, the virtual interface consists of a set of inner and boundary virtual faces (Figure 4). The boundary virtual faces are not linked to the cells on the opposite interface and must be treated in accordance with the chosen boundary conditions. The inner virtual faces create a link between the adjacent cells on the interface.

Consider the modification of the diffusive (4) and convective (5) members with a transition from the original to the virtual faces for the case of cell P belonging to the unmatched interface (Figure 5).

A discrete equivalent of convective member (4) for cell P belonging to the unmatched interface is expressed as:

$$C_P={\displaystyle \sum _{f=face(P_{inner},P_{bound})}\mathsf{\rho }{\mathsf{\phi }}_f(u_f\cdot S_f)}+{\displaystyle \sum _{f=face(P_{interface})}{\displaystyle \sum _v\mathsf{\rho }{\mathsf{\phi }}_{f^v}(u_{f^v}\cdot S_{f^v})}},$$

(21)

where ${\mathsf{\phi }}_{f^v}$ and $u_{f^v}$ are the scalar quantity and velocity on virtual face f^v created on original face f at the current time step, and $S_{f^v}$ is the area of virtual face f^v normal to original face f.

A discrete equivalent of diffusive member (5) for cell P belonging to the unmatched interface is expressed as:

$$D_P={\displaystyle \sum _{f=face(P_{inner},P_{bound})}\mathsf{\Gamma }(\nabla {\mathsf{\phi }}_f\cdot S_f)}+{\displaystyle \sum _{f=face(P_{interface})}{\displaystyle \sum _v\mathsf{\Gamma }(\nabla {\mathsf{\phi }}_{f^v}\cdot S_{f^v})}},$$

(22)

where $\nabla {\mathsf{\phi }}_{f^v}$ is the gradient of the scalar quantity on virtual face f^v.

The values of unknowns ${\mathsf{\phi }}_{f^v}^n$, $u_{f^v}^n$ in (21) belong to the center of virtual face f^v (Figure 5) and can be calculated by any of the known schemes [20]. For the approximation scheme CD, the values of ${\mathsf{\phi }}_{f^v}^{}$ and $u_{f^v}^{}$ are determined on face f^v assuming the linear variation in ${\mathsf{\phi }}_{f^v}^{}$ and $u_{f^v}^{}$ between the centers of cells P and N_v (Figure 5) and are defined by

$${\mathsf{\phi }}_{f^v}={\mathsf{\lambda }}_{f^v}{\mathsf{\phi }}_P+(1-{\mathsf{\lambda }}_{f^v}){\mathsf{\phi }}_{N_v},$$

(23)

$$u_{f^v}={\mathsf{\lambda }}_{f^v}{u}_P+(1-{\mathsf{\lambda }}_{f^v})u_{N_v},$$

(24)

where the value of the geometric interpolation coefficient ${\mathsf{\lambda }}_{f^v}$ is defined as the ratio of distances f N_v and P N_v projected on face normal n_f for cell P adjacent to cell N_v through face f^v (Figure 5), and is given by

$${\mathsf{\lambda }}_{f^v}=\frac{\left|n_f\cdot d^{N_v{f}^v}\right|}{\left|n_f\cdot d^{N_v{f}^v}\right|+\left|n_f\cdot d^{P{f}^v}\right|},$$

(25)

where $d^{N_v{f}^v}$ and $d^{P_v{f}^v}$ are the vectors constructed between the cell centers and face center f (Figure 5).

Considering the nonorthogonality correction, the value of gradient $\nabla {\mathsf{\phi }}_{f^v}^n$ on virtual face f^v of cell P (Figure 1) at the n-th iteration of the current time step equals:

$$\nabla {\mathsf{\phi }}_{f^v}^n=({\mathsf{\phi }}_{N_v}^n-{\mathsf{\phi }}_P^n)\frac{S_{f^v}}{S_{f^v}\cdot d_{P{N}_v}}+\overline{\nabla {\mathsf{\phi }}_{f^v}^{n-1}}-(\overline{\nabla {\mathsf{\phi }}_{f^v}^{n-1}}\cdot d_{P{N}_v})\frac{S_{f^v}}{S_{f^v}\cdot d_{P{N}_v}},$$

(26)

where $\overline{\nabla {\mathsf{\phi }}_{f^v}^{n-1}}={\mathsf{\lambda }}_{f^v}\nabla {\mathsf{\phi }}_P^{n-1}+(1-{\mathsf{\lambda }}_{f^v})\nabla {\mathsf{\phi }}_{N_v}^{n-1}$ is the interpolated gradient of the scalar quantity on face f^v calculated at the previous iteration.

The value of gradient $\nabla$φ_P at the center of cell P belonging to the unmatched interface is calculated by the Green–Gauss algorithm (the same is for cells N_v):

$$\nabla {\mathsf{\phi }}_P=\frac{1}{V_P}\left({\displaystyle \sum _{f=face(P_{inner},P_{bound})}{\mathsf{\phi }}_f{S}_f}+{\displaystyle \sum _{f=face(P_{interface})}{\mathsf{\phi }}_f^{av}{S}_f}\right),$$

(27)

where ${\mathsf{\phi }}_f^{av}=\frac{{\displaystyle \sum _v{\mathsf{\phi }}_{f^v}{S}_{f^v}}}{{\displaystyle \sum _v{S}_{f^v}}}$ is the averaged value of the scalar quantity on interface face f.

As a result of the above transformations, additional terms incorporating the link between the adjacent cells through the virtual faces are included in the discrete equivalent of the equation of transfer (11) for each cell P belonging to the unmatched interface (Figure 5). In the general case, for the grid model composed of two unmatched areas A and B (Figure 6), a system of linear algebraic equations, (28), is developed, which can be solved by one of the iteration methods [22]:

$$\left(\begin{array}{cccccccc}{a}_{A_1}& \cdots & 0& a_{A_1,B_1}& a_{A_1,B_2}& \cdots & 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & a_{A_N}& 0& 0& \cdots & a_{A_N,B_{M-1}}& a_{A_N,B_M}\\ a_{A_1,B_1}& \cdots & 0& a_{B_1}& \cdots & \cdots & 0& 0\\ a_{A_1,B_2}& \cdots & 0& \cdots & a_{B_2}& \cdots & 0& 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ 0& \cdots & a_{A_N,B_{M-1}}& 0& 0& \cdots & a_{B_{M-1}}& 0\\ 0& \cdots & a_{A_N,B_M}& 0& 0& \cdots & 0& a_{B_M}\end{array}\right)(\mathsf{\phi })=(b),$$

(28)

where $a_{A_i}$, i $\in$ [1…N] are the diagonal coefficients of cells A_i, $a_{B_j}$, j $\in$ [1…M] are the diagonal coefficients of cells B_j, $a_{A_i{B}_j}$ are the nondiagonal coefficients relating cells A_i to cells B_j j through the virtual faces of the unmatched interface, (φ) is the vector of the sought scalar quantity, and (b) is the vector of the right-hand member.

To simplify the presentation, we do not show the nondiagonal coefficients relating the cells through the inner faces in the system of linear algebraic Equation (28). Let us write the diagonal coefficient in the system of linear algebraic Equation (28) for the cells in area A without considering the boundary conditions at the n-th iteration of the current time step (the same is for the cells in area B):

$$a_{A_i}=\mathsf{\rho }\frac{{\mathsf{\phi }}_{A_i}^n}{\Delta t}{V}_{A_i}+S_{A_i}^{inner}+S_{A_i}^{interface},$$

(29)

where $S_{A_i}^{inner}$ is the term created due to the adjacency of cell A_i with the cells in area A through the inner faces, just as in (14), and $S_{A_i}^{interface}$ is the term created due to the adjacency of cell A_i with cells B_i through the virtual faces of the unmatched interface:

$$S_{A_i}^{interface}={\displaystyle \sum _{f=face(A_{interface})}{\displaystyle \sum _v\left(\mathsf{\rho }{\mathsf{\lambda }}_{f^v}{u}_{f^v}^n\cdot S_{f^v}-\mathsf{\Gamma }\frac{S_{f^v}}{n_f\cdot d^{A_i{B}_i}}\right)}}.$$

(30)

Nondiagonal coefficients $a_{A_i{A}_j}$ in the system of linear algebraic Equation (28) relating the cells in area A through the inner faces at the n-th iteration of the current time step are determined identically to (16).

The nondiagonal coefficients in the system of linear algebraic Equation (28) relating cells A_i to cell B_j through the virtual faces of the unmatched interface at the n-th iteration of the current time step are defined by

$$a_{A_i{B}_j}={\displaystyle \sum _{f=face(A_{interface})}{\displaystyle \sum _v\left(\mathsf{\rho }(1-{\mathsf{\lambda }}_{f^v})u_{f^v}^n\cdot S_{f^v}+\mathsf{\Gamma }\frac{S_{f^v}}{n_f\cdot d^{A_i{B}_j}}\right)}},$$

(31)

where B_j, j $\in$ [1…M] are the cells bordering on cell A_i through virtual faces f^v.

The coefficients on the right side of the system of linear algebraic Equation (12) for the cells in area A (the same is for the cells in area B) without considering the boundary conditions at the n-th iteration are given by

$$b_{A_i}^n=\mathsf{\rho }\frac{{\mathsf{\phi }}_{A_i}^{t-1}}{\Delta t}{V}_{A_i}+b_{A_i}^{inner}+b_{A_i}^{interface},$$

(32)

where $b_{A_i}^{inner}$ is determined identically to (19), and $b_{A_i}^{interface}$ is given by

$$b_{A_i}^{interface}={{\displaystyle \sum _{f=face(P_{interface})}{\displaystyle \sum _v\left[\mathsf{\Gamma }{S}_{f^v}\left(\overline{\nabla {\mathsf{\phi }}_{f^v}}\cdot n_f-\frac{\overline{\nabla {\mathsf{\phi }}_{f^v}}\cdot d^{A_i{B}_j}}{n_f\cdot d^{A_i{B}_j}}\right)\right]}}}^{n-1}.$$

(33)

This method combines adjacent unmatched grid fragments into a single domain by means of created virtual interfaces and can be used for all equations in (1). It does not require the modification of the original grid and considers connections between adjacent cells through a set of virtual faces, which creates additional terms in the system of linear algebraic equations of the computational model. Thus, the calculations on models composed of unmatched computational domains by this method are equivalent to the calculations on matched computational models in terms of their computational expenses but are much less expensive in terms of the grid model construction.

In the section below, we illustrate the performance of our algorithm using simulations of a turbulent flow in a circular pipe with a sudden restriction, and in a round diffuser.

3. Numerical Experiments

Our method was implemented in the Russian software package LOGOS which is intended for three-dimensional heat and mass transfer simulations [16]. The performance of the method was demonstrated by numerical experiments which simulated widely known tests based on empirical calculations [24]. LOGOS had already been successfully verified by a series of test CFD calculations [25], so the numerical experiments were only aimed at comparing the results of calculations with identical models on the grids composed of matched and unmatched grid fragments.

3.1. Turbulent Incompressible Fluid Flow in a Circular Pipe with a Sudden Restriction

Let us consider how the method behaves in the numerical simulations of a turbulent incompressible fluid flow in a circular pipe with a sudden restriction [11]. The geometry of the computational model and the locations of control sections are shown in Figure 7. To test the method, the calculations were run on a matched and unmatched grid model.

The unmatched grid model consisted of two unmatched fragments constructed independently by the pre- and postprocessor of the LOGOS software package. The unmatched grid fragments in the restriction area are shown in Figure 8. The matched grid model was constructed by the geometric merging of the unmatched grid fragments in the LOGOS pre- and postprocessor.

The governing parameters in the computational models included a fluid density equal to 998.2 kg/m³ and a dynamic viscosity equal to 0.00101 kg/(m∙s). Turbulence was modeled by the RANS SST turbulence model [26] with a 0.05 percent turbulence intensity and a 0.001 m mixing length. The convective members in the calculations were approximated by the first-order upwind scheme UD [18].

Figure 9 shows the boundary conditions of the computational models. The outer boundaries were impermeable no-slip rigid walls. The inlet flow rate was V = 0.4227 m/s. The outlet pressure was P = 0 Pa. Under set parameters, the flow Reynolds number was 20925. To match the solution in the unmatched model, we used the interfaces implemented in LOGOS in accordance with the method described above (Figure 9). For the nonintersecting interface area, the impermeable no-slip wall boundary condition was used.

The calculations were performed by LOGOS in a steady setup on 80 parallel processors up to a convergence in mass in the order of 10⁻⁶ kg. The iteration convergence graph is shown in Figure 10.

To analyze the results, we compared the profiles of the velocity module (V) and pressure (P) along the horizontal line passing through the middle of the computational model (Figure 11), and the profiles of velocity and pressure at the control sections (Figure 12 and Figure 13).

The resulting profiles of velocity and pressure in the longitudinal and transverse sections demonstrated the close agreement between the results calculated by the matched and unmatched grid models. The convergence graphs indicate that the presence of the unmatched interface has a minor effect on the iteration convergence rate.

3.2. Turbulent Incompressible Fluid Flow in a Circular Diffuser

Let us consider how the method behaves in the numerical simulations of a turbulent flow in a circular pipe with a conical diffuser [24]. The geometry of the computational model and the locations of control sections are shown in Figure 14. To test the method, the calculations were run on a matched and unmatched grid model.

The unmatched grid model consisted of three unmatched areas (upstream pipe fragment, diffuser, and downstream pipe fragment). Grids in each of the areas were generated independently. A grid fragment in the diffuser area is shown in Figure 14. The matched grid model was constructed by an original mesh transformation in the LOGOS pre- and postprocessor to transform the unmatched interfaces into matched ones.

The governing parameters in the computational models included fluid density equal to 772.9 kg/m³ and a dynamic viscosity equal to 0.002682 kg/(m∙s). Under set parameters, the flow Reynolds number was 49989. Turbulent mixing was modeled by the SST turbulence model [26] with a 0.05 percent turbulence intensity and a 0.001 m mixing length. The convective members in the calculations were approximated by the first-order upwind scheme UD [18].

Figure 15 shows the boundary conditions of the computational models. The outer boundaries were impermeable no-slip rigid walls. The inlet velocity was constant and controlled by a constant flow rate of Q = 3.159662 kg/s. To match the solution in the unmatched model, we used the interfaces implemented in LOGOS in accordance with the method described above.

The calculations were performed by LOGOS in a steady setup on 80 parallel processors up to a convergence in mass in the order of 10⁻⁶ kg. The iteration convergence graph is shown in Figure 16.

To analyze the results, we compared the profiles of the velocity module and pressure along the horizontal line passing through the middle of the computational model (Figure 17), and the profiles of velocity and pressure at the control sections (Figure 18 and Figure 19).

The resulting profiles of velocity and pressure in the longitudinal and transverse sections demonstrated a close agreement between the results calculated by the matched and the unmatched grid models. The convergence graphs indicate that the presence of the unmatched interface has a minor effect on the iteration convergence rate.

4. Conclusions

In the paper, we described a numerical method, which considers specific aspects of solving the Navier–Stokes equations in viscous incompressible flow simulations in the vicinity of interfaces between unmatched unstructured grid fragments. The method combines adjacent unmatched fragments of an unstructured grid into a single domain by means of virtual interfaces considering the connections between adjacent cells through virtual faces. The algorithm of this method was implemented in the software package LOGOS within its CFD solver, which was used to run test calculations. Results of the test calculations by this method on grid models with matched and unmatched interfaces were presented. The presence of the unmatched interfaces was shown to have no significant effect on the flow pattern and solution convergence rate.