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

: The direct strong-form collocation method with reproducing kernel approximation is introduced to efﬁciently and effectively solve the natural convection problem within a parallelogrammic enclosure. As the convection behavior in the ﬂuid-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 ﬂexible 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.


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

Governing Equations
As shown in Figure 1a, the parallelogram is formed by an angle θ 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: where K is the permeability of the isotropic medium; µ, p, and ρ are the dynamic viscosity, pressure, and density of the fluid, respectively; and g is the gravitational constant. The stream function Ψ is introduced to relate the flux as given by: where T β is the thermal expansion coefficient ( 0 T β ≥ ) and ν is the kinematic viscosity of the fluid. Based on Equations (11) and (12), the Darcy-modified Rayleigh number T Ra can be defined as: which is also the non-dimensional parameter.

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 0 ψ = . As such, the non-dimensional components of velocity x q and y q are zero at boundary walls, while they can be evaluated by using Equations (5) and (6). The explicit boundary conditions for T and T n ∂ ∂ together with 0 ψ = are described as follows: T y (14) , 0, For the sake of convenience, Equations (1) and (2) are normalized as: where the non-dimensional parameters x = X H , y = Y H , and ψ = Ψ ρα are used. α 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 nondimensional temperature is given by: 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: where a steady state is assumed in the above two equations. Introducing Equations (5) and (6) to Equation (9) gives: For natural heat convection in a non-porous medium, Rayleigh number Ra and Darcy number Da are: where β T is the thermal expansion coefficient (β T ≥ 0) and ν is the kinematic viscosity of the fluid. Based on Equations (11) and (12), the Darcy-modified Rayleigh number Ra T can be defined as: which is also the non-dimensional parameter.

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 ψ = 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 ∂T/∂n together with ψ = 0 are described as follows: ψ cos θ H/L , y = 0, Along the inclined bottom side: Along the inclined upper side: In Equations (16) and (17), n denotes the non-dimensional normal to the inclined wall.

Reproducing Kernel Approximation for Two-Phase Variables
To efficiently solve the nonlinear system of equations, the RK shape functions are introduced for two variables ψ and T at N s source points as follows: where a ψ I and a T I are the generalized coefficients. Ψ I (x) is the RK shape function consisting of a monomial basis H(x − x I ), coefficient vector b(x), and kernel function ϕ a (x − x I ) of the following form: where b(x) is found by the reproducing conditions listed below: where α is the multi-index defined as |α| = α 1 + α 2 , and p is the order of H(x − x I ). The detailed derivation of RK shape functions is referred to previous work [19]. The explicit form of RK shape functions is: with the moment matrix M(x) given by: 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: where the normalized nodal distance s = x − x I /a, and a denotes the radius of a circular support. For the pth order monomial basis, the RK support size is chosen to be proportional to the average nodal distance h, i.e., a = (p + δ)h with δ > 0 for any discretization. The two equations in Equation (18) can be combined as: with: Substituting Equation (24) into Equations (8) and (10) leads to:  (14)- (17) gives the following equations: along x = 0, and: along x = cos θ H/L , and: for the inclined bottom side, and: for the inclined upper side.

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 . . . . . . N p ) for governing equations in the domain, x q (q = 1 . . . . . . N q ) for Neumann boundary conditions, and x r (r = 1 . . . . . . 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: where: with: 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, ψ n , T n , ψ, and T are assumed in the following derivation. Therefore: 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: in which the superscript k denotes the step of iteration and J is the Jacobian matrix defined as J = ∂A ∂a . Thus, the Jacobian matrix is derived as: in which: with: and: 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]: where ∆a k denotes the increment of vector a k at kth step. The explicit expressions for matrices A and J are detailed in Appendix A. The stopping criterion for iteration is given by: which is adopted in the present study.

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].

Numerical Examples
To validate the proposed method in solving the nonlinear system coupled by twofield 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 3 a h = 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.

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  Figure 4 compares the number of iterations obtained by three refinements; Figure 5 shows the corresponding contour plots of ψ and T . Clearly, the present method converges fast and stably in 7 iterative steps for all discretizations; the contour plots of both ψ and T also exhibit the same pattern. The convergence of approximation is validated through the refinement of discretization. In the following study,  Figure 7, where a similar distribution of ψ 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.5h to 4h . 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 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 ( ) 0 0 ,T ψ 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. By taking the boundary conditions described in Equations (14)- (17) or Equations (30)-(33) as an example, the total collocation points are 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.

Numerical Examples
To validate the proposed method in solving the nonlinear system coupled by twofield 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.

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 Ra T = 100 with initial conditions ψ 0 , T 0 = (0, 0.5). Figure 4 compares the number of iterations obtained by three refinements; Figure 5 shows the corresponding contour plots of ψ and T. Clearly, the present method converges fast and stably in 7 iterative steps for all discretizations; the contour plots of both ψ 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 Ra T and ψ 0 , T 0 , 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 ψ and T are shown in Figure 7, where a similar distribution of ψ and T to those obtained by uniform discretization is observed.    horizontal directions for 1 T Ra = ; as T Ra increases to 60 and 100, the contours of ψ gradually showed anti-symmetric distribution along the line inclined at 45  to the horizontal direction as referred to Figures 5b and 12b. As for the contour of T , with increasing T Ra , its distribution moved from a vertical direction toward a curved orientation as shown in Figures 5e and 12c,d.
(a) (b)        Based on the quadratic basis adopted in constructing the RK shape functions, the RK support size generally ranges from 2.5 h to 4 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 Ra T = 100 with ψ 0 , T 0 = (0, 0.5). 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 ψ 0 , T 0 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 ψ 0 le f t , ψ 0 right , T 0 le f t , T 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 Ra T governs the nonlinearity, and various values of Ra T such as 1 and 60 are considered herein. The comparison of the number of iterations obtained by various Ra T is made in Figure 11; it was found that more iteration steps are needed for a larger value of Ra T . For illustration purpose, the contour plots of ψ and T obtained by Ra T = 1 and Ra T = 60 are presented in Figure 12. From Figure 12a, the contour of ψ is almost symmetric in both vertical and horizontal directions for Ra T = 1; as Ra T increases to 60 and 100, the contours of ψ gradually showed anti-symmetric distribution along the line inclined at 45 • to the horizontal direction as referred to Figures 5b and 12b. As for the contour of T, with increasing Ra T , its distribution moved from a vertical direction toward a curved orientation as shown in Figures 5e and 12c,d.

Diamond Domain
A unit diamond with 30 θ =  is discretized uniformly and non-uniformly as shown in the schematic diagrams Figure 13a Figure 14; the corresponding contour plots of ψ and

Diamond Domain
A unit diamond with 30 θ =  is discretized uniformly and non-uniformly as shown in the schematic diagrams Figure 13a Figure 14; the corresponding contour plots of ψ and

Diamond Domain
A unit diamond with θ = 30 • 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 Ra T = 100 and ψ 0 , T 0 = (0, 0.5), 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 ψ 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 ψ and T are shown in Figure 17, where similar patterns of ψ and T to those obtained by uniform discretization are retrieved, while the smoothness of the contours is reduced.
sults are shown in Figures 18 and 19, with Figure 19a illustrating the arrangement of ( ) Figures 18 and 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 T Ra , Figure 20 compares the number of iterations obtained by various T Ra ranging from 1 to 100; the larger value T Ra is, the more iterations it requires. The selected contour plots of ψ and T obtained by 1 T Ra = and 60 T Ra = are shown in Figure 21.
As T Ra increases, the contour of ψ changes from symmetry to anti-symmetry with respect to the line inclined at 45  in the horizontal direction, and its magnitude increases as well, as shown in Figures 15b and 21a,b. The contour of T gradually changes the orientation from left to right with increasing nonlinearity, as can be observed in Figures 15e  and 21c,d.       The influence of initial conditions are investigated by considering ψ 0 , T 0 for entire domain and ψ 0 le f t , ψ 0 right , T 0 le f t , T 0 right for two half domains, and the corresponding results are shown in Figures 18 and 19, with Figure 19a illustrating the arrangement of Figures 18 and 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 Ra T , Figure 20 compares the number of iterations obtained by various Ra T ranging from 1 to 100; the larger value Ra T is, the more iterations it requires. The selected contour plots of ψ and T obtained by Ra T = 1 and Ra T = 60 are shown in Figure 21. As Ra T increases, the contour of ψ changes from symmetry to anti-symmetry with respect to the line inclined at 45 • in the horizontal direction, and its magnitude increases as well, as shown in Figures  15b and 21a,b. The contour of T gradually changes the orientation from left to right with increasing nonlinearity, as can be observed in Figures 15e and 21c,d.

Parallelogram Domains I and II
As shown in Figure 22a,b, two parallelograms (L/H = 4) with inclined angles 30 • and −30 • 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 Ra T = 100 and ψ 0 , T 0 = (0, 0.5), three refinements of uniform discretization with N c = 10 × 40, 20 × 50, 30 × 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 ψ and T obtained by N c = 20 × 50 are shown in Figure 24, which are consistent with the results given in [7]. The distributions of both ψ and T in two domains are dissimilar due to different inclination; in particular, the contour of ψ in domain II with a negative inclined angle has two small vortices at the two ends.

Parallelogram Domains I and II
As shown in Figure 22a Figure 24, which are consistent with the results given in [7].
The distributions of both ψ and T in two domains are dissimilar due to different inclination; in particular, the contour of ψ 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  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 3 s = is presented in Figure 30, and the corresponding contour plots of ψ and T are given in Figure 31. Similar trends in the distribution of both contours ψ 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 4 s = is further investigated, as shown in Figure 32, which has not been discussed in the literature. It is found that the RK support size   The influence of initial conditions on the solutions is investigated by concerning ψ 0 , T 0 and ψ 0 le f t , ψ 0 right , T 0 le f t , T 0 right for Ra T = 100 and N c = 20 × 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 Ra T is presented in Figure 26, where more iterations are needed for larger Ra 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 Ra T = 1, Ra T = 60, and Ra T = 100 are shown in Figures 24, 27 and 28. For both domain I and domain II, the magnitude of ψ increases with larger Ra T , and the order of contour (nonlinearity) of T increases with larger Ra T as well. Specifically, for domain II, the two vortices in the contour of ψ become more pronounced as Ra T increases.    The disturbance of discretization is added to test the stability of the method. By using Ra T = 100 and N c = 20 × 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 ψ and T are given in Figure 31. Similar trends in the distribution of both contours ψ 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 ψ and T are shown in Figure 33 with corresponding iteration given in Figure 34. The patterns of the contours are ensured.

Conclusions
The Newton-Raphson collocation method with reprod is proposed to solve the natural convection problem within a

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 Ra T controls the nonlinearity of the system. In all examples, as the value of Ra T increases, the magnitude of ψ increases, and so does the nonlinearity of contour T. Even so, the present method yields the desired accuracy for a large value of Ra 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 threephase coupling problems, i.e., double-diffusive natural convection problems, is left to be investigated by future works.