Numerical Solution of Natural Convection Problems Using Radial Point Interpolation Meshless (RPIM) Method Combined with Artificial-Compressibility Model

Abstract

A numerical method is used to solve the thermal analysis of natural convection in enclosures. This paper proposes the use of an implicit artificial-compressibility model in conjunction with the Radial Point Interpolation Meshless (RPIM) method to mimic laminar natural convective heat transport. The technique couples the pressure with the velocity components using an artificial compressibility model. The RPIM is used to discretize the spatial terms of the governing equation. We solve the semi-algebraic system implicitly in backward Euler pseudo-time. The proposed method solves two test problems—natural convection in the annulus of concentric circular cylinders and trapezoidal cavity. Additionally, the results are validated using experimental and numerical data available in the literature. Excellent agreement was seen between the numerical results acquired with the suggested method and those obtained through the standard techniques found in the literature.

1. Introduction

Nowadays, the development of meshless methods for numerical modeling, especially heat transfer by natural convection, has attracted much attention from researchers. This method is flexible because it does not require a mesh so that any domain shape can be modeled. Coordinate points spread across the space domain are needed as a substitute for a mesh, so this method is much simpler than methods that require a mesh, such as the Finite Element, Finite Volume, Finite Difference, and Lattice Boltzmann methods [1,2]. Sadat and Couturier [3] developed a meshless method called the Diffuse Approximation Method (DAM) for modeling natural convection. Sophy et al. [4] successfully developed this method for modeling three-dimensional convection. At almost the same time, Liu and Wu [5] developed a meshless method called the Point Interpolation Method (PIM); this method uses low-order polynomial approximations such as the DAM method. The results of this method were proven to be suitable only for calculations with low Rayleigh Numbers. Wijayanta and Pranowo [6] used a localized meshless with Radial Basis Function (RBF) as the interpolation function; the calculation results showed this method’s high accuracy.

The method above solves the model equation in “strong form”; this approach is entirely accurate but sensitive to changes in numerical parameters, so it is less robust. Several researchers proposed solving the model equation in “weak form” by integrating it according to the Galerkin procedure. This method was adopted from the Finite Element method (FEM). The weak-form meshless method was initially successfully developed by Belytschko et al. [7,8] for solid mechanics modeling, especially for calculating crack propagation. This method is called element-free Galerkin (EFG), which received much attention from researchers and was developed to solve various problems in the engineering field. Manzin and Bottauscio [9] used the EFG method to calculate the scattering of electromagnetic waves in a non-homogeneous medium. Singh et al. [10] expanded the use of the EFG method for calculating nonlinear conduction heat transfer. The use of the EFG method for numerical calculations of fluid mechanics was carried out by Singh and Jain [9]; they developed a parallel algorithm for the numerical solution of the Steady Navier–Stokes Equations without involving convection terms. Vlastelica et al. [11] explored the use of the EFG method for numerical calculations of incompressible viscous fluid flows. They observed the effect of the number of integration points on the accuracy of the EFG method.

The development of the EFG method continues to be carried out significantly to increase its accuracy; Zhang et al. [12,13,14] developed the multiscale EFG method. This method divides the numerical solution into two parts, namely the global level solution and the local level solution. Local-level solutions help increase accuracy, especially in parts of the domain that contain sharp gradients. This method has been successfully used to calculate numerical solutions to the Stokes equation [12], incompressible fluid flow [11], and heat transfer by natural convection in complex-shaped domains [12].

Recently, some fluid flow and heat convection issues have been taken into account by the meshless local Petrov–Galerkin (MLPG) approach. The MLPG method was successfully developed to simulate incompressible fluid flow using the stream function-vorticity formulation [15,16,17]. However, this formulation is difficult to apply to complex geometric shapes because the boundary conditions for the flow function and vorticity at the wall are difficult to determine. Apart from that, this method also cannot solve three-dimensional problems. Efforts to solve incompressible fluid flow using primitive variables were pioneered by Atluri and colleagues [18,19]. Najafi et al. [20,21] succeeded in further develo** the MLPG method based on primitive variables for convection heat transfer problems.

The main weakness of the EFG method is that the Dirichlet boundary conditions are difficult to apply directly because the interpolation function uses a least square error approach with low-order polynomial interpolation. Applying Dirichlet boundary conditions requires using the Lagrange Multiplier approach, so the results are not exact. Based on the description above, this research proposes using a meshless numerical method with a “weak form” approach and basis functions using RBF for modeling natural convection in porous media and air. Using RBF basis functions makes it easier to implement Dirichlet boundary conditions. This method is called the Radial Point Interpolation Meshless (RPIM) method [22].

2. The Governing Equations

The model equations used are the two-dimensional Navier–Stokes and energy conservation equations. If the temperature differences that occur are small, then the Boussinesq approach is applicable for modeling the buoyancy force. Thus, the two-dimensional Navier–Stokes equations and energy equations can be written as follows [23]:

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

(1)

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{\partial p}{\partial x}+\frac{Pr}{R{a}^{0.5}}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$$

(2)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{\partial p}{\partial y}+\frac{Pr}{R{a}^{0.5}}\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+Pr\theta$$

(3)

$$\frac{\partial \theta }{\partial t}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}=\frac{1}{R{a}^{0.5}}\left(\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}\right)$$

(4)

where the two-dimensional spatial coordinates are (x,y), the time is t, the pressure is p, the temperature is θ, and the velocity components are (u,v). The dimensionless equations are obtained by dividing dimension variables with the reference variables. The reference variable for length is L_r = H, for velocity V_r = (α/H) Ra^0.5, for the time t_r = (H²/α) Ra^−0.5, for the temperature (θ) is defined as follows: (T − T_c)/(T_h − T)_c.

3. The Radial Point Interpolation Meshless (RPIM) Method

This section briefly introduces the RPIM approximation for function $u\left(x\right)$ utilizing just nodal function value (u_i) at a collection of local nodes. Let $x_i$ be a set of randomly distributed points in a domain (i = 1, 2,..., n), where n is the number of nodes in the support domain. The expression for the function $u\left(x\right)$ approximation can be written as follows [22,24]:

$$u\left(x\right)={\displaystyle \sum _{i=1}^n{\varphi }_i{u}_i}$$

(5)

The above equation becomes:

$$u\left(x\right)=\phi \left(x\right)U_s$$

(6)

$$\phi \left(x\right)=R_Q^T\left(x\right)R_Q^{-1}\left(x_i\right)$$

(7)

when expressed in matrix form.

$$x^T=\left[x,y\right]$$

(8)

$$U_s={\left[u_1,u_2,u_3,\dots ,u_n\right]}^T$$

(9)

$$R_Q\left(x_i\right)=\left[\begin{array}{cccc}{R}_1\left(x_1\right)& R_2\left(x_1\right)& \dots & R_n\left(x_1\right)\\ R_1\left(x_2\right)& R_2\left(x_2\right)& \dots & R_n\left(x_2\right)\\ \dots & \dots & \dots & \dots \\ R_1\left(x_n\right)& R_2\left(x_n\right)& \dots & R_n\left(x_n\right)\end{array}\right]$$

(10)

Whereas R is the matrix holding the radial basis function (RBF) interpolation, U_s is the vector containing the value u at the node point (Figure 1). The vector $x$ represents spatial coordinates. The function chosen for interpolation is the multiquadric RBF function:

$$R_i\left(x,y\right)={\left(r_i^2+C^2\right)}^{0.5}={\left[{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+C^2\right]}^{0.5}$$

(11)

The constants C are the shape parameters.

In Equation (6), $\phi \left(x\right)$ is RPIM shape functions corresponding to the nodal value and given by

$$\phi \left(x\right)={\left[{\varphi }_1\left(x\right),{\varphi }_2\left(x\right),{\varphi }_3\left(x\right),\dots ,{\varphi }_n\left(x\right)\right]}^T$$

(12)

in which

$${\varphi }_k\left(x\right)={\displaystyle \sum _{i=1}^n{R}_i\left(x\right)G_{\left(i,k\right)}}$$

(13)

where $G_{\left(i,k\right)}$ is the element of matrix $R_Q^{-1}$.

After the shape functions are determined, the shape function derivatives can be calculated as follows:

$$\begin{array}{l}\frac{\partial {\varphi }_k}{\partial x}\left(x\right)={\displaystyle \sum _{i=1}^n\frac{\partial R_i}{\partial x}\left(x\right)G_{\left(i,k\right)}}\\ \frac{\partial {\varphi }_k}{\partial y}\left(x\right)={\displaystyle \sum _{i=1}^n\frac{\partial R_i}{\partial y}\left(x\right)G_{\left(i,k\right)}}\end{array}$$

(14)

Similar expressions can be used for the primitive variables in the model equations (u,v,p,θ):

$$\begin{array}{l}u\left(x\right)={\displaystyle \sum _{i=1}^n{\varphi }_i{u}_i};v\left(x\right)={\displaystyle \sum _{i=1}^n{\varphi }_i{v}_i}\\ p\left(x\right)={\displaystyle \sum _{i=1}^n{\varphi }_i{p}_i};\theta \left(x\right)={\displaystyle \sum _{i=1}^n{\varphi }_i{\theta }_i}\end{array}$$

(15)

The RPIM method still requires a mesh background for spatial integration with the Gauss quadrature method, so this method is not a pure meshless method. The shape of the elements in the background mesh can be arbitrary as long as Gauss quadrature integration can be carried out. This study chose the background mesh as a quadrilateral quadrature domain because it is simple (Figure 1).

After the compact support and quadrature domains are determined, the next step is to carry out Galerkin integration using a weight function equal to the shape function. Equation (4) is used as an example of integration:

$${\displaystyle \underset{{\Omega }_Q}{\int }{\varphi }_i\left(\frac{\partial \theta }{\partial t}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}-\frac{1}{R{a}^{0.5}}\left(\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}\right)\right)}d\Omega =0$$

(16)

Substitute the temperature expansion (Equation (14)) into Equation (15):

$$M_{ij}\frac{\partial {\theta }_j}{\partial t}+C_{ij}{\theta }_j+\frac{1}{R{a}^{0.5}}{K}_{ij}{\theta }_j=0$$

(17)

$$M_{ij}={\displaystyle \underset{{\Omega }_Q}{\int }{\varphi }_i{\varphi }_j}d\Omega$$

(18)

$$C_{ij}={\displaystyle \underset{{\Omega }_Q}{\int }{\varphi }_i}\left(u_j\frac{\partial {\varphi }_j}{\partial x}+v_j\frac{\partial {\varphi }_j}{\partial y}\right)d\Omega$$

(19)

$$K_{ij}={\displaystyle \underset{{\Omega }_Q}{\int }\left(\frac{\partial {\varphi }_i}{\partial x}\frac{\partial {\varphi }_j}{\partial x}+\frac{\partial {\varphi }_i}{\partial y}\frac{\partial {\varphi }_j}{\partial y}\right)}d\Omega$$

(20)

Equation (17) is a global semi-algebraic equation. All spatial integration calculation matrices are carried out using the Galerkin procedure in each support domain and then assembled into a global semi-algebraic equation (Equation (17)).

The discretization of Equations (2) and (3) is carried out in the same way; the results are as follows:

$$M_{ij}\frac{\partial u_j}{\partial t}+C_{ij}{u}_j+\frac{Pr}{R{a}^{0.5}}{K}_{ij}{u}_j-{D^x}_{ij}{p}_i=0$$

$$M_{ij}\frac{\partial v_j}{\partial t}+C_{ij}{v}_j+\frac{Pr}{R{a}^{0.5}}{K}_{ij}{v}_j-{D^y}_{ij}{p}_i-Pr{M}_{ij}{\theta }_j=0$$

The main difficulty in solving the incompressible flow equation is due to the pressure being decoupled from the velocity field. Chorin proposed the artificial compressibility method to alleviate this problem [25]. Due to their ease of use and effectiveness, the artificial compressibility for incompressible flows has been widely used by numerous researchers [26,27,28,29]. The pressure field is linked to Equation (1) in the following way using the artificial compressibility method:

$$\frac{\partial p}{\partial t}+\beta \left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0$$

(21)

The derivative of pressure over time in the equation above is a pseudo-variable because the equation above actually only applies to compressible flow. Thus, the transient solution of the system of equations has no physical meaning; only the steady state has physical meaning. Meanwhile, β is an artificial compressibility factor whose value is constant. The β constant is equivalent to the square of the sound speed of the medium, so in this paper, the β is 10⁶.

Discretization of the RPIM method for Equation (19) produces the discrete equation:

$$M_{ij}\frac{\partial p_j}{\partial t}+\beta \left({D^x}_{ij}{u}_j+{D^y}_{ij}{v}_j\right)=0$$

(22)

Time integration of the semi-algebraic equations is carried out using the implicit Euler method so that the calculations are stable:

$$\left[\begin{array}{cccc}{M}_{ij}& \Delta t\beta {D^x}_{ij}& \Delta t\beta {D^y}_{ij}& 0\\ -\Delta t{D^x}_{ij}& A_{ij}+\frac{Pr}{R{a}^{0.5}}\Delta t{K}_{ij}& 0& -\Delta t{M}_{ij}Pr\\ -\Delta t{D^x}_{ij}& 0& A_{ij}+\frac{Pr}{R{a}^{0.5}}\Delta t{K}_{ij}& 0\\ 0& 0& 0& A_{ij}+\frac{1}{R{a}^{0.5}}\Delta t{K}_{ij}\end{array}\right]\left[\begin{array}{c}{p}_j^{n+1}\\ u_j^{n+1}\\ v_j^{n+1}\\ {\theta }_j^{n+1}\end{array}\right]=\left[\begin{array}{c}{M}_{ij}{p}_j^n\\ M_{ij}{u}_j^n\\ M_{ij}{v}_j^n\\ M_{ij}{\theta }_j^n\end{array}\right]$$

(23)

$$A_{ij}=M_{ij}+\Delta t{C}_{ij}$$

(24)

The RPIM algorithm for solving natural convection problems is described as follows [30]:

• Specify maximum pseudo-time for the calculation and time step, the initial condition, boundary condition, primitive variables (u, v, p, θ), node coordinates, basis function parameters, background mesh, Support domain and quadrature domain.
• Calculate the derivative operator matrices in Equations (17)–(20).
• Pseudo-time marching starts.
• Calculate the global matrix and the right side vector in Equation (2).
• Solve the global matrix to obtain the primitive variables ${\left(u,v,p,\theta \right)}^{n+1}$
• Update the pseudo time (n + 1).
• Check if it has reached the maximum time limit, if not back to step 2. If it is go to step 8.
• Finish.

4. Results and Discussion

In the framework of the artificial-compressibility approach, the current improved primitive variable-based RPIM method is presently used to evaluate steady-state two-dimensional laminar natural convection issues. The natural convection heat transfer in the annulus of concentric circular cylinders and the trapezoidal enclosures are used to demonstrate the capability and accuracy of the proposed method.

4.1. Natural Convection in the Annulus of Concentric Circular Cylinders

In this section, we simulate the natural convection heat transfer in a horizontal annulus of concentric circular cylinders at different temperatures using the RPIM model that was previously described. We selected this problem because, besides having curved boundaries and a wealth of experimental and numerical data available in the literature, it is a classical heat transfer problem that can be used to verify the accuracy of the current RPIM method.

In order to match the experimental conditions given in [24], we set the radius of the outer cylinder r_o = 1.625 and the inner one r_i = 0.625. The number of nodes, which are clustered near the walls to capture the boundary layer, is 161 × 25 for all simulations. The Prandtl number is 0.717, time step Δt = 1, and the fixed artificial compressibility factor equals 1 × 10⁶, while the Rayleigh (Ra) numbers varied from 1 × 10³ to 2.0 × 10⁵. The temperature and velocity field are set to zero as the initial conditions. The inner cylinder is heated θ = 1, and the outer cylinder is cooled θ = 0 for the temperature boundary conditions. The Figure 2 shows the domain geometry, boundary conditions and the node distribution.

The temperature distribution obtained numerically is compared with the experimental temperature distributions from the interferograms under the same conditions (Figure 3). The left of Figure 3 is the isotherm plot obtained numerically. In contrast, the right of Figure 3 is an interferogram obtained experimentally. Within the cylinders, the fringes stand in for isotherms. Comparison of isotherm plots shows good agreement between numerical calculations and experimental results.

Figure 4, Figure 5 and Figure 6 show the isothermal, and velocity fields, and stream function at steady state conditions for various Ra = 2.38 × 10³, 3.2 × 10⁴, and 1.02 × 10⁵. The buoyant force causes the fluid at the top of the inner cylinder to flow upwards and forms 2 symmetrical circulation eddies (Figure 4c, Figure 5c and Figure 6c). The Figure 4b, Figure 5b and Figure 6b show the flow direction and also indicate the thickness of the boundary layer. For low Ra number values, the speed of fluid movement is slow due to the influence of buoyancy, and the boundary layer also appears thick. The contribution of the convection term in the governing equations is weak, while the influence of the diffusion term is strong. As a result, the temperature distribution also resembles the distribution of conduction heat transfer. However, if Ra increases, the influence of the convection term becomes stronger, the flow of fluid movement becomes faster, and the thickness of the boundary layer decreases. The (Figure 4b, Figure 5b and Figure 6b). This fast fluid movement also carries heat, so the temperature distribution is more even.

The average Nusselt number (Nu) for various Rayleigh numbers derived using the current RPIM method is compared with those experiments of Kuehn and Goldstein [32], Finite Difference Lattice Boltzmann method of Shi et al. [33], Discontinuous Galerkin method of Pranowo and Deendarlianto [34], and Galerkin Radial Basis Function of Ho-Minh et al. [35], Figure 7. The current results agree with the aforementioned references, as this figure illustrates.

Based on

Table 1. The node convergence of the natural convection in the annulus for Ra = 2.38 × 10³.

Nodes	Nu
	RPIM (Present Study)	Kuehn and Goldstein [32]
101 × 15	1.494	1.38
121 × 20	1.478	1.38
161 × 25	1.402	1.38

, it can be seen that increasing the number of nodes produces a Nusselt number that is closer to the experimental results.

4.2. Natural Convection in a Trapezoidal Enclosure

The second test case, which involves natural convection in a trapezoidal enclosure for various Rayleigh numbers, is utilized to assess the proposed method. The short vertical wall’s height (H) is the characteristic length, and the cavity length (L) is four times H. Meanwhile, the upper cavity’s fixed inclination is 15 degrees. The boundary conditions are as follows: the left wall is heating, and the right wall is cooling [29]. The bottom and top walls are adiabatic walls. Meanwhile, the boundary condition for the velocity component at the wall is zero. For the values Ra = 10³, 10⁴, and Ra = 10⁵, the shape parameter value is C = 20, and the number of nodes is 141 × 41. Meanwhile, for Ra = 10⁶, the C value used is C = 25, and the number of nodes is 151 × 51. Figure 8 illustrates the geometry and the node distribution.

The calculation results of natural convection in the trapezoidal enclosure are presented in the form of isothermal and stream functions. The airflow forms circulation, flowing up from the right due to the buoyant effect due to heating, then flowing to the left and down due to cooling, and finally flowing to the right. An increase in the Ra value causes the dominance of the convective term to strengthen and the speed to increase, so the isothermal forms a pattern that tends increasingly to the left at the top. Conversely, the isothermal tends increasingly to the right at the bottom (Figure 9a, Figure 10a, Figure 11a and Figure 12a).

An increase in the Ra value causes the air circulation pattern to tend to slope because the speed increases (Figure 9b, Figure 10b and Figure 11b); apart from that, it also causes the calculation time to become longer to reach convergent conditions. This can be seen in Figure 13a–d. At the beginning of the calculation, the Nusselt number has a high value due to the effects of sudden heating and cooling, which then decreases and increases slightly. After reaching a steady state, the Nu value remains unchanged again. For the case of Ra = 10⁶, the calculation until the pseudo time t = 500 is not enough. Steady conditions have not been reached; the flow pattern continues to change. This can be seen in Figure 12, where the flow circulation pattern should be gentle, but flow oscillations still occur.

The numerical calculation results of the RPIM method were compared with those of the numerical calculations of the Finite Volume method carried out by Moukalled and Darwish [36] to test the accuracy of the RPIM method. The results are presented in

Table 2. The comparison result of natural convection in trapezoidal enclosure other method results.

Ra	Nu
	Moukalled and Darwish [36]	RPIM (Present Study)
103	0.6153	0.7205
104	1.9920	1.9993
105	4.4310	4.4419
106	8.8400	8.2320

; the comparison shows good agreement.

5. Conclusions

This paper uses the artificial compressibility model and the RPIM method to solve the two-dimensional natural convection heat transfer problems. To assess the method’s accuracy, two benchmark problems involving natural convection in the annulus of concentric circular cylinders and trapezoidal cavities are used as sample calculations. The suggested method’s straightforward methodology makes solving natural convection in an irregular domain easy. There is also reasonable agreement when comparing the isotherms of the temperature fields of the RPIM in the annulus with the experimental interferogram. The proposed method’s outstanding accuracy is demonstrated by comparison with solutions from other numerical methods found in the literature, including the meshless Galerkin-RBF, Finite Difference Lattice Boltzmann (FDLBM), finite volume (FVM), and discontinuous Galerkin finite element (DGFEM) methods.

We plan to apply this technology in future studies to fuel cell modeling [37], two-phase thermofluid flow [38,39], and nanofluid simulation [40]. We also intend to employ parallel processing to accelerate the numerical calculating process.