Direct Collocation with Reproducing Kernel Approximation for Two-Phase Coupling System in a Porous Enclosure

Abstract

The direct strong-form collocation method with reproducing kernel approximation is introduced to efficiently and effectively solve the natural convection problem within a parallelogrammic enclosure. As the convection behavior in the fluid-saturated porous media involves phase coupling, the resulting system is highly nonlinear in nature. To this end, the local approximation is adopted in conjunction with Newton–Raphson method. Nevertheless, to unveil the performance of the method in the nonlinear analysis, only single thermal natural convection is of major concern herein. A unit square is designated as the model problem to investigate the parameters in the system related to the convergence; several benchmark problems are used to verify the accuracy of the approximation, in which the stability of the method is demonstrated by considering various initial conditions, disturbance of discretization, inclination, aspect ratio, and reproducing kernel support size. It is shown that a larger support size can be flexible in approximating highly irregular discretized problems. The derivation of explicit operators with two-phase variables in solving the nonlinear system using the direct collocation is carried out in detail.

1. Introduction

In the literature, there is a bunch of numerical studies on double-diffusive natural convection in enclosures filled with fluid-saturated porous media [1]. It was not until the investigation by Costa that the enclosure of a parallelogram was numerically and theoretically studied in detail [2]. In the practical application, a parallelogrammic enclosure is referred to as the diode of heat and mass transfer owing to the fact that its inclined angle, aspect ratio, and boundary conditions can be designed to meet different purposes concerning various heat and mass transfer characteristics. Nevertheless, the coupling feature of three variables, including temperature, stream function, and concentration, in the problem of double-diffusive natural convection in a parallelogrammic porous enclosure inevitably leads to a highly nonlinear system of partial differential equations (PDEs).

Among the numerical methods for solving PDEs governed by diffusion-convection in parallelogrammic enclosures, the mesh-based methods were adopted early, while the meshfree methods were introduced in recent years. For instance, the finite difference method was employed in the analysis of double-diffusive convection in a porous enclosure with a rectangular shape in 2002 [3]. Around the same time, the finite element method was used to solve the two-dimensional double-diffusive natural convection problem [2]. With the prosperous development of numerical methods in computational mechanics, there was a tendency towards boundary-type methods, in addition to the domain-type methods. To name but a few, the boundary domain integral method was applied to study the double diffusive natural convection in porous media [4]; the boundary element method was introduced to simulate three-dimensional double-diffusive natural convection in porous media [5]. As the meshfree methods get rid of mesh constraint and the strong form methods discard domain integral, a growing number of studies appear and are reported as follows: the local radial basis collocation method was used to solve the two-dimensional double diffusion and natural convection problem in 2012 [6]; the exponentially convergent scalar homotopy algorithm was introduced as the nonlinear solver. Later, in 2018, the generalized finite difference method combined with the Newton–Raphson method was proposed to solve both single thermal natural convection in parallelogrammic enclosures and double-diffusive natural convection in parallelogrammic enclosures filled with fluid-saturated porous media [7].

By examining the work done by Fan et al. [6] and Li et al. [7], the major features in common is the local approximation adopted in the meshfree methods, in which a stable system of equations can be established for nonlinear analysis. As motivated by this idea, the present study proposes the reproducing kernel collocation method in conjunction with the Newton–Raphson method to efficiently and effectively solve the natural convection problems with parallelogrammic enclosures. As there exist several parameters affecting the convergence of the proposed method in solving nonlinear PDEs, only single thermal natural convection in parallelogrammic enclosures filled with fluid-saturated porous media is considered in this study.

In the collocation methods, both global and local functions can be used in the approximation. The commonly adopted global function is the radial basis function (RBF) defined by the Euclidean norm; nevertheless, the shape parameter of RBF cannot be determined uniquely by a given formula so far. As the resulting system is often ill-conditioned with a large condition number, numerical instability might be raised; for more information, see [8,9,10]. In contrast, the reproducing kernel (RK) shape function is a local function; although it maintains an algebraic convergence rate, the resulting system is more stable to yield promising results [11,12,13]. Other types of approximation adopted in the collocation method, such as global expansion in Bernstein polynomials, exponential functions, and Taylor matrix, can be referred to the studies [14,15,16] for solving nonlinear ODE systems. Regarding the two-phase coupling problems, the concept of u-p formulation for describing displacement–pressure behavior was developed early in the Biot theory for the analysis of poromechanics in fluid-saturated media [17]. Later, the u-p formulation was further introduced in the weak-form reproducing kernel particle method (RKPM) to perform the hydro-mechanical analysis in porous media coupled by solid stiffness matrix and compressibility matrix [18]. Unlike the u-p formulation for describing displacement and pressure, the present study first introduces the strong-form reproducing kernel collocation method (RKCM) to the analysis of two-phase coupling system governed by thermal-convection equations.

The structure of this study is the following: the mathematical formulation of the two-phase coupling in a porous enclosure is re-visited in Section 2. The RKCM for the nonlinear system coupled by two variables is introduced with detailed formulation in direct collocation in Section 3. Four numerical examples are given in Section 4. Section 5 concludes the present study.

2. Mathematical Formulation

2.1. Governing Equations

As shown in Figure 1a, the parallelogram is formed by an angle $\theta$ with two sides of length $L$ and $H$. Assuming that the fluid-saturated porous medium has two phases, the components of flux or velocity $Q_x$ and $Q_y$ along the $X$ and $Y$ directions are described by Darcy’s law with the following proportionality to pressure drop:

$$Q_x=-\frac{K}{\mu }\left(\frac{\partial p}{\partial X}\right),$$

(1)

$$Q_y=-\frac{K}{\mu }\left(\frac{\partial p}{\partial Y}+\rho g\right),$$

(2)

where $K$ is the permeability of the isotropic medium; $\mu$, $p$, and $\rho$ are the dynamic viscosity, pressure, and density of the fluid, respectively; and $g$ is the gravitational constant. The stream function $\Psi$ is introduced to relate the flux as given by:

$$\rho Q_x=\frac{\partial \Psi }{\partial Y},$$

(3)

$$\rho Q_y=-\frac{\partial \Psi }{\partial X},$$

(4)

For the sake of convenience, Equations (1) and (2) are normalized as:

$$q_x=\frac{Q_x}{\alpha /H}=\frac{H}{\rho \alpha }\frac{\partial \Psi }{\partial Y}\equiv \frac{\partial \psi }{\partial y},$$

(5)

$$q_y=\frac{Q_y}{\alpha /H}=-\frac{H}{\rho \alpha }\frac{\partial \Psi }{\partial X}\equiv -\frac{\partial \psi }{\partial x},$$

(6)

where the non-dimensional parameters $x=\frac{X}{H}$, $y=\frac{Y}{H}$, and $\psi =\frac{\Psi }{\rho \alpha }$ are used. $\alpha$ is the thermal diffusivity of the porous medium; nevertheless, it is assumed to be the thermal diffusivity of the fluid alone herein, i.e., the effective thermal diffusivity [1,2]. The non-dimensional temperature is given by:

$$T=\frac{T_p-T_{pL}}{T_{pH}-T_{pL}},$$

(7)

where the subscript $p$ in Equation (7) represents the quantity in the physical domain; the subscripts $H$ and $L$ refer to the corresponding higher and lower values. For the non-solute transferring problem, the governing equations are described by the balance laws of linear momentum and thermal energy:

$$\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}+R{a}_T\frac{\partial T}{\partial x}=0,$$

(8)

$$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}-\frac{\partial }{\partial x}\left(q_{x}T\right)-\frac{\partial }{\partial y}\left(q_{y}T\right)=0,$$

(9)

where a steady state is assumed in the above two equations. Introducing Equations (5) and (6) to Equation (9) gives:

$$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}-\frac{\partial \psi }{\partial y}\frac{\partial T}{\partial x}+\frac{\partial \psi }{\partial x}\frac{\partial T}{\partial y}=0,$$

(10)

For natural heat convection in a non-porous medium, Rayleigh number $Ra$ and Darcy number $Da$ are:

$$Ra=\frac{g{\beta }_T\left(T_{pH}-T_{pL}\right)H^3}{\nu \alpha },$$

(11)

$$Da=\frac{K}{H^2},$$

(12)

where ${\beta }_T$ is the thermal expansion coefficient (${\beta }_T\ge 0$) and $\nu$ is the kinematic viscosity of the fluid. Based on Equations (11) and (12), the Darcy-modified Rayleigh number $R{a}_T$ can be defined as:

$$R{a}_T=Da\times Ra=\frac{g{\beta }_T\left(T_{pH}-T_{pL}\right)KH}{\nu \alpha },$$

(13)

which is also the non-dimensional parameter.

2.2. Boundary Conditions

Referring to the general model depicted in Figure 1b, the boundary walls of the porous medium in a parallelogram are assumed to be impermeable, which indicates $\psi =0$. As such, the non-dimensional components of velocity $q_x$ and $q_y$ are zero at boundary walls, while they can be evaluated by using Equations (5) and (6). The explicit boundary conditions for $T$ and $\partial T/\partial n$ together with $\psi =0$ are described as follows:

$$\begin{array}{l}\psi \left(0,y\right)=0,\\ T\left(0,y\right)=1,\end{array}$$

(14)

$$\begin{array}{l}\psi \left(\frac{\cos \theta }{H/L},y\right)=0,\\ T\left(\frac{\cos \theta }{H/L},y\right)=0,\end{array}$$

(15)

Along the inclined bottom side:

$$\begin{array}{l}\psi =0,\\ \frac{\partial T}{\partial n}=0,\end{array}$$

(16)

Along the inclined upper side:

$$\begin{array}{l}\psi =0,\\ \frac{\partial T}{\partial n}=0.\end{array}$$

(17)

In Equations (16) and (17), $n$ denotes the non-dimensional normal to the inclined wall.

3. RKCM for a Nonlinear Coupling System

3.1. Reproducing Kernel Approximation for Two-Phase Variables

To efficiently solve the nonlinear system of equations, the RK shape functions are introduced for two variables $\psi$ and $T$ at $N_s$ source points as follows:

$$\begin{array}{l}\psi \left(x\right)={\displaystyle \sum _{I=1}^{N_s}{\Psi }_I\left(x\right)} a_I^{\psi },\\ T\left(x\right)={\displaystyle \sum _{I=1}^{N_s}{\Psi }_I\left(x\right)} a_I^T,\end{array}$$

(18)

where $a_I^{\psi }$ and $a_I^T$ are the generalized coefficients. ${\Psi }_I\left(x\right)$ is the RK shape function consisting of a monomial basis $H\left(x-x_I\right)$, coefficient vector $b\left(x\right)$, and kernel function ${\phi }_a\left(x-x_I\right)$ of the following form:

$${\Psi }_I\left(x\right)=H^{\mathrm{T}}\left(x-x_I\right)b\left(x\right){\phi }_a\left(x-x_I\right),$$

(19)

where $b\left(x\right)$ is found by the reproducing conditions listed below:

$${\displaystyle \sum _{I=1}^{N_s}{\Psi }_I\left(x\right)x_I^{\alpha }}=x^{\alpha },\left|\alpha \right|\le p,$$

(20)

where $\alpha$ is the multi-index defined as $\left|\alpha \right|={\alpha }_1+{\alpha }_2$, and $p$ is the order of $H\left(x-x_I\right)$. The detailed derivation of RK shape functions is referred to previous work [19]. The explicit form of RK shape functions is:

$${\Psi }_I\left(x\right)=H^{\mathrm{T}}\left(0\right)M^{-1}\left(x\right)H\left(x-x_I\right){\phi }_a\left(x-x_I\right),$$

(21)

with the moment matrix $M\left(x\right)$ given by:

$$M\left(x\right)={\displaystyle \sum _{I=1}^{N_s}H\left(x-x_I\right)H^{\mathrm{T}}\left(x-x_I\right){\phi }_a\left(x-x_I\right)},$$

(22)

To ensure continuity, the high order B-spline kernel function is usually required in strong form collocation methods. The quintic B-spline kernel function used herein is given by:

$${\phi }_a\left(s\right)=\{\begin{array}{cc}\frac{11}{20}-\frac{9s^2}{2}+\frac{81s^4}{4}-\frac{81s^5}{4},\hfill & \hfill 0\le s<\frac{1}{3}\\ \frac{17}{40}+\frac{15s}{8}-\frac{63s^2}{4}+\frac{135s^3}{4}-\frac{243s^4}{8}+\frac{81s^5}{8},\hfill & \hfill \frac{1}{3}\le s<\frac{2}{3}\\ \frac{81}{40}-\frac{81s}{8}+\frac{81s^2}{4}-\frac{81s^3}{4}+\frac{81s^4}{8}-\frac{81s^5}{40},\hfill & \hfill \frac{2}{3}\le s<1\\ 0,\hfill & \hfill s\ge 1\end{array},$$

(23)

where the normalized nodal distance $s=\Vert x-x_I\Vert /a$, and a denotes the radius of a circular support. For the $p$th order monomial basis, the RK support size is chosen to be proportional to the average nodal distance $h$, i.e., $a=\left(p+\delta \right)h$ with $\delta >0$ for any discretization.

The two equations in Equation (18) can be combined as:

$$\left[\begin{array}{c}\psi \left(x\right)\\ T\left(x\right)\end{array}\right]=\left[\begin{array}{cc}{\Psi }^{\mathrm{T}}\left(x\right)& 0\\ 0& {\Psi }^{\mathrm{T}}\left(x\right)\end{array}\right]a,$$

(24)

with:

$$a=\left[\begin{array}{c}{a}^{\psi }\\ a^T\end{array}\right],$$

(25)

$${\Psi }^{\mathrm{T}}\left(x\right)=\left[\begin{array}{cccc}{\Psi }_1\left(x\right)& {\Psi }_2\left(x\right)& \cdots & {\Psi }_{N_s}\left(x\right)\end{array}\right],$$

(26)

$$a^{\psi }=\left[\begin{array}{c}{a}_1^{\psi }\\ a_2^{\psi }\\ ⋮\\ a_{N_s}^{\psi }\end{array}\right],a^T=\left[\begin{array}{c}{a}_1^T\\ a_2^T\\ ⋮\\ a_{N_s}^T\end{array}\right],$$

(27)

Substituting Equation (24) into Equations (8) and (10) leads to:

$$\left[{\Psi }_{,xx}^{\mathrm{T}}\left(x\right)+{\Psi }_{,yy}^{\mathrm{T}}\left(x\right)\right]a^{\psi }+R{a}_T{\Psi }_{,x}^{\mathrm{T}}\left(x\right)a^T=0,$$

(28)

$$\left[{\Psi }_{,xx}^T\left(x\right)+{\Psi }_{,yy}^T\left(x\right)\right]a^T-{\Psi }_{,y}^{\mathrm{T}}\left(x\right)a^{\psi }\cdot {\Psi }_{,x}^{\mathrm{T}}\left(x\right)a^T+{\Psi }_{,x}^{\mathrm{T}}\left(x\right)a^{\psi }\cdot {\Psi }_{,y}^{\mathrm{T}}\left(x\right)a^T=0,$$

(29)

Similarly, substituting Equation (24) into Equations (14)–(17) gives the following equations:

$$\begin{array}{l}{\Psi }^{\mathrm{T}}\left(x\right)a^{\psi }=0,\\ {\Psi }^{\mathrm{T}}\left(x\right)a^T=1,\end{array}$$

(30)

along $x=0$, and:

$$\begin{array}{l}{\Psi }^{\mathrm{T}}\left(x\right)a^{\psi }=0,\\ {\Psi }^{\mathrm{T}}\left(x\right)a^T=0,\end{array}$$

(31)

along $x=\frac{\cos \theta }{H/L}$, and:

$$\begin{array}{l}{\Psi }^{\mathrm{T}}\left(x\right)a^{\psi }=0,\\ {\Psi }_{,n}^{\mathrm{T}}\left(x\right)a^T=0,\end{array}$$

(32)

for the inclined bottom side, and:

$$\begin{array}{l}{\Psi }^{\mathrm{T}}\left(x\right)a^{\psi }=0,\\ {\Psi }_{,n}^{\mathrm{T}}\left(x\right)a^T=0.\end{array}$$

(33)

for the inclined upper side.

3.2. Newton–Raphson Collocation Method for Nonlinear System

From Equations (28) and (29), it is obvious that the system is nonlinear and coupled by two phases. To this end, the Newton–Raphson method is introduced under the collocation framework, which is designated as a Newton–Raphson collocation method in this work.

As a beginning of the formulation, three sets of collocation points are defined: $x_p(p=1\dots \dots N_p)$ for governing equations in the domain, $x_q(q=1\dots \dots N_q)$ for Neumann boundary conditions, and $x_r(r=1\dots \dots N_r)$ for Dirichlet boundary conditions. Note that the same discretization of collocation points $N_c$ and source points $N_s$ is adopted herein; for two field variables, the total number of collocation equations is $2N_c$, which is equal to $N_p+N_q+N_r$. By using the direct collocation, Equations (28)–(33) can be recast as follows:

$$A=\left[\begin{array}{c}{A}_p\left(x_p\right)\\ A_q\left(x_q\right)\\ A_r\left(x_r\right)\end{array}\right]=0,$$

(34)

where:

$$A_p\left(x_p\right)={\left[\begin{array}{c}{A}_p^{\psi }\left(x_p\right)\\ A_p^T\left(x_p\right)\end{array}\right]}_{N_p\times 1},$$

(35)

with:

$$A_p^{\psi }\left(x_p\right)=\left[{\Psi }_{,xx}^{\mathrm{T}}\left(x_p\right)+{\Psi }_{,yy}^{\mathrm{T}}\left(x_p\right)\right]a^{\psi }+R{a}_T{\Psi }_{,x}^{\mathrm{T}}\left(x_p\right)a^T,$$

(36)

$$\begin{array}{l}{A}_p^T\left(x_p\right)=\left[{\Psi }_{,xx}^T\left(x_p\right)+{\Psi }_{,yy}^T\left(x_p\right)\right]a^T-{\Psi }_{,y}^{\mathrm{T}}\left(x_p\right)a^{\psi }\cdot {\Psi }_{,x}^{\mathrm{T}}\left(x_p\right)a^T\\ +{\Psi }_{,x}^{\mathrm{T}}\left(x_p\right)a^{\psi }\cdot {\Psi }_{,y}^{\mathrm{T}}\left(x_p\right)a^T\end{array}$$

(37)

Nevertheless, to obtain the general expression of boundary collocation equations for wider application, the prescribed non-zero values, but not limited to, of both Neumann and Dirichlet boundary conditions are presented. That is, ${\overline{\psi }}_n$, ${\overline{T}}_n$, $\overline{\psi }$, and $\overline{T}$ are assumed in the following derivation. Therefore:

$$A_q\left(x_q\right)=\left[\begin{array}{c}{A}_q^{\psi }\left(x_q\right)\\ A_q^T\left(x_q\right)\end{array}\right]={\left[\begin{array}{c}{\Psi }_n^{\mathrm{T}}\left(x_q\right)a^{\psi }-{\overline{\psi }}_n\\ {\Psi }_n^{\mathrm{T}}\left(x_q\right)a^T-{\overline{T}}_n\end{array}\right]}_{N_q\times 1},$$

(38)

$$A_r\left(x_r\right)=\left[\begin{array}{c}{A}_r^{\psi }\left(x_r\right)\\ A_r^T\left(x_r\right)\end{array}\right]={\left[\begin{array}{c}{\Psi }^{\mathrm{T}}\left(x_r\right)a^{\psi }-\overline{\psi }\\ {\Psi }^{\mathrm{T}}\left(x_r\right)a^T-\overline{T}\end{array}\right]}_{N_r\times 1},$$

(39)

As the collocation system is coupled in matrix $A_p$, the nonlinear version of direct collocation is established by using the Newton–Raphson method given by:

$$a^{k+1}=a^k-{\left(J^k\right)}^{-1}{A}^k,$$

(40)

in which the superscript $k$ denotes the step of iteration and $J$ is the Jacobian matrix defined as $J=\frac{\partial A}{\partial a}$. Thus, the Jacobian matrix is derived as:

$$J=\left[\begin{array}{c}{J}_p\\ J_q\\ J_r\end{array}\right]=\left[\begin{array}{cc}\frac{\partial A_p}{\partial a^{\psi }}& \frac{\partial A_p}{\partial a^T}\\ \frac{\partial A_q}{\partial a^{\psi }}& \frac{\partial A_q}{\partial a^T}\\ \frac{\partial A_r}{\partial a^{\psi }}& \frac{\partial A_r}{\partial a^T}\end{array}\right],$$

(41)

in which:

$$J_p={\left[\begin{array}{cc}\frac{\partial A_p^{\psi }}{\partial a^{\psi }}& \frac{\partial A_p^{\psi }}{\partial a^T}\\ \frac{\partial A_p^T}{\partial a^{\psi }}& \frac{\partial A_p^T}{\partial a^T}\end{array}\right]}_{N_p\times 2N_s}=\left[\begin{array}{cc}{\Psi }_{,xx}^{\mathrm{T}}\left(x_p\right)+{\Psi }_{,yy}^{\mathrm{T}}\left(x_p\right)& R{a}_T{\Psi }_{,x}^{\mathrm{T}}\left(x_p\right)\\ J_{p1}& J_{p2}\end{array}\right],$$

(42)

with:

$$\begin{array}{l}{J}_{p1}=-{\Psi }_{,y}^{\mathrm{T}}\left(x_p\right)T_{,x}\left(x_p\right)+{\Psi }_{,x}^{\mathrm{T}}\left(x_p\right)T_{,y}\left(x_p\right),\\ J_{p2}={\Psi }_{,xx}^{\mathrm{T}}\left(x_p\right)+{\Psi }_{,yy}^{\mathrm{T}}\left(x_p\right)-{\psi }_{,y}\left(x_p\right){\Psi }_{,x}^{\mathrm{T}}\left(x_p\right)+{\psi }_{,x}\left(x_p\right){\Psi }_{,y}^{\mathrm{T}}\left(x_p\right),\end{array}$$

(43)

and:

$$J_q={\left[\begin{array}{cc}\frac{\partial A_q^{\psi }}{\partial a^{\psi }}& \frac{\partial A_q^{\psi }}{\partial a^T}\\ \frac{\partial A_q^T}{\partial a^{\psi }}& \frac{\partial A_q^T}{\partial a^T}\end{array}\right]}_{N_q\times 2N_s}=\left[\begin{array}{cc}{\Psi }_n^{\mathrm{T}}\left(x_q\right)& 0\\ 0& {\Psi }_n^{\mathrm{T}}\left(x_q\right)\end{array}\right],$$

(44)

$$J_r{}_{}={\left[\begin{array}{cc}\frac{\partial A_r^{\psi }}{\partial a^{\psi }}& \frac{\partial A_r^{\psi }}{\partial a^T}\\ \frac{\partial A_r^T}{\partial a^{\psi }}& \frac{\partial A_r^T}{\partial a^T}\end{array}\right]}_{N_r\times 2N_s}=\left[\begin{array}{cc}{\Psi }^{\mathrm{T}}\left(x_r\right)& 0\\ 0& {\Psi }^{\mathrm{T}}\left(x_r\right)\end{array}\right],$$

(45)

However, from Equation (40), it can be observed that the computation of the inverse of the Jacobian matrix at each iteration is computationally expansive. As such, the following two-step Newton–Raphson method is adopted [7]:

$$J^k\Delta a^k=-A^k,$$

(46)

$$a^{k+1}=a^k+\Delta a^k.$$

(47)

where $\Delta a^k$ denotes the increment of vector $a^k$ at $k$th step. The explicit expressions for matrices $A$ and $J$ are detailed in Appendix A. The stopping criterion for iteration is given by:

$$\max \left|\Delta a^k\right|<{10}^{-9}.$$

(48)

which is adopted in the present study.

3.3. Remarks

For illustration purpose, a schematic diagram of discretization shown in Figure 2 is considered, where the numbers of collocation points in the parallelogram are marked by $N_v$, $N_h$, and $N_d$ for vertical, inclined, and domain points, respectively. It is noted that the four corner points are included in the present method without causing any difficulty in the numerical simulation, while they were excluded in the generalized finite difference method as reported in [7].

By taking the boundary conditions described in Equations (14)–(17) or Equations (30)–(33) as an example, the total collocation points are $2N_c=2N_d+4N_v+4N_h$ with $N_p=2N_d$, $N_q=2N_h$, and $N_r=4N_v+2N_h$. This indicates that the Jacobian matrix is a square matrix, and the resulting system is determined. As such, the present study adopts the direct collocation method. For the problem in consideration, the same location and number of source points and collocation points are utilized herein, i.e., N_s = N_c. Unlike the weighted collocation method with least-squares formulation in the previous studies [11,12,13,19], there are typically more collocation points than source points in the collocation system, and boundary conditions are imposed by weights to balance the errors in the domain and on the boundary. With the RK approximation in the direct collocation method, a sparse system of Jacobian matrix is reached, thereby making the Newton–Raphson collocation method efficient.

4. Numerical Examples

To validate the proposed method in solving the nonlinear system coupled by two-field variables, four numerical examples with various parallelogrammic geometries are provided in the following study. In the reproducing kernel collocation method, the RK shape functions are constructed by using the quadratic basis with a quintic B-spline kernel function, and the support size of RK shape functions is chosen as $a=3h$ without loss of generality. In particular, the parametric study is conducted in each example to thoroughly understand the effect of nonlinearity and the convergence of approximation.

4.1. Square Domain

The schematic diagrams of the uniform and non-uniform discretizations in a unit square are depicted in Figure 3a,b, respectively. In the uniform discretization, three refinements with $N_s=N_c={20}^2,{30}^2,{40}^2$ are considered in the approximation for $R{a}_T=100$ with initial conditions $\left({\psi }^0,T^0\right)=\left(0,0.5\right)$. Figure 4 compares the number of iterations obtained by three refinements; Figure 5 shows the corresponding contour plots of $\psi$ and $T$. Clearly, the present method converges fast and stably in 7 iterative steps for all discretizations; the contour plots of both $\psi$ and $T$ also exhibit the same pattern. The convergence of approximation is validated through the refinement of discretization. In the following study, $N_c={30}^2$ is used. In Figure 3b, the non-uniformly discretized domain is obtained by adding 3% (disturbing factor $s=3$) disturbance to Figure 3a. By using the same values of $R{a}_T$ and $\left({\psi }^0,T^0\right)$, the number of iterations for various non-uniform discretizations obtained by $s=1,2,3$ are given in Figure 6. For $s=3$, the contour plots of $\psi$ and $T$ are shown in Figure 7, where a similar distribution of $\psi$ and $T$ to those obtained by uniform discretization is observed.

Based on the quadratic basis adopted in constructing the RK shape functions, the RK support size generally ranges from $2.5 \mathrm{h}$ to $4 \mathrm{h}$. To understand the influence of support size for the nonlinear system, the number of iterations corresponding to various support size is summarized in Figure 8 for $R{a}_T=100$ with $\left({\psi }^0,T^0\right)=\left(0,0.5\right)$. It is shown that the number of iterations is 7 for all cases regardless of the support size. As for the Newton–Raphson method, the choice of initial conditions may affect the convergence of approximation. To this end, several sets of initial conditions $\left({\psi }^0,T^0\right)$ are tested, and the corresponding iteration is shown in Figure 9; all approximate solutions converge within 7 iterative steps, although their initial errors in the first step may be quite different. For wider application, various initial values, such as $\left({\psi }_{left}^0,{\psi }_{right}^0,T_{left}^0,T_{right}^0\right)$, are considered as depicted in Figure 10a; the corresponding iteration is shown in Figure 10b. It is found that the initial conditions might affect the path of convergence, and more iterations might be needed to meet the stopping criterion. Nevertheless, most initial conditions yield convergent approximation.

For the nonlinear system in consideration, the parameter $R{a}_T$ governs the nonlinearity, and various values of $R{a}_T$ such as 1 and 60 are considered herein. The comparison of the number of iterations obtained by various $R{a}_T$ is made in Figure 11; it was found that more iteration steps are needed for a larger value of $R{a}_T$. For illustration purpose, the contour plots of $\psi$ and $T$ obtained by $R{a}_T=1$ and $R{a}_T=60$ are presented in Figure 12. From Figure 12a, the contour of $\psi$ is almost symmetric in both vertical and horizontal directions for $R{a}_T=1$; as $R{a}_T$ increases to 60 and 100, the contours of $\psi$ gradually showed anti-symmetric distribution along the line inclined at ${45}^{\circ }$ to the horizontal direction as referred to Figure 5b and Figure 12b. As for the contour of $T$, with increasing $R{a}_T$, its distribution moved from a vertical direction toward a curved orientation as shown in Figure 5e and Figure 12c,d.

4.2. Diamond Domain

A unit diamond with $\theta ={30}^{\circ }$ is discretized uniformly and non-uniformly as shown in the schematic diagrams Figure 13a,b, respectively. The analysis is performed by following the previous example as described in Section 4.1. For $R{a}_T=100$ and $\left({\psi }^0,T^0\right)=\left(0,0.5\right)$, the number of iterations obtained by uniform discretizations with $N_c={20}^2,{30}^2,{40}^2$ are shown in Figure 14; the corresponding contour plots of $\psi$ and $T$ are given in Figure 15. It is shown that only 7 steps are required in the iteration process. In addition, the convergence of approximation is assured by the presence of consistent patterns in the contours, as also shown in the previous work [7]. The non-uniform discretization with $N_c={30}^2$ in Figure 13b is obtained by adding 3% disturbance to the uniform discretization. For various percentage of disturbance $s=1,2,3$, the number of iterations are given in Figure 16. For large disturbance $s=3$, more iterations are required due to a longer path of convergence; the corresponding contour plots of $\psi$ and $T$ are shown in Figure 17, where similar patterns of $\psi$ and $T$ to those obtained by uniform discretization are retrieved, while the smoothness of the contours is reduced.

The influence of initial conditions are investigated by considering $\left({\psi }^0,T^0\right)$ for entire domain and $\left({\psi }_{left}^0,{\psi }_{right}^0,T_{left}^0,T_{right}^0\right)$ for two half domains, and the corresponding results are shown in Figure 18 and Figure 19, with Figure 19a illustrating the arrangement of $\left({\psi }_{left}^0,{\psi }_{right}^0,T_{left}^0,T_{right}^0\right)$. From Figure 18 and Figure 19b, it is observed that the initial conditions may affect the path of convergence, i.e., the number of iterations. As for the influence of $R{a}_T$, Figure 20 compares the number of iterations obtained by various $R{a}_T$ ranging from 1 to 100; the larger value $R{a}_T$ is, the more iterations it requires. The selected contour plots of $\psi$ and $T$ obtained by $R{a}_T=1$ and $R{a}_T=60$ are shown in Figure 21. As $R{a}_T$ increases, the contour of $\psi$ changes from symmetry to anti-symmetry with respect to the line inclined at ${45}^{\circ }$ in the horizontal direction, and its magnitude increases as well, as shown in Figure 15b and Figure 21a,b. The contour of $T$ gradually changes the orientation from left to right with increasing nonlinearity, as can be observed in Figure 15e and Figure 21c,d.

4.3. Parallelogram Domains I and II

As shown in Figure 22a,b, two parallelograms ($L/H=4$) with inclined angles ${30}^{\circ }$ and $-{30}^{\circ }$ are considered, in which the former is denoted as domain I, while the latter is denoted as domain II. For clarity, when comparing the results of two parallelograms in domain I and domain II, the corresponding results are depicted in terms of blue and red lines as will be shown in the following figures. For $R{a}_T=100$ and $\left({\psi }^0,T^0\right)=\left(0,0.5\right)$, three refinements of uniform discretization with $N_c=10\times 40,20\times 50,30\times 60$ are considered. The path of convergence expressed in terms of iteration is depicted in Figure 23. In each domain, the approximation converges with the same number of iterations regardless of whether different points are used in the discretization; moreover, domain I shows slower convergence than domain II. The contour plots of $\psi$ and $T$ obtained by $N_c=20\times 50$ are shown in Figure 24, which are consistent with the results given in [7]. The distributions of both $\psi$ and $T$ in two domains are dissimilar due to different inclination; in particular, the contour of $\psi$ in domain II with a negative inclined angle has two small vortices at the two ends.

The influence of initial conditions on the solutions is investigated by concerning $\left({\psi }^0,T^0\right)$ and $\left({\psi }_{left}^0,{\psi }_{right}^0,T_{left}^0,T_{right}^0\right)$ for $R{a}_T=100$ and $N_c=20\times 50$; the corresponding paths of iteration are shown in Figure 25a,b, respectively. In general, domain I needs more iterations to reach the desired accuracy in approximation as compared with domain II. Nevertheless, domain II might require more iterations under certain initial conditions. The influence of $R{a}_T$ is presented in Figure 26, where more iterations are needed for larger $R{a}_T$; in addition, domain I shows slower convergence than domain II. To further unveil the nonlinearity of the problem with different inclined angles, the contour plots obtained by $R{a}_T=1$, $R{a}_T=60$, and $R{a}_T=100$ are shown in Figure 24, Figure 27 and Figure 28. For both domain I and domain II, the magnitude of $\psi$ increases with larger $R{a}_T$, and the order of contour (nonlinearity) of $T$ increases with larger $R{a}_T$ as well. Specifically, for domain II, the two vortices in the contour of $\psi$ become more pronounced as $R{a}_T$ increases.

The disturbance of discretization is added to test the stability of the method. By using $R{a}_T=100$ and $N_c=20\times 50$ with various $s=1,2,3$, the convergence paths of two domains are compared in Figure 29; again, domain I shows slower convergence than domain II, and a larger $s$ requires more iterations. For illustration purpose, the non-uniform discretization obtained by $s=3$ is presented in Figure 30, and the corresponding contour plots of $\psi$ and $T$ are given in Figure 31. Similar trends in the distribution of both contours $\psi$ and $T$ are observed while the smoothness of the contours is reduced to some extent in comparison with those obtained by uniform discretization. Lastly, a higher level of disturbance $s=4$ is further investigated, as shown in Figure 32, which has not been discussed in the literature. It is found that the RK support size $a=3h$ is not able to yield satisfactory results. Nevertheless, the flexibility of RKCM enables a larger support size in the approximation. When using $a=4h$, the contour plots of $\psi$ and $T$ are shown in Figure 33 with corresponding iteration given in Figure 34. The patterns of the contours are ensured.

5. Conclusions

The Newton–Raphson collocation method with reproducing kernel approximation is proposed to solve the natural convection problem within a parallelogrammic enclosure. Particularly, a nonlinear system is established on the basis of local approximation. In addition to the benchmark problems given in the literature, a unit square enclosure is designed to unveil the parameters in RKCM when nonlinearity is of importance in such a porous enclosure. From the parametric study, the convergence of the method is stable and efficient, regardless of various initial conditions, disturbance of discretization, inclination, aspect ratio, and the RK support size. Additionally, it is observed that the parameter $R{a}_T$ controls the nonlinearity of the system. In all examples, as the value of $R{a}_T$ increases, the magnitude of $\psi$ increases, and so does the nonlinearity of contour $T$. Even so, the present method yields the desired accuracy for a large value of $R{a}_T$. Although the RK support size is not sensitive to uniformly discretized domains, for highly irregular discretization, it is shown that a larger support size in RKCM offers flexibility in yielding accurate approximation. The feasibility of the method for analyzing two-phase coupling problems is, therefore, demonstrated. To keep the article an adequate length, the extension to three-phase coupling problems, i.e., double-diffusive natural convection problems, is left to be investigated by future works.