Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 2, pp. 469-483, Warsaw 2014 THE METHOD OF FUNDAMENTAL SOLUTIONS FOR STATIONARY FLOW THROUGH AN AXISYMMETRIC CYLINDRICAL FIBROUS FILTER Jan Adam Kołodziej, Paweł Fritzkowski Poznan University of Technology, Institute of Applied Mechanics, Poznań, Poland e-mail: jan.kolodziej@put.poznan.pl; pawel.fritzkowski@gmail.com A problem of steady-state incompressible fluid flow through a fibrous cylindrical filter is considered. The pressure field is obtained by applying the method of fundamental solutions which gives continuous function in the filter region. The components of filtration velocity are calculated from the appropriate derivatives. In numerical examples, various types of the filter are considered and some computational issues are discussed. A simple algorithm for achieving the optimalpseudo-boundary location is used,within the frameworkof themethod of fundamental solutions, byminimizing the maximum absolute boundary error. Optimiza- tion results for various numbers of source points and collocation points are compared. The variation of total discharge with the inlet size is shown. Keywords: filtration, Darcy equation, method of fundamental solutions, cylindrical fibrous filter 1. Introduction Fibrous filters are one of the cost-effective means used to remove particulate matter suspended in gas stream.Most works on the filtration flow and collection efficiency have involved a single- fiber or microscopic point of view. This paper, by contrast, attempts to consider a filter as a whole. In order to calculate the collection efficiency, detailed information about a flow field in the filter is necessary. From themathematical point of view, whenmodelling a filter as a porous medium, the pressure field or pressure squared field satisfies Laplace equation in the domain andmixed boundary conditions. Now, there aremany numerical methods for solution of such a boundary value problem, e.g. FDM, FEM, FVM. However, these ones are mesh-basedmethods which means that the solution is given only at certain discrete points or is approximated by low-degree polynomials. In some cases, derivatives of the solution are most desirable. During the last two decades, plenty ofmeshlessmethods have been developed and effectively applied to solvemany problems in science and engineering. One of the techniques is themethod of fundamental solutions (MFS), which provides an approximated solution being a continuous function with continuous derivatives. As yet, the MFS has been successfully applied to solve boundary value problems in awide variety of disciplines.Most of authors consider planar doma- ins which lead to particularly easy implementation of themeshless technique. However, in some cases spacial problems can be reduced to the two-dimensional ones. For instance, Karageorghis and Fairweather (1998, 1999, 2000) used the MFS to solve many axisymmetric problems in heat transfer, elasticity and acoustics. Also Ramachandran and Gunjal (2009) presented some axisymmetric heat transfer problems to compare various boundary collocation methods. As it was shown, the fundamental solution of the axisymmetric Laplace equation can be expressed in terms of complete elliptic integrals. In our work, a problem of steady-state incompressible flow through a cylindrical filter is considered. As in many practical applications, it is assumed that the filter is filled with a fi- brous porous medium and Darcy’s law is employed to describe the phenomenon. The Laplace 470 J.A. Kołodziej, P. Fritzkowski type governing equation is solved by means of the MFS, which provides an approximation of the pressure field. Hence, it is easy to obtain velocity of the fluid flow and evaluate the total discharge. The knowledge of the filtration velocity field can be used for calculation of transport of contaminant particles, and their possible removal by the fibers can be investigated (e.g. see Dunnett and Clement, 2006). However, in this paper we restrict our study to determination of the velocity field only. In this paper, the Darcy filtration equation is assumed to be the governing one. It is well known that this equation is justified for low Reynolds numbers. This condition is fulfilled in many cases of the filtration flow, because of low average flow velocity as well as small pore size. An extended discussion on applicability of Darcy’s law can be found e.g. in the paper of Zeng andGrigg (2006). The authors revised two types of criteria, the Reynolds number and the Forchheimer number, and gave their critical values which relate to the situation when the so called ‘non-Darcy effect’ appears and hence Darcy’s law stops being applicable. An outline of this paper is as follows. In Section 2, we present mathematical description of the axisymmetric problem. Section 3 is devoted to the MFS formulation. In Section 4, we de- monstrate some numerical experiments and discuss the results. Finally, conclusions and remarks are given in Section 5. 2. Mathematical formulation 2.1. Boundary value problem Perhaps the simpliest type of a cylindrical fibrous filter comprises a central inlet and outlet (see Fig. 1a). Let us consider the stationary fluid flow through such a filter whose internal geometry forms a cylinder of radius c and height h. Size of the inlet and outlet is specified by the radii a and b, respectively. Assume that the fluid is incompressible. It is exposed to the pressure pin at the inlet and flows past the fibrousmaterial within the filter. Finally, the outlet pressure equals pout. If we treat the fibrous material as an isotropic porous medium, then filtration velocity can be specified according to Darcy’s law in terms of fluid pressure p q=− κ µ gradp (2.1) where κ is the permeability of the medium and µ is the dynamic viscosity of the fluid. Conse- quently, a continuity equation for the incompressible fluid reduces to Laplace’s equation ∇2p=0 (2.2) Due to axial symmetry of the problem, the potential p is independent of the angular co- ordinate ϕ. Thus, the domain can be reduced to a rectangular region Ω with boundary Γ whose revolution about the z-axis could form the original cylinder (see Fig. 1b). Now, governing equation (2.2) takes the form ∂2p ∂r2 + 1 r ∂p ∂r + ∂2p ∂z2 =0 (2.3) whereas the velocity vector has only two non-zero components qr =− κ µ ∂p ∂r qz =− κ µ ∂p ∂z (2.4) The method of fundamental solutions for stationary flow... 471 Fig. 1. Geometry of the filter: (a) 3-D cylindrical domain; the housing is filled with a fibrous porous medium, (b) the domain reduced to a rectangle Fig. 2. Non-dimensional description of the domain and boundary conditions for the problem Using c as a characteristic dimension, one can introduce the following non-dimensional va- riables (see Fig. 2) R= r c Z = z c A= a c B= b c H = h c Moreover, the pressure p can be scaled according to the formula P = p−pout pin−pout so that Pin = 1 and Pout = 0. Now, for the dimensionless pressure field P(R,Z), the velocity is given by QR = cµ κ∆p qr =− ∂P ∂R QZ = cµ κ∆p qz =− ∂P ∂Z (2.5) where ∆p= pin−pout. Equation (2.3), in turn, can be rewritten in the following form ∂2P ∂R2 + 1 R ∂P ∂R + ∂2P ∂Z2 =0 (2.6) As it can be seen in Fig. 2, the boundary Γ = 5 ⋃ m=1 Γm 472 J.A. Kołodziej, P. Fritzkowski and one can specify the boundary conditions as P = f1(R,Z) on Γ1 ∂P ∂Z = g4(R,Z) on Γ4 ∂P ∂Z = g2(R,Z) on Γ2 P = f5(R,Z) on Γ5 ∂P ∂R = g3(R,Z) on Γ3 (2.7) where f1, f5 and g2, g3, g4 are prescribed functions for Dirichlet and Neumann boundary conditions, respectively. In the given case, we have f1 = 1, f5 = 0 and g2 = g3 = g4 = 0. Additionally, the symmetry condition should be taken into consideration, thatmeans zero deri- vative on the axis ∂P ∂R =0 for R=0 (2.8) To sum up, in this axisymmetric problem we seek for the function P(R,Z) which satis- fies partial differential equation (2.6) in the domain Ω, together with conditions (2.7) on the boundary Γ and condition (2.8) on the axis. 2.2. Rate of fluid flow After solving the problem, one can easily obtain the velocity field fromEq. (2.5) and evaluate the total discharge in the next step. It is sufficient to calculate the flow rate at the filter inlet w=2π a ∫ 0 qzr dr (2.9) Hence w=2πc κ∆p µ A ∫ 0 QZRdR=2πc κ∆p µ W where W is dimensionless discharge, which may be specified as W = A ∫ 0 QZRdR (2.10) When dealing with non-dimensional problem (2.6)-(2.8), the permeability of the porous medium is unimportant. However, its role is crucial in the calculation of dimensional velocity field anddischarge evaluation.Thepermeability usually is expressed as some functionof porosity of thefibrousmediumtimes fibersdiameter squared.Let us introduce theporosity term ε, which means a void fraction in the material. In literature the cell models proposed by Happel (1959) or Kuwabara (1959) and the improved ones by other authors (e.g. see Kołodziej et al., 1998), are relatively popular. These models concern a two-dimensional arrangement of parallel fibers. In real media, a three-dimensional arrangement of fibers exists. Assuming that the arrangement of fibres is random and 0.4¬ ε¬ 0.8, one can use the experimental formula provided by Rahli et al. (1996) κ=0.0606d2f π 4 ε5.1 1−ε (2.11) where df denotes the average fiber diameter. The method of fundamental solutions for stationary flow... 473 3. Numerical solution procedure Let us consider boundary value problem (2.6)-(2.8). If P =(R,Z) is a point in Ω, whereas the point Pj =(Rj,Zj) does not belong to Ω, and D2 =(R+Rj) 2+(Z−Zj) 2 k2 = 4RRj D2 then the fundamental solution of axisymmetric Laplace equation (2.6) is given by Φj(R,Z)= 4K(k) D (3.1) and its partial derivatives are determined as ∂Φj ∂R = 2{D2[E(k)− (1−k2)K(k)]−2R(R+Rj)E(k)} RD3(1−k2) ∂Φj ∂Z =− 4(Z−Zj)E(k) D3(1−k2) (3.2) where K(k) and E(k) are the complete elliptic integrals of thefirst and secondkind, respectively K(k)= π/2 ∫ 0 1 √ 1−k2 sin2θ dθ E(k) = π/2 ∫ 0 √ 1−k2 sin2θ dθ (3.3) In the MFS, the solution to the problem is approximated by a linear combination of the fundamental solutions Φj P(R,Z)= N ∑ j=1 cjΦj(R,Z) (3.4) where cj areunknowncoefficients.Thepoints {Pj} N j=1 are singularities (or sourcepoints) located outside the solutiondomain. Inpractice, theyareusuallyplaced on somecontourwhich is similar to the boundary Γ and lies at a distance S from Γ . Additionally, M collocation points Pi are chosen on theboundary.Now, the coefficients cj aredetermined in suchaway that theboundary conditions are satisfied at the collocation points (Kolodziej and Zielinski, 2009, pp. 15-17). Fundamental solution (3.1) ensures the fulfillment of symmetry condition (2.8), thus we use the set of conditions (2.7) only N ∑ j=1 cjΦj(Ri,Zi)= 1 for all (Ri,Zi)∈Γ1 N ∑ j=1 cj ∂ ∂Z Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ2 N ∑ j=1 cj ∂ ∂R Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ3 N ∑ j=1 cj ∂ ∂Z Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ4 N ∑ j=1 cjΦj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ5 (3.5) The number of the collocation points is given by M = 5 ∑ m=1 Mm 474 J.A. Kołodziej, P. Fritzkowski where Mm denotes thenumberof points on Γm.When M =N, linear algebraic system(3.5) can be solvedwith use of theGauss eliminationmethod.Otherwise, if M >N, the over-determined system is solved by the least squares approach. As values of the unknown coefficients are found, one can evaluate the pressure field from (3.4) and the velocity field according to (2.5) and (3.2). 4. Numerical experiments 4.1. Experiment 1 Nowwe turn to application of the presented solution procedure to a test problem.Consider a flowofwater through afilterwhose dimensions are specified as follows: c=40mm,h=140mm, a= b=10mm.Assume the inlet and outlet pressure: pin =400kPa, pout =200kPa.Moreover, for the porous material we take ε=0.6 and df =0.15mm. The supposed viscosity of water is µ=0.001Pa·s. Firstly, we focus on determining the optimal placement of the source points. As mentioned above, theymay lie on some pseudo-boundary, at the distance S from the boundary. In fact, the distance affects the solution accuracy, and choosing right values of S plays akey role in theMFS. According to the concept presentedbyKarageorghis (2009), weminimize the absolutemaximum error eMAX on Γ ; such a technique relies on the maximum principle for harmonic functions. Thus,we choose a set of M∗ boundarypoints {(Rl,Zl)} M∗ l=1, different fromthe collocation points. Next, for consecutive values of S, we evaluate adifferencebetween the solution (or its derivative) and the assumed boundary values el =      P(Rl,Zl)−f(Rl,Zl) for Dirichlet boundary conditions ∂P ∂n (Rl,Zl)−g(Rl,Zl) for Neumann boundary conditions where ∂/∂n denotes the outward normal derivative at the boundary point (Rl,Zl). Finally, the maximum absolute error is given by eMAX = max l=1,...,M∗ |el| (4.1) Thus, one can find the optimal distance Sopt by minimization of the function eMAX(S). As can be seen in Fig. 3, in this problem, the objective function has a pseudo-random appearance beyond the initial range. Consequently, we use a simple search algorithm for the optimization, that means S is increased linearly to calculate the error. Fig. 3. Variation of the error with distance of source points The analyzed results were obtained for various numbers M and N. However, two other symbols are used: MC and NC, which denote the numbers of collocation points and source The method of fundamental solutions for stationary flow... 475 points per unit length C. On each boundary segment, the points are distributed proportionally and uniformly. Figure 4 illustrates the maximum error when keeping MC fixed and varying NC. In each case, we searched for Sopt in the interval 〈0.01,0.7〉 with a step ∆S=0.001. The graphs reveal that the improvement of the accuracy for increasing NC is not so evident. Although the error is relatively high, it seems to stabilize for greater values of MC within the approximate range 0.6¬NC/MC ¬ 0.8. Doubtless, the extreme case, when NC =MC, should be omitted. Fig. 4. Maximum error on the boundary Fig. 5. Dimensionless total discharge Actually, the key question is how the error impacts on values of the total discharge. The next four graphs (Fig. 5) show the non-dimensional quantity W for the given MC and NC. As with eMAX, we observe similar intervals of NC/MC where very small fluctuations of the 476 J.A. Kołodziej, P. Fritzkowski discharge appear. It seems that the high error is local and does not disturb the velocity field considerably. Table 1. Selected results of the error minimization MC NC M N Sopt eMAX eRMS W w [l/min] eW 16 13 89 72 0.051 3.56E-02 1.06E-02 6.21E-02 37.048 2.51E-02 14 78 0.075 2.90E-02 8.16E-03 6.02E-02 35.931 1.09E-02 20 14 111 78 0.074 4.20E-02 9.78E-03 6.05E-02 36.105 6.90E-03 15 83 0.129 4.79E-02 9.98E-03 5.89E-02 35.126 6.60E-03 16 89 0.393 5.28E-02 1.43E-02 5.95E-02 35.489 2.80E-03 17 94 0.039 3.33E-02 9.51E-03 6.32E-02 37.693 3.03E-02 30 22 166 122 0.043 3.98E-02 7.61E-03 6.29E-02 37.549 1.04E-02 23 127 0.065 3.93E-02 7.01E-03 6.18E-02 36.849 5.30E-03 24 133 0.078 4.80E-02 8.12E-03 6.11E-02 36.462 1.80E-03 25 138 0.026 3.32E-02 9.14E-03 6.39E-02 38.104 4.28E-02 26 144 0.026 2.30E-02 7.38E-03 6.39E-02 38.121 1.28E-02 40 30 221 166 0.031 4.46E-02 6.91E-03 6.39E-02 38.095 1.04E-02 31 171 0.040 3.07E-02 5.48E-03 6.31E-02 37.649 6.50E-03 32 177 0.055 3.68E-02 5.58E-03 6.24E-02 37.257 2.60E-03 33 182 0.019 3.95E-02 1.02E-02 6.42E-02 38.315 1.30E-03 34 188 0.020 3.26E-02 7.48E-03 6.44E-02 38.401 1.20E-02 35 193 0.035 2.70E-02 4.80E-03 6.33E-02 37.793 8.00E-03 Table 1 presents selected results of the discussed optimization. The examples show how the numbers MC and NC correspond to the total number of collocation and source points (M and N). Apart from the maximum absolute error, the root mean square one is given, according to the formula eRMS = √ √ √ √ 1 M∗ M∗ ∑ l=1 e2l (4.2) As it can be seen, the total discharge w, computed for the given dimensions and properties, fluctuates slightly and its approximate value is 37-38l/min. Additionally, the table includes values of the relative error between the total inflowand outflow. If the former quantity is treated as the reference one, the error can be expressed in the following way eW = |Wout−Win| Win (4.3) where the subscripts “in” and “out” allow one to distinguish between the discharge at the filter inlet and outlet. As the results indicate, eW does not exceed 5%, and inmany cases is even less than 1%. Figure 6, in turn, illustrates the distribution of the non-dimensional pressure and velocity in the Z direction obtained for MC =40, NC =35 and S =Sopt. Also, for the same parameters, we examined variation of the discharge W with the radius A assuming that B is constant (see Fig. 7). Obviously, the discharge value grows due to an increase of the inlet size. As mentioned above, relatively high values of the boundary error seem to be a local effect. Indeed, the rootmean square error is oneorder ofmagnitude smaller than themaximumabsolute error (see Table 1). Presumably, the maximum values of the boundary error appear in the neighborhood of the so called boundary singularities: where the boundary condition suddenly changes from P = f to ∂P/∂n = g. For instance, Fig. 8 shows the error distribution on the The method of fundamental solutions for stationary flow... 477 Fig. 6. Field of pressure (left) and Z-velocity (right) for MC =40,NC =35 and S=Sopt =0.035 Fig. 7. Total discharge as a function of radius A Fig. 8. Absolute error e on the upper boundary: (a) for 0¬R¬A; (b) for A¬R¬C upper boundary (Z =H) of the filter: Γ1∪Γ2. The greatest values occur around the singularity for R=A. Moreover, a high error arises near by the corner (R=C). 4.2. Experiment 2 Let us now consider another type of cylindrical filter: assume that there is a circumferential outlet instead of the central one. We suppose that the filter housing, containing the fibrous material, is constructed in such a way that it does not disturb the axial symmetry of the fluid flow through the outlet. The domain of such a problem is illustrated inFig. 9. As can be seen, in this case, B denotes the outlet half-width. Also, the system of equations (3.5) should be slightly reformulated due to a change in boundary conditions (2.7) 478 J.A. Kołodziej, P. Fritzkowski Fig. 9. The domain and boundary conditions for the filter with a circumferential outlet N ∑ j=1 cjΦj(Ri,Zi)= 1 for all (Ri,Zi)∈Γ1 N ∑ j=1 cj ∂ ∂Z Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ2 N ∑ j=1 cj ∂ ∂R Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ3 N ∑ j=1 cjΦj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ4 N ∑ j=1 cj ∂ ∂Z Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ5 For computations, assume the following dimensions: c = 45mm, h = 200mm, a = 15mm, b=20mm.Moreover, we take pin =400kPa, pout =250kPa and µ=0.001Pa·s. For the porous material ε=0.6 and df =0.15mm. Similarly to the first experiment, we applied the optimization algorithm for various values of MC and NC, when S ∈ 〈0.01,0.7〉 and ∆S=0.001. Table 2 presents selected results. Accor- dingly, the total discharge w is approximately equal to 48l/min. Again, its fluctuations comes from locally high boundary error, particularly near the singularity at R=A. Theminimization procedure allows one to reduce eMAX below 3%.The discharge error, eW in turn, reaches 1-2%; only for the lowest MC the error is greater than 5%. Table 2. Selected results of the error minimization MC NC M N Sopt eMAX eRMS W w [l/min] eW 16 13 104 84 0.051 2.90E-02 8.01E-03 9.36E-02 47.125 1.68E-02 14 91 0.089 3.95E-02 7.46E-03 9.08E-02 45.712 5.49E-02 20 16 130 104 0.038 3.18E-02 9.55E-03 9.50E-02 47.811 2.10E-03 17 110 0.058 3.30E-02 6.82E-03 9.30E-02 46.823 1.85E-02 23 149 0.051 3.16E-02 5.06E-03 9.48E-02 47.698 9.10E-03 30 25 194 162 0.027 3.66E-02 7.24E-03 9.62E-02 48.434 5.80E-03 26 168 0.026 2.87E-02 6.95E-03 9.60E-02 48.303 1.19E-02 30 194 0.064 3.78E-02 4.05E-03 9.47E-02 47.678 1.52E-02 40 34 259 220 0.020 3.24E-02 6.64E-03 9.66E-02 48.619 5.10E-03 35 226 0.020 2.70E-02 6.11E-03 9.69E-02 48.777 1.87E-02 Figure 10 shows distribution of the non-dimensional pressure and axial velocity computed for MC = 40, NC = 35 and S = Sopt. The change with the radius A in the discharge W has the same character as before (see Fig. 11). The method of fundamental solutions for stationary flow... 479 Fig. 10. Field of pressure (left) and Z-velocity (right) for MC =40,NC =35 and S=Sopt =0.02 Fig. 11. Total discharge as a function of radius A 4.3. Experiment 3 Finally, consider a filter with both circumferential inlet and outlet. Herein, one should also assume that the structure of the filter housing does not disturb the radial flow through such an inlet and outlet. The resulting domain of the filtration problem is presented in Fig. 12. Here, B denotes the outlet half-width, whereas A is the inlet half-width. Because of the different geometry, the system of equations (3.5) takes the form N ∑ j=1 cj ∂ ∂Z Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ1 N ∑ j=1 cjΦj(Ri,Zi)= 1 for all (Ri,Zi)∈Γ2 N ∑ j=1 cj ∂ ∂R Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ3 N ∑ j=1 cjΦj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ4 N ∑ j=1 cj ∂ ∂Z Φj(Ri,Zi)= 0 for all (Ri,Zi)∈Γ5 What is more, the formulas related to the total dischargemust bemodified. For the circum- ferential inlet we have 480 J.A. Kołodziej, P. Fritzkowski w=2πc h ∫ h−2a qr dz=2πc κ∆p µ H ∫ H−2A QR dZ =2πc κ∆p µ W (4.4) where the dimensionless discharge W is given by W = H ∫ H−2A QR dZ (4.5) Fig. 12. The domain and boundary conditions for the filter with a circumferential outlet In this example, suppose the dimensions: c=45mm, h=200mm, a=25mm, b=20mm. Moreover, assume that pin = 400kPa, pout = 200kPa, µ = 0.001Pa·s, ε = 0.5 and df =0.15mm. The making use of the optimization procedure (S ∈ 〈0.01,0.7〉 and ∆S =0.001) produced the results shown in Table 3. As can be seen, the dimensional discharge w ≈ 35l/min. In this case, the highest values of the boundary error occur at two singular points: at Z =0and Z =H as R=C. Taking S = Sopt can reduce the maximum error below 2%. Furthermore, it should be noticed that the error eW reaches smaller values than in the two previous experiments. Table 3. Selected results of the error minimization MC NC M N Sopt eMAX eRMS W w [l/min] eW 16 13 105 84 0.493 1.96E-02 4.72E-03 1.64E-01 34.662 4.00E-03 14 91 0.477 1.87E-02 6.06E-03 1.66E-01 35.217 1.00E-04 20 14 130 91 0.443 1.87E-02 5.24E-03 1.64E-01 34.763 2.90E-03 15 97 0.466 1.72E-02 4.68E-03 1.65E-01 34.846 1.04E-02 22 142 0.321 1.70E-02 3.74E-03 1.66E-01 35.260 3.00E-04 30 23 195 149 0.275 1.49E-02 4.64E-03 1.66E-01 35.153 3.00E-03 24 155 0.306 1.64E-02 4.17E-03 1.66E-01 35.254 7.80E-03 31 200 0.207 1.75E-02 3.75E-03 1.66E-01 35.217 1.60E-03 40 32 259 207 0.200 1.74E-02 3.44E-03 1.65E-01 34.956 3.20E-03 34 220 0.201 1.52E-02 3.12E-03 1.66E-01 35.126 6.40E-03 The distribution of the pressure P and the velocity QZ is presented in Fig. 13. The W(A) function, in turn, is shown inFig. 14. One can observe that the dependency is of different nature for the circumferential inlet and outlet. The method of fundamental solutions for stationary flow... 481 Fig. 13. Field of pressure (left) and Z-velocity (right) for MC =40,NC =32 and S=Sopt =0.2 Fig. 14. Total discharge as a function of radius A 5. Final remarks In this work, we have presented the application of the method of fundamental solutions to steady-state fluid flow through a cylindrical filter. Due to axial symmetry of the domain and boundary conditions, the problem has been reduced to a two-dimensional one. This formulation has led to very easy implementation of the method, which should be emphasized. In this paper,we assumed that the fluid is incompressible. For a gas flow through a filter, one can take into account the compressibility. The continuity equation for a steady axisymmetric compressible flow has the form 1 r ∂ ∂r (rρqr)+ ∂ ∂z (ρqz)= 0 (5.1) where ρ is density. According to the equation of state for the perfect gas ρ = p/(RT) and Darcy’s law for isothermal condition, the continuity equation takes the form 1 r ∂ ∂r ( rp ∂p ∂r ) + ∂ ∂z ( p ∂p ∂z ) =0 (5.2) or the equivalent form 1 r ∂ ∂r ( r ∂p2 ∂r ) + ∂ ∂z (∂p2 ∂z ) =0 (5.3) 482 J.A. Kołodziej, P. Fritzkowski The last equation has the same form as equation (2.3), if instead of pressure p, we have pres- sure squared p2. This equation must be solved with the same boundary conditions as for the incompressible flow (see e.g. Uściłowska and Kołodziej, 2006). In such a case, the methodology of solution of the nondimensional problem is practically the same as for the incompressible case. Differences essentially exist in calculation of dimensional parameters. The presented numerical examples have been related to various types of the cylindrical filter: diverse positions of the inlet and outlet have been considered. In all the cases, we used a simple algorithm for achieving the optimal pseudo-boundary location by minimizing the maximum absolute boundary error. The optimization results for various number of source points and collocation points havebeencompared. Inour experiments, themaximumerror is relatively high, however, it turnsout that it doesnot significantly affect thepressureandvelocity distribution.As ithasbeen shown, themaximumerror appears in theneighborhoodof theboundarysingularities, due to a sudden change in boundary conditions. Since a non-zero difference between total inflow and outflow is a purely numerical effect, it can be also treated as a measure of the solution quality. The obtained results indicate that in most cases the discharge error is lower than 2% and has the smallest values as both the inlet and outlet of the filter are circumferential. Moreover, we examined the variation of total discharge with the inlet size. The function is qualitatively different when dealing with the filter comprising a circumferential inlet and outlet. All in all, with the typical advantages of the meshless methods (e.g. no discretization of the domain), themethod of fundamental solutions allows one to obtain reasonable results: the field of pressure and velocity as well as the total discharge within acceptable tolerance. References 1. Dunnett S.J., Clement C.F., 2006, A numerical study of the effects of loading from diffusive deposition on the efficiency of fibrous filters,Aerosol Science, 37, 1116-1139 2. Happel J., 1959, Viscous flow relative to arrays of cylinders,AIChE Journal, 5, 174-177 3. Karageorghis A., 2009, A Practical Algorithm for Determining the Optimal Pseudo-Boundary in the Method of Fundamental Solutions, Advances in Applied Mathematics and Mechanics, 1, 510-528 4. Karageorghis A., Fairweather G., 1998, The method of fundamental solutions for axisym- metric acoustic scattering and radiation problems, Journal Acoustical Society of America, 104, 3212-3218 5. Karageorghis A., Fairweather G., 1999, The Method of Fundamental Solutions for axisym- metric potential problems, International Journal for Numerical Methods in Engineering, 44, 1653- 1669 6. Karageorghis A., Fairweather G., 2000, The Method of Fundamental Solutions for axisym- metric elasticity problems,Computational Mechanics, 25, 524-532 7. Kołodziej J.A., Dzięcielak R., Kończak Z., 1998, Permeability tensor for heterogeneous porous medium of fibrous type,Transport in Porous Media, 32, 1-19 8. Kolodziej J.A, Zielinski A.P., 2009, Boundary Collocation Techniques and their Application in Engineering, WIT Press, Southampton 9. Kuwabara S., 1959, The forces experienced by randomly distributed parallel circular cylinders or spheres in viscous flow at small Reynolds numbers, Journal of Physical Society of Japan, 14, 527-532 10. Rahli O., Tadrist L., Miscevic M., 1996, Experimental analysis of fibrous porous media per- meability,AIChE Journal, 42, 3547-3549 The method of fundamental solutions for stationary flow... 483 11. RamachandranP.A.,GunjalP.R., 2009,Comparisonof boundary collocationmethods for sin- gular and non-singular axisymmetric heat transfer problems,Engineering Analysis with Boundary Elements, 33, 704-716 12. Uściłowska A., Kołodziej J.A., 2006, Solution of the nonlinear equation for isothermal gas flows inporousmediumbyTrefftzmethod,ComputerAssistedMechanics andEngineering Science, 13, 445-456 13. ZengZ-W.,GriggR., 2006,A criterion for non-Darcyflow inporousmedia,Transport in Porous Media, 63, 57-69 Manuscript received September 23, 2011; accepted for print November 8, 2013