Engineering Analysis with Boundary Elements The method of fundamental solutions for the Stokes ﬂow with the subdomain technique

The collocation version of the Method of Fundamental Solutions (MFS) with subdomains is introduced in the present work for the solution of the 2D Stokes ﬂow in backward-facing-step geometry, including Dirichlet and Neumann boundary conditions. The motivation for the present work is the inability of the MFS to solve such problems and the problems with slits and cracks due to the discretization of a single domain. The inability stems from the artiﬁcial boundary that is diﬃcult or impossible to properly geometrically set in such cases. The solution for such problems is found by splitting such domains into subdomains. The MFS equations for the equilibrium conditions at the collocation points on the interface between the adjacent subdomains are derived for the Stokes equation. A matrix that simultaneously solves the collocation problem on all the subdomains is formed and solved. A sensitivity study of the MFS results is performed by comparing the relative root mean square error with the reference solution obtained by the classical mesh-based ﬁnite volume method on a very ﬁne mesh. The subdomain technique is veriﬁed by dividing the domain into 2, 3 and 5 subdomains. The velocity, vorticity and pressure compare very well with the reference solution in all three cases while the solution for the single domain approach is outstandingly poor and inappropriate. The paper shows that the proposed subdomain technique maintains the simplicity, true meshless character and accuracy of the MFS for the Stokes ﬂow in cases where the domain topology requires the use of the subdomain technique.


Introduction
The Stokes flow is a type of steady flow where body forces can be disregarded and the inertial forces are negligible compared to the viscous forces [1] . This can happen when the movement of the fluid is very slow (creeping flow), the characteristic length scale is very small (microfluidic flow) or the fluid dynamic viscosity is very high (e.g., honey). The equations that describe such types of flows represent only the diffusive part of the steady Navier-Stokes equations. The resulting Stokes flow equations have been solved in the past using many numerical techniques such as the following: domain discretization methods (finite difference method [2] , finite element method [3] , and control volume method [4] ), boundary discretization methods (boundary element method [5][6][7] ) and meshless methods (radial basis functions [8] , method of fundamental solutions [9 , 10] , and method of regularized sources [11][12][13] ), among others.
We focus on the implementation of the method of fundamental solutions (MFS) for the Stokes flow in special geometrically nontrivial situations in the present paper. The main advantage of the MFS is that the approach does not use a computational mesh (polygonization in 2D), numerical integration or evaluation of singular functions. It is a truly system matrix are of the same order of magnitude, and the matrix becomes ill conditioned. In either of these two extreme cases, the accuracy of the results is compromised. In the present implementation of the MFS, we place the source points in the outward normal direction on the physical boundary at a distance of 6.5 times the spacing of the collocation points on the physical boundary [14] . It is true that more sophisticated methods have been recently developed for placing the source points; however, they require additional computational effort [21][22][23][24][25] . The arrangement of source points though the defined simple strategy can lead to poor or useless solutions for some geometries of the computational domain, such as sharp corners, slits or cracks, among others. The main idea of this paper is to divide such problematic domains into two or more regular subdomains. This is a strategy that overcomes the identified problem in an ordered way. In addition, the coupling of different subdomains with different material properties is a first step in the discretization of the dispersed multiphase Stokes flow, as appears in femtosecond crystallography.
The division of domains into subdomains, known as the subdomain technique, and the application of interface conditions have been successfully used by various authors. Ravnik et al. [26] used this approach for the solution of vorticity and kinematic equations in the case of a 3D laminar viscous flow using the boundary element method. Panzer [27] modeled the typical vented loudspeaker enclosure using the boundary element method and subdomain technique. The MFS and domain decomposition (subdomain technique) for the Laplace equation in the case of surface seepage flow in layered soil were successfully applied in [28 , 29] . We have recently successfully implemented the domain decomposition approach in connection with the modified method of regularized sources for potential flow [30] .
There are two approaches to domain decomposition: with overlapping subdomains and without overlapping subdomains. For overlapping subdomains, the problem is solved iteratively. On the interface between adjacent subdomains, we consider the values of the variables from the adjacent (overlapping) subdomain as a boundary condition. This approach makes sense when the leading equations are nonlinear, e.g., Navier-Stokes equations. Additionally, this approach is suitable for parallelizing code with the message passing interface (MPI) protocol. In the case of a linear governing equation, e.g., Laplace or Stokes equations, with overlapping subdomains and iterative solving, we would effectively nonlinearize the problem. Therefore, it is better to use the nonoverlapping subdomains and apply interface condition equations to collocation points that are common to adjacent subdomains. In this case, the system of linear equations becomes sparse, and one of the dedicated direct solvers can be used to solve it.
The remainder of this paper is organized as follows. In Section 2 , the problem and the solution method are described. In Section 3 , the backward-facing-step flow is used as a representative numerical example. In Section 4 , a comparison of the results from the proposed method and benchmark solution is discussed. The concluding remarks and foreseen developments are given in the last section.

Governing equations
Consider a fixed two-dimensional domain Ω with boundary Γ filled with a fluid whose steady velocity and pressure fields are described by the following Stokes system of equations: where v = ( v x ,v y ) represents a 2D velocity vector, P represents the pressure field and is the dynamic viscosity. A two-dimensional Cartesian coordinate system with base vectors i x and i y and coordinates p x and p y is introduced. Eqs. (1) and (2) are subsequently rewritten as − + − + To solve the system of Eqs. (1) and (2) in the domain Ω, the boundary conditions of the Dirichlet type = ̄ on Γ D or Neumann type ∕ = ̄ on Γ N have to be set, where the domain boundary Γ is structured into unnecessary connected parts of Γ D and Γ N ; however, Γ = Γ D ∪Γ N .
According to the MFS formulation of the Stokes flow [31] , the velocity and pressure at point p = i x p x + i y p y can be expressed as a sum of N coefficients and , which are not known a priori, multiplied by the point forces (Stokeslets) placed at N source points = + where = 1 , 2 , ..., : The 2D Stokeslets, i.e., the fundamental solutions of the Stokes equation, are [32] : , The Stokeslets derivatives required for the interface condition equations at the subdomain boundaries of the domain are: *  6) and (7) are determined by collocating the solution of the problem with the known values of Dirichlet or Neuman boundary conditions at the collocation points p n where n = 1, 2, .., N . This results in the solution of the 2 N × 2 N system of linear equations, which is solved by the DGESV routine from the LAPACK software package [33] , Intel implementation (MKL).
After the coefficients are known, the velocity (6)-(7) and the velocity derivatives at any location p can be obtained explicitly: These derivatives can also be used for the direct calculation of the vorticity = ∕ − ∕ or for checking the numerically calculated value of the incompressibility condition (1).

Problematic geometry layout
In the MFS, some domain topologies might have a problem with the location of the source points that form the artificial boundary. An example of such a problematic geometry is shown in Fig. 1 , where a slit with a width of the same order as the distance between the collocation points in the domain is present. If we follow the recommendation [21 , 22] of the optimal distance between the source and the collocation point r n , then the source points are positioned inside the domain, although they should be outside of the domain. The solution of the problem with the source point inside the domain results in the erroneous singularity of the  solution at the position of the source. The same problem occurs when cracks (slits with zero thickness) are present.
Another example of problematic geometry is shown in Fig. 2 and represents the so-called l -shaped boundary where two adjacent straight edges form an internal angle of 270°In this case, some of the source points s n , which belong to collocation points p n on the vertical edge and the collocation points p n on the horizontal edge, lie on the same line and vice versa. If the point distribution on edges is uniform, the source points s i and s j can even share the same location. As we will see later, the MFS solution in such situations is not acceptable.

Subdomains -interface conditions
Cases with problematic geometry can be resolved in the framework of the MFS by using the subdomain technique. This means that the computational domain is subdivided into M nonoverlapping subdomains Ω = Ω 1 ∪Ω 2 ∪…Ω M . By doing so, interfaces between these subdomains appear, denoted by Γ I and Γ between adjacent subdomains I and J , respectively. The situation is depicted in Fig. 3 . At each collocation point on the interface, four additional unknown coefficients appear in Eqs. (6)-(8) : and for point p i , which belong to subdomain I ; and and for point p j , which belong to subdomain J . To close the system of linear equations, we need to add four more equations for each collocation point on the interface. These equations represent interface conditions [34] : -equality of velocity (kinematic condition) -equality of normal stress -equality of tangential stress The symmetric deformation tensor is The unit normal n and unit tangent t vectors at the interfacial collocation point are: n = i x n x + i y n y and t = − i x n y + i y n x , respectively. On the interface, the normal and tangent for subdomain J are opposite to the normal and tangent for subdomain I .
The MFS formulations of the interface conditions (25)-(28) are: where N I and N J are the number of collocation points of the I th and J th adjacent subdomains, respectively; and p i ≡ p j . The system matrix is fully populated in the case of a single domain. In the case of division into two subdomains, zeros appear in the matrix, as shown in Fig. 4 . On the diagonal are two nonzero blocks representing contributions from Stokeslets for the collocation points of each subdomain. Below the 1 st block are the rows belonging to the interface conditions (30)-(31) while above the 2 nd block are the rows arising from conditions (32)- (33). As the number of subdomains increases, the matrix becomes increasingly sparse, and a special solver for band matrices can be used instead of a direct Gauss solver.

Numerical example 1
The backward-facing-step flow is one of the standard test cases for the validation of numerical schemes [35][36][37] . A 2D channel of height h with a sudden unilateral extension to a height of 2 h is considered, as shown in Fig. 5 . The corner with a 270°internal angle at point (0, 0) poses a major challenge for any numerical method, particularly the MFS and boundary element methods.    lutions from the analytical solutions was on the order of 10 − 4 m/s, from which we conclude that the Stokes current is fully developed at location p x = 4 h .
-At channel walls, the no-slip condition v x = v y = 0 is prescribed. Here, the step height is h = 1 m, the average inlet velocity is U = 2 m/s and the dynamic viscosity is = 1 Pa ⋅s. A comparison of the results is made at points lying on the vertical cross-section at the location ∕ ℎ = 0 . 1 immediately behind the step, where we expect large changes in the velocity and pressure due to the sudden expansion of the channel.

Reference solution
The control volume method (CVM) solver [38] on a very fine mesh with 1,440,000 square cells of size Δ∕ ℎ = 0 . 0025 was employed as a reference solution with which we compare the MFS results. Because the     CVM code has no specialized Stokes solver, we set the fluid density = 10 − 20 kg/m 3 to minimize the contribution of the convective term ( v ⋅ ∇ ) v in the 2D Navier-Stokes equations. Fig. 6 shows the velocity field, and Fig. 7 shows the normalized pressure contours, where ̄ in and ̄ out are average pressures at the inlet and the outlet boundaries, respectively.

Single domain
The backward-facing-step flow case is first solved with the MFS using a single domain, denoted as MFS 1 . The number of boundary collocation points with equidistant spacing Δ∕ ℎ = 0 . 025 is N = 560. The velocity, pressure and vorticity are evaluated at 3462 interior points, which form a 100 × 40 uniform arrangement. The interior points coincide with the cell center nodes of the CVM mesh.
As Fig. 8 shows, the MFS 1 solution is very poor, especially near the step and corners. The solution for the pressure field, as shown in Fig. 9 , is even worse. Large pressure oscillations at the boundary can be observed with maximum pressure P max = 114,643,347 Pa at the outlet  boundary and minimum pressure P min = 114,563,794 Pa at the inlet boundary, which is physically unrealistic for the given boundary conditions. The number of conditions of the system matrix for this case is Cond(A) = 7.6 × 10 18 . The reason for such a poor performance of the MFS in the configuration shown is explained in Section 2.

Subdomains
Since the MFS 1 solution was inaccurate, we employ the MFS with the subdomain technique, denoted as MFS S , to solve the problem. The Table 3 Maximum absolute divergence and number of conditions of the system matrix.   domain is divided into N s = 2, 3 and 5 subdomains. Fig. 10 shows the collocation points for the cases with the domain subdivision. We must emphasize that each collocation point on the interface has only one adjacent subdomain, as we can see in detail in the picture with 5 subdomains, where subdomains 2, 3, 4 and 5 meet.
The pointization, which is the representation of (sub)domain boundaries with meshless collocation points, was performed by equidistantly distributing points on each subdomain boundary. Table 1 shows the number of collocation points N, where N I represents the interface between the subdomains. Fig. 12 and Fig. 17 .
We performed a sensitivity study with uniform, increasingly dense spacing between collocation points Δ∕ ℎ = 0 . 05 , 0 . 025 , 0 . 0125 and division of the domain into 2 subdomains. The results are compared with the reference solution.
The relative root mean square error (RRMSE) [39] was used as a measure of the deviation of the MFS solution from the reference values. The RRMSE of variable is defined as where K is the number of points, i is the reference solution at point i and num is the numeric (MFS 1 or MFS S ) solution at the same point. The RRMSE of the velocity and vorticity for the considered point spacing for MFS S with 2 subdomains are presented in Table 2 .        The dimensionless velocity and vorticity profiles along the vertical line at ∕ ℎ = 0 . 1 for MFS 1 with 2 subdomains are compared in Fig. 11 -Fig. 13 . The figures show that the RRMSE values and profiles are very close for the medium-and fine-density point arrangement while the solution with the coarse-density point arrangement deviates slightly. A  medium point arrangement is used in further numerical studies in this paper. Fig. 14 and Fig. 15 show the velocity vectors and normalized pressure, respectively, in the case of division of the domain into 5 subdomains. The velocity vectors and pressure contours are the same as in the reference solution in Fig. 6 and Fig. 7 .

Comparison with the reference solution
MFS 1 and MFS S are compared with the reference solution. Table 3 shows the maximum absolute value of the divergence and number of conditions of the system matrix. The table shows that divergence for MFS S solutions is 7-8 orders of magnitude smaller than in the MFS 1 case. Additionally, the number of conditions for the MFS S cases is 5 orders of magnitude lower than that for the MFS 1 case. The computer memory usage and code execution time, as reported by the system, are shown in Table 4 . The table shows that the need for memory increases as the number of subdomains increases; however, the execution time decreases. The reason for the latter is that the system matrix in the MFS becomes increasingly sparse as the number of subdomains increases, despite the increased size. Table 5 shows the pressure difference between the inlet and the outlet. The MFS S results are close to the reference value while the MFS 1 solution is completely wrong. Table 6 shows the RRMSEs for the velocity, vorticity and pressure. The values for MFS 1 are very poor while the errors for the MFS S are small and very close.
The comparison of the normalized v x velocity, v y velocity and vorticity profiles with the reference solution along the vertical cross-section at ∕ ℎ = 0 . 1 are presented in Fig. 16 -Fig. 18 . The MFS 1 profile deviates greatly from the reference profile while the MFS S values practically coincide with the reference profile.

Numerical example 2
The Stokes flow in a 2D channel of height h and length of 4 h , where a "thin wall " obstacle (crack) is placed at a distance of 1.5 h from the inlet boundary serves as a second numerical example, as shown in Fig. 19 . The Dirichlet boundary conditions are prescribed on all boundaries: a fully developed velocity profile = 6 ( ∕ ℎ − 2 ∕ ℎ 2 ) , v y = 0 with an average velocity U = 2 m/s at the inlet and outlet and the no-slip condition v x = v y = 0 at solid walls.
The case is solved with the MFS using a single domain approach and by using a subdomain technique with 3 subdomains. The collocation points are placed on boundaries equidistantly with spacing Δ∕ ℎ = 0 . 05 , as shown in Fig. 20 . Fig. 21 and Fig. 22 show the velocity vectors for the MFS 1 and MFS S solutions. We see that the solution is wrong if we use a single domain approach. Conversely, if we use the subdomain technique, the result for the velocity field is as expected.

Conclusions
In some cases of domain geometry topology, such as slits, cracks or 270°internal angles, the MFS does not provide a useful solution. The well-known related test case represents the backward-facing-step flow. In this particular case, the cause of the problem is the locations of the source points. Since some points lie on the same line as the collocation point, the system matrix is very ill conditioned, so even a small (rounding or truncation) error in the fundamental solution evaluation can cause a large deviation of the results from the correct solution.
The solution of the problem is demonstrated in the present paper in the form of incorporation of the subdomain technique into the original MFS. As shown in the article, by dividing the problematic domain into two or more subdomains, we can present an MFS solution that matches very well with the reference solution. We also see that the velocity profiles, vorticity profiles and RRMSEs with division into 2, 3 or 5 subdomains do not differ significantly. This finding leads to the recommendation that it is best to use the minimum number of subdomains to eliminate the problematic geometry. In the case of backward-facing-step flow, this means 2 subdomains.
However, if we want to use parallelized solvers for sparse matrices due to the structure of the MFS S system matrix, it is better to divide the domain into as many subdomains as possible. This is especially true for 3D geometry, where fully populated system matrices can be huge. The presented subdomain technique is directly applicable in 3D space by replacing the 2D fundamental solution and its derivatives with their 3D equivalents.
Subdomains and the correct implementation of interface conditions are also prerequisites for handling multiphase flows using the MFS. In this case, the viscosity of the fluids and the surface tension force must also be included in the interface conditions. Additionally, the shared edge between the two adjacent subdomains is not fixed but deforms accordingly due to the balance of the forces on either side of the interface.
In the continuation of this work, we intend to apply the represented numerical procedure for efficient meshless simulation of the gas-liquid interface and flow conditions of a two-phase gas-liquid system [40] in the case of a flow-focusing device at the micro scale [41] .

Dedication
The authors dedicate this work to late Professor Frank J. Rizzo, whose pioneering work on fundamental solution-based numerical methods, among other great personalities in the field, inspired the present authors to start to work on the subject in the mid-Eighties; and the authors are still passionate about this beautiful research area.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.