Impact of Artificial Compressibility on the Numerical Solution of Incompressible Nanofluid Flow

Nomenclature

1  Introduction

Heat transfer problems commonly occur in many industries, most notably in the strategic oil and gas and petrochemical industries. In these problems, heat is mainly transferred via convection, conduction, and radiation [1], with convectional heat transfer between fluids and solids being particularly important. Flow inside the cavity has been widely used as a case study to evaluate the performance of new numerical methods [2]. Upwind numerical methods for compressible flows have been widely developed. These methods determine the convective terms at the boundary of the cells by analyzing the propagation of sound waves. For incompressible flows, the speed of sound is infinite, hence these methods cannot be used. On the other hand, in solving the governing equations of incompressible flows, the absence of the pressure term in the continuity equation is problematic. The artificial compressibility method can simultaneously resolve each of these issues. This method changes the type of equations by adding a pressure gradient term to the continuity equation. Adding this term creates a kind of artificial compressibility that allows the definition of virtual sound speed. Therefore, it is possible to apply the various numerical methods used for compressible flows. These benefits motivate the authors to use this scheme. There is no specific restriction for using this method. To determine the pressure at each step of the numerical solution, the Navier-Stokes equations are coupled via an expression connected to the continuity equation, and the governing equations are transformed from elliptical to hyperbolic [3]. The added expression contains a pressure gradient to time and is referred to as the artificial compressibility term. The effect of the artificial compressibility term decreases throughout the numerical solution, and eventually, when the steady-state solution is reached, this term should equal zero, as the gradient of each parameter to time is zero in a steady-state condition. By introducing artificial compressibility, the propagation speed of virtual waves can be determined. In fact, the propagation speed of the waves in an ideal incompressible flow is infinite [4]. The artificial compressibility method, on the other hand, limits the wave propagation speed to a value given by the compressibility coefficient function. The behavior of the boundary layer in viscous fluids is highly sensitive to pressure changes in the flow direction, especially at the separationpoint [5]. If the flow separation occurs, the propagating pressure wave at a finite velocity can change the local pressure gradient and affect the separation point [6]. This phenomenon delays convergence to a steady-state.

Artificial compressibility methods have widely been used by several researchers. Qin et al. [7] compared the SIMPLE and artificial compressibility methods based on their convergence speed and accuracy. Two-dimensional incompressible cavity flow was simulated using these methods in the finite volume framework. The results revealed that the proposed method has a high degree of accuracy and rapid convergence. The optimum value of the artificial compressibility factor was not surveyed by the authors. Ntouras et al. [8] used artificial compressibility to simulate the free surface flow. The simulation was performed to demonstrate the ability of artificial compressibility on unsteady flow with a free surface. The authors showed that the governing equations have hyperbolic nature in a pseudo-time state. When the results of the artificial compressibility approach were compared to experimental and analytical solutions, they showed good accuracy. A constant artificial compressibility factor was used in their study. Some researchers combined artificial compressibility with other methods. For instance, Manzanero et al. [9] combined the artificial compressibility with the discontinuous Galerkin spectral element method to solve incompressible Navier-Stokes equations. A mathematical entropy function was defined to consider the energy terms. The third-order Runge-Kutta method was applied to time marching. The third-order Runge-Kutta method was applied to time marching. As benchmarks, Kovasznay flow, a lid-driven cavity, and the inviscid Taylor–Green vortex were simulated, and good accuracy of the proposed scheme was reported.

Mousa et al. [10] also integrated artificial compressibility with the line approach and utilized the resulting scheme to simulate two-dimensional incompressible flow inside a cavity. The proposed method was used to model magneto-hydrodynamic free convection, and the results were validated against other experimental and numerical results, with good agreement observed. The influence of the artificial compressibility factor was not studied in this work. Jiang et al. [11] combined an artificial compressibility scheme with the Crank-Nicolson Leapfrog method to solve incompressible laminar flows. The results demonstrated that this integrated method is long-term stable and is a second-order scheme. Nagata et al. [12] conducted a theoretical examination of the artificial compressibility system. Following that, a novel explicit artificial compressibility method was proposed, and a two-dimensional flow inside a squared cavity was simulated using the proposed method. Theoretically, they calculated the pseudo speed of sound and used it in their simulations. Reynolds number was considered between 100 and 10,000. The steady and unsteady flows were simulated in this work, and the obtained results demonstrated the accuracy of the proposed method in both cases. Dupuy et al. [13] simulated turbulent flow using an artificial compressibility scheme. To achieve numerical discretization, the finite difference approach in a staggered system was applied. The artificial mach number was defined in this work, and the findings indicated that by extrapolating the artificial mach number, the convergence rate may be accelerated. Additionally, the results demonstrated that the method used to acquire convective fluxes has an effect on the accuracy of the results. Due to the fact that the proposed technique was an explicit mechanism, low memory was required.

Zhang et al. [14] created a new scheme by combining the artificial compressibility method and the spectral collocation approach. They used this scheme to solve three-dimensional incompressible flows numerically. The proposed numerical scheme’s findings were compared to an analytical solution. The combined scheme offered a high degree of accuracy, and the accuracy exponentially increased with time marching. Various computational schemes have been used in the last few decades to discretize virtual governing equations, including averaging methods [15]. According to the history of computing schemes, characteristic methods were initially employed to analyze compressible equations. When the virtual compressibility was added to the incompressible governing equations, the equations were examined using the characteristics approaches as well. Adibi et al. [16,17] studied the Chorin compressibility coupled with a new characteristic base method for two-dimensional flows in order to achieve faster convergence. Their results showed the robustness and high-speed convergence of the proposed scheme. The cavity flow and flow over a cylinder were simulated, and the flow was considered laminar and incompressible. They considered both free and mixed convection in their studies. The artificial compressibility method was considered 5 by the authors.

The literature survey demonstrates the importance of artificial compressibility methods for researchers in different numerical fields. Different researchers also worked on the accuracy of available schemes to improve their accuracy. Marti et al. [18] worked on a finite element-based scheme to improve its accuracy. Three fluid flows with moving boundaries were simulated by the enhanced particle finite element method, and the proposed scheme offered more accuracy. Goona et al. [19] attempted to solve the three-dimensional Poisson equation numerically. By removing error sources, a finite element approach was improved, and a novel scheme was proposed. Numerical simulations of the previous and new proposed schemes were performed, and the results demonstrated that the maximum error can reduce from 8.15% to 0.00091% when the novel scheme is used. Baeza et al. [20] proposed a new scheme for third-order weighted basically nonoscillatory systems. The study revealed that even at a critical time, optimal accuracy occurred for smooth data. As a result, their proposed scheme was based on smoothing the data by adding some nodes with stencil data. Several numerical simulations were conducted using the proposed method, and the results indicated that this scheme is more efficient in error reduction against CPU time. Numerical methods are widely used to investigate flow patterns and heat transfer in a variety of industries due to their advantages, including accessibility to a variety of solutions, low cost, and sensitivity analysis [21–25]. In this context, different numerical methods have been considered to increase the accuracy of the results, reduce the computational costs, and couple the results with analytical methods [26–30]. In these numerical studies, a method has been introduced to solve linear and nonlinear differential equations, which requires a high level of mathematical ability [31–33]. Researchers have combined nanofluids [34–40] and MHD [41–46] with mechanical engineering to improve efficiency in various applications. As a result of these approaches, optimal states for various issues have been found [47–50]. It is even possible for specialists to work and solve common problems in different fields such as medicine [51–53], biotechnology [54,55], and micro-technology [56–58] as a result of advances in numerical calculations and the use of the recent experiences mentioned above. These numerical methods offer acceptable accuracy in solving problems related to thermo-elastic and porous media [59–61]. The accuracy and convergence speed of every numerical scheme is very crucial. Computational fluid dynamics software uses techniques with high convergence speed and good accuracy. In numerical simulations, the artificial compressibility scheme was widely employed. However, the artificial compressibility factor has an impact on its convergence speed and accuracy. Therefore, it is critical to determine the optimal artificial compressible coefficient. Although the artificial compressibility approach is used in various numerical simulations, previous research has not examined the effect of the compressibility factor on the accuracy of results and convergence speed. The effect of the artificial compressibility coefficient on the accuracy of results and the speed of convergence is investigated in the present study. Another novelty of this study is the combination of the artificial compressibility factor with the fifth-order Runge-Kutta method. This novel scheme is formulated in a finite volume framework and convective and viscous fluxes are obtained in cell borders. Numerical simulations are performed for a wide range of Grashof and Reynolds numbers. Consequently, the results of this paper could be used in the numerical solution of natural, mixed, and forced convection in different applications such as heat exchangers. This research allows researchers in the field of numerical methods to ensure the speed of convergence as well as the accuracy of the results by choosing the optimum value for the artificial compressibility factor. It also can be useful for improving the performance of commercial numerical simulation software including ANSYS and COMSOL.

The second section of this work describes the governing equations and the novel numerical technique. In the third section, the problem-solving process and several code components are discussed. In the subsequent stage, the results of numerical simulations are obtained and compared to find suitable answers to the research questions.

2  Problem Formulation

The following equations are applied to two-dimensional and incompressible fluid flow with heat transfer. Density and viscosity are assumed constant for Newtonian fluids with laminar flows. The first equation represents continuity, and the second and third ones signify momentum equations. The fourth equation expresses the energy equation [62].

∂u∂x+∂v∂y=0(1)

∂u∂t+u∂u∂x+v∂v∂y=−1ρ∂p∂x+v(∂2u∂x2+∂2v∂y2)(2)

∂u∂t+u∂u∂x+v∂v∂y=−g−1ρ∂p∂y+v(∂2u∂x2+∂2v∂y2)(3)

∂T∂t+u∂T∂x+v∂T∂y=1ρCp(u∂p∂x+v∂p∂y)+kρCp(∂2T∂x2+∂2T∂y2)+μρCpϕ(4)

where ϕ=[∂u∂x+∂v∂y]2+2[(∂u∂x)2+(∂v∂y)2].

A virtual compressibility sentence (1β∂p∂t) is added to the equation, which tends to zero in the final steady-state solution of this governing equation. In this case, the continuity equation will change as follow:

1β∂p∂t+∂u∂x+∂v∂y=0(5)

For many thermos-flows, a faster convergence can be achieved using the Boussinesq model instead of defining density as a function of temperature.

(ρ−ρo)g≅−ρoβ(T−To)g(6)

where ρo is the (constant) flow density, To is operating temperature, and β is the coefficient of thermal expansion. The above Eq. (4) is obtained using the Boussinesq approximation to remove the buoyancy expression from the momentum equation. This approximation is accurate as long as the change in density is small. To calculate time-dependent thermos flow in confined environments, the Boussinesq approximation is required. As long as the temperature change is small, this model can also be used for steady-state problems. The Boussinesq model cannot be utilized to solve problems involving chemicals or combustion. The governing equation in dimensionless vector form is given by

∂Q∂t+∂F∂x+∂G∂y=∂R∂x+∂S∂y+H(7)

Q¯=[puvT],F=[βup+u2uv(T−EcP)u],G=[βvuyp+u2(T−EcP)v],R=1Re[0∂u∂x∂v∂x1Pr∂T∂x](8)

S=1Re[0∂u∂y∂v∂y1Pr∂T∂y],H=[00GrRe2TEcReϕ],ϕ=[∂u∂x+∂v∂y]2+2[(∂u∂x)2+(∂v∂y)2](9)

In the above relations, Q is the vector of the flow variables, and F, G, R, and S are the convective and viscous flux vectors, respectively; β is the virtual compressibility parameter; Re represents the Reynolds number, and H is the source vector. The set of Eqs. (7)–(9) is obtained using Eq. (10) and the * sign is omitted from all parameters for simplicity. The dimensionless pressure concludes the pressure and the gravitational force of the earth. In other words, the static pressure is obtained first by adding the gravity of the earth to the pressure, and then the dimensionless static pressure is obtained.

u∗=uUref,v∗=vUref,x∗=xUref,y∗=yUref,p∗=p+ρgy−p¥ρrefU2ref,t∗=tUrefLref,T∗=T−Tref1Tref2−Tref1(10)

In Eq. (4), the numbers Gr, Ec, Pr, and Re are dimensionless numbers that are defined by the following relations:

Re=ρrefUrefLrefμ,Pr=Cpμk,Ec=u2refCp(Tref2−Tref1),Gr=βexg(Tref2−Tref1)L3refυ2(11)

where Uref is the upper wall velocity, Lref is the length of the upper wall, and Tref1,Tref2 are the high and low temperatures, respectively. By integrating the governing equations on the cell surfaces, the following equation is obtained:

∂∂t∬QdA+∬(∂(F−R)∂x+∂(G−S)∂y)=∬HdA(12)

The following equation is produced by applying the Green theorem and converting the double integral on the cell surface to the integral on the cell boundaries.

∂∂t∬QdA+∮((F−R)dy−(G−s)dx)=∬HdA(13)

The discretization of the resulting equation in a two-dimensional space with quadrilateral cells leads to the following relation for cells (i, j).

Aij∂Qijk∂t=[∑m=14((Fijk−Rijk)Δy−(Gijk−Sijk)Δx)m−AijHijk],k=1,2,3,4,(14)

The convective terms are obtained by the second-order averaging scheme.

λ∗=(λi,j+λi+1,j)/2→λ=p,u.v−and−T(15)

The fourth-order averaging scheme is used to calculate viscous fluxes on the boundary of secondary cells.

∂η∂x|AB=1S∬s∂η∂xds=1S∮ηdy=1S∑k=14ηΔyk→η=u,v,w(16)

The fifth-order Range-Kutta has been used for time discretization Eq. (12). The main reason behind using this method is that it has a higher convergence rate and larger stability range compared to simple time discretization. Additionally, because the amounts collected in the previous stage and the initial stage are used to determine the quantities in each step, there is no need to store the quantities in the middle stages.

Qijk(p)=Qijk(n)−F(Qijk(p−1)),α=14,16,38,12,1→p=1,2,…,5(17)

In the case of the upper index, P signifies the five steps of the Rang-Kota method, and the index n represents the time step.

A square rectangular cavity with a movable lid is chosen in this study as a benchmark, illustrated in Fig. 1. The fluid contained within the cavity is referred to as an incompressible nanofluid (water-Tio2 nanofluid). Upper and lower walls move at a constant speed of U, while the remaining walls are fixed and insulated. The no-slip law is used to determine the velocity and temperature boundary conditions whereas a second-order extrapolation is used to calculate the pressure.

Figure 1: The cavity and its boundary conditions

3  Materials and Methods

The current work is performed in several steps. First, the dimensionless governing equations are discretized in a finite volume framework. Then viscous fluxes and convective fluxes are obtained by a fourth-order and second-order averaging scheme, respectively. Time marching is done by a fifth-order Runge-Kutta method and the velocity, pressure, and temperature fields are obtained. The grid independence and validation of the results are done. In the next step, simulations are repeated for different artificial compressibility factors. Finally, the accuracy and convergence speed are compared for different artificial compressibility factors to find its optimum value. A code has been written in FORTRAN software to perform the above steps. This code contains different subroutines and the main part. The subroutines are then called in the main part. Gird generation is done in the first subroutine. Initial and boundary conditions are given in the second subroutine. Convective fluxes and viscous fluxes are obtained using the third and fourth subroutines. Charts are obtained using these files’ data via TECPLOT software. The flowchart of this work is shown in Fig. 2.

Figure 2: Different steps of current work

4  Results and Discussion

Simulations were run in four distinct grids in order to ensure the accuracy of the numerical solution. The simulations are carried out for β = 5, Gr = 40, Pr = 6.161, Re = 20, and the Nusselt number is calculated, as shown in Fig. 3.

Figure 3: Grid independence (Nusselt number at β=5, Gr = 40, Pr = 6.161, Re = 20)

To find the optimum grid, high accuracy and less computational time are considered. As seen in Fig. 3, the grids 40 × 40, 60 × 60, and 80 × 80 produce approximately similar results. In a 20 × 20 grid, the dimensions of the cells become larger, the error becomes much higher, and thus the grid independence is not obtained. The grid with 40 × 40 (=1600 cells) is chosen for this study because it provides the optimal balance of accuracy and computing time. Finally, validation has been done with the results of Iwatsu et al. [63], as shown in Fig. 4a. The velocity in the vertical centerline of the cavity is calculated and compared to that of Iwatsu et al. [63]. A good agreement is seen between the present and Iwatsu’s results. Different schemes could be used in numerical simulations. Our focus is on the artificial compressibility method in this work. But we compared three different schemes in the first step. The convergence speed of the artificial compressibility method is compared with that of SIMPLE and Lattice Boltzmann schemes in Fig. 4b. This comparison is done in the optimum compressibility factor. Results show better convergence speed for the artificial compressibility method. The error is obtained by

Error=∑i=1Nox∑j=1Noy(ξi,jk+1−ξi,jk)Nox.Noy(18)

Figure 4: a) Comparison of present and Iwatsu’s [63] results in terms of the velocity on the vertical centerline of the cavity, b) Comparison of convergence speed of three different schemes

One of the well-known problems of viscous flow is the flow inside a square cavity. Depending on the Reynolds number, one or more vortices have been formed inside the cavity. As the Reynolds number increases, the number of vortices increases. The incompressible and viscous steady flow inside the cavity is a well-known example of a case test used to validate a novel numerical scheme. The simple and the clustered quadrilateral grid have been used in this work. The natural and mixed convection nanofluid flows are considered in this paper. First, natural convection (Gr = 0) is considered and then mixed convections with non-zero Grashof numbers are simulated. In the first case, the continuity and momentum equations are solved simultaneously and the velocity and pressure fields are determined. By placing the velocity and pressure fields in the energy equation, the temperature field is determined. Nevertheless, in mixed convection, the momentum equation is a function of Grashof number and temperature. Therefore, all the equations should be solved simultaneously to determine the pressure, velocity and temperature fields. In numerical techniques, the Courant number is taken into account and is defined as follows:

C=uΔtΔx(19)

Courant number is a factor used to determine the convergence of the numerical method. Convergence speed improves as the Courant number is increased. However, the numerical method diverges with higher Courant numbers. Mathematical approaches can be used to determine the range of courant numbers required for convergence in simpler cases with simple governing equations. Nevertheless, in more complex equations, this range is determined by numerical experiments. To find this range, the written code for the numerical scheme is run for different Courant numbers until the code output diverges and the highest achievable value for Courant numbers is determined. Convergence history for various Courant numbers with β=5 is derived and displayed in Fig. 5a. As seen in Fig. 5a, increasing the courant number speeds up convergence. The simulation has been run for different artificial compressibility factors (0.1, 1, 5, and 20) at a constant Courant number and the results are depicted in Fig. 5b. The convergence speed is seen to be high for small artificial compressibility factors. For each artificial compressibility factor, the simulation has been repeated to the maximum Courant number (Fig. 5b). The results indicate that there is no substantial difference in the rate of convergence.

Figure 5: a) Convergence history for various artificial compressibility factor (Gr = 0, Re = 100, Pr = 6.161), b) Convergence history for different artificial compressibility factors and maximum possible Courant number (Gr = 0, Re = 100, Pr = 6.161)

To assess the effect of the artificial compressibility coefficient on the results, the friction coefficients at the bottom plate are calculated and compared for various artificial compressibility coefficients, as illustrated in Fig. 6a. In the figure, the coefficient of friction gradually increases for the bottom plate and the page, changes behavior and the coefficient of friction decreases. As seen in Fig. 6a, the compression ratio has no effect on the physical outcomes. This is because this coefficient appears as a pressure-time gradient coefficient in the continuity relation. Because these problems are steady and have a steady-state solution, the pressure gradient with respect to time is zero, and the artificial compressibility coefficient impact is eliminated during the last step of numerical solutions. Additionally, the equations for unsteady flows are utilized to obtain the numerical solution for steady-state flows. In other words, the pseudo time marching is done on unsteady equations to reach a solution for steady-state flow problems. This approach begins with an arbitrary initial condition for the numerical solution of the equations and determines the solutions in each step based on the preceding step data. The temperature variation along the horizontal centerlines of the cavity is obtained by numerical solution of the energy equation. Obtained temperature is shown in Fig. 6b at different artificial compressible factors. The results show that obtained temperature is the function of the artificial compressibility coefficient. In other words, the results show that the effect of the artificial compressibility coefficient has not been completely eliminated in the final step of numerical simulation. This implies that the results are converged in such a way that the pressure gradient coefficient is not zero in the final steps of the numerical solution.

Figure 6: a) Friction coefficient variation at the top plate with different artificial compressibility (Gr = 0, Re = 100, Pr = 6.16), b) Temperature variation in the horizontal centerline of the cavity with different artificial compressibility factors (Gr = 0, Re = 100, Pr = 6.161)

It is noted that the deviation from the results occurs when the artificial compressibility coefficient is changed abruptly. With 200 times the artificial compressibility coefficient, the resulting change is observed (Fig. 6b). Therefore, researchers who use the artificial compression method for numerical simulations must first eliminate the artificial compressibility factor effect on the results before proceeding with other simulations. In the other words, the simulation must be performed using different artificial compressibility factors, and once the effect of the compressibility factor is eliminated from the results, the subsequent simulations must use this artificial compressibility factor. An artificial compressibility factor of 1 to 5 can be a good start for simulations, according to the results shown in Fig. 6b.

The governing equations are numerically solved, and velocity, pressure and temperature fields are obtained. The resulting temperature field is used to display the isotherms (Fig. 7) using the “Tecplot” software. Fig. 7 demonstrates that the isotherms are perpendicular to the left and right walls. This is because the left and right walls are insulated. In addition, the magnitude of the artificial compressibility factor shows a minimal effect on isotherms. Mixed convection with a high Grashof number is simulated in this step. The momentum equation is a function of the Grashof number and the governing equations are solved simultaneously. The pressure, velocity and temperature fields are obtained by solving the governing equations. The convergence history for different artificial compressibility factors is observed and shown in Fig. 8a. Convergence speed is found slower when the artificial compressibility factor is low (0.1) or high (20). Therefore, the fastest convergence occurs for the artificial compressibility factors of 1 to 5. The results of temperature changes are analyzed in different directions. Fig. 8b shows the temperature variation along the horizontal centerline that passes through the center of the cavity with different β. Temperature changes along the screen are decreased and the majority of changes are associated with beta = 20, with the exception of divergence at beta = 20. To survey the obtained results in Fig. 8b, the simulations are conducted with a grid of 6400 cells and the results are shown in Fig. 9a. As seen in the figure, the temperature variations are not the same, which is not due to the number of cells used in the computational domain. These differences are mainly because of the artificial compressibility term. This term does not tend to zero in some cases, according to the results. The Nusselt number variation is illustrated in Fig. 9b for the top and bottom walls. Except for β=1, the figure reveal that the Nusselt number variation is independent of the artificial compressibility factor.

Figure 7: Isotherms with different artificial compressibility factors (Gr = 0, Re = 100, Pr = 6.161)

Figure 8: a) Convergence history for different artificial compressibility factors (Re = 20, Gr = 40, Pr = 6.161), b) Temperature variation along the horizontal centerline of the cavity with different artificial compressibility factors (Re = 20, Gr = 40, Pr = 6.161)

Figure 9: a) Temperature variation along the horizontal centerline of the cavity with different artificial compressibility factors with 6400 cells (Re = 20, Gr = 40, Pr = 6.161) b) Nusselt number variation on the top wall for different artificial compressibility factor (Re = 20, Gr = 40, Pr = 6.161)

5  Conclusion

In this study, a numerical simulation of a two-dimensional incompressible flow with heat transfer inside a square cavity was investigated using an artificial compressibility method. The effect of the artificial compressibility coefficient on convergence speed and accuracy of results was also assessed. It is found that the accuracy of the results is affected if the artificial compression ratio is either too small or too large. In other words, the artificial compressibility term may not tend to zero at the final steps of time marching. The results also demonstrated that the optimal values of the artificial compressibility factor lie between 1 and 5 for viscous incompressible flows. Using the artificial compressibility method necessitates a measure comparable to network independence, as indicated by the results. In other words, simulations should initially be conducted with varying values of the synthetic compressibility factor, and if the results do not change, simulations should continue. By determining the optimal number for the artificial compressibility factor, researchers in the field of numerical techniques can ensure the speed of convergence as well as the accuracy of the results. This study employed the fifth-order Rung-Kutta method in conjunction with the artificial compressibility method. Simulations could utilize more numerical approaches. The convective fluxes are obtained using the method of averaging. The characteristics-based method could be used.

Funding Statement: The authors extend their appreciation to the Deanship of Scientific Research at King Khalid University for funding this work through the Large Groups Project under grant number RGP. 2/235/43.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.