An implicit meshless scheme for the solution of transient non-linear Poisson-type equations

Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

Contents lists available at SciVerse ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

An implicit meshless scheme for the solution of transient non-linear Poisson-type equations G.C. Bourantas a, V.N. Burganos b,n a Applied Mathematics and Computational Science and Earth and Environmental Sciences and Engineering, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia b Institute of Chemical Engineering Sciences—Foundation for Research and Technology, Stadiou, Platani, Patras 26504, Greece

art ic l e i nf o

a b s t r a c t

Article history: Received 2 November 2012 Accepted 4 April 2013 Available online 23 May 2013

A meshfree point collocation method is used for the numerical simulation of both transient and steady state non-linear Poisson-type partial differential equations. Particular emphasis is placed on the application of the linearization method with special attention to the lagging of coefficients method and the Newton linearization method. The localized form of the Moving Least Squares (MLS) approximation is employed for the construction of the shape functions, in conjunction with the general framework of the point collocation method. Computations are performed for regular nodal distributions, stressing the positivity conditions that make the resulting system stable and convergent. The accuracy and the stability of the proposed scheme are demonstrated through representative and well-established benchmark problems. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Meshfree point collocation method MLS Nonlinear Poisson equation Linearization method Lagging coefficients

1. Introduction Non-linearity is inherent in the majority of physical phenomena and engineering processes, resulting in non-linear ordinary or partial differential equations. In fact, these non-linear equations are usually difficult to solve, since no general technique works globally and, thus, each individual equation has to be studied as a separate problem. Of particular interest are the non-linear Poisson-type equations that are encountered in transport problems, notably heat conduction, mass transport and flow in different geometries, including pore structures. Solutions to this type of problems are usually required in non-regular two- and three-dimensional geometries with nonuniform boundary conditions. The non-linear character of the partial differential equations that govern these problems renders the analytical solution either difficult or impossible, especially when dealing with irregular geometries and complex boundary conditions. The use of conventional numerical procedures, such as finite differences [1], finite elements [2], boundary element method [3], etc., to solve these problems necessitates high levels of discretization, resulting in large computational time. In fact, the difficulty of the finite difference method to deal with problems with irregular geometry is a major drawback of the method. On the other hand, the finite

n

Corresponding author. Tel.: +30 2610 965215. E-mail address: [email protected] (V.N. Burganos).

0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.04.003

element methods do handle irregular geometries, yet the refinement procedure becomes a major task. To overcome these drawbacks, new sophisticated numerical methods have appeared over the last years [4,5]. Meshless (or meshfree) methods emerged as potential candidates to replace traditional numerical methods for problems where the latter failed to provide accurate or fast results. Thus, the challenges that these methods have to face are both the efficiency and the accuracy of computations. In a general sense, a numerical method must be robust and accurate in order to be able to provide numerical results in real time when dealing with real world applications. Additionally, they have to deal with the complexity of the phenomena under consideration, which makes the computations even more challenging. In fact, in cases where no analytical solutions are available, the refinement procedure becomes a necessity and, therefore, an easy and efficient refinement procedure must be available. Despite the tremendous improvement of the meshless methods, they are still at their early stages of development and, thus, several test studies must be conducted in order for these methods to establish themselves as more general practical tools ([6] and references there in). Among the existing meshless approximation/interpolation methods, both Moving Least Squares (MLS) and Radial Basis Functions (RBF) have been widely used to solve numerous challenging problems [7]. Meshfree methods based on MLS have a clear advantage over other meshfree methods thanks to the simplicity and stability of the MLS method in field variable approximation/interpolation [4,5]. The meshless methods provide accurate numerical results

1118

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

when dealing with certain Poisson-type problems, linear and nonlinear, in two and three dimensions [8–13]. The governing equations can be solved in their weak form [10–12] or in their strong formulation [8,9,13]. Despite their accuracy, the former methods deal with the computation of integrals and may prove very time demanding and inaccurate, whereas the latter lead to dense matrices that can be difficult to solve when dealing with real world applications with increased number of degrees of freedom. Despite their computational cost, thanks to their key advantage of being applicable to arbitrarily irregular geometries, they are attractive for solving non-linear Poisson type problems. In this work, a meshfree point collocation method is used for the numerical simulation of both transient and steady state nonlinear Poisson-type partial differential equations with the objective to extend the applicability of the method to a wider range of non-linear differential equations. The strong form equations are solved with the collocation method. The implementation is easy and straightforward and, the resulting linear system is sparse and positive definite. Thus, the numerical solution of the harmonic operator can be accurate and fast, taking advantage of the sparsity of the matrix. The article is organized as follows. In Section 2, the governing equations are presented along with a complete description of the MLS approximation scheme. Section 3 presents the meshless point collocation transient, coupled solver along with a θ-weighted time discretization approach, which is employed to solve numerically the partial differential equations. Numerical issues are discussed in Section 4, where typical transient and steady state physical problems are solved and the validity of the proposed meshless techniques is demonstrated. Finally, in Section 5, the main conclusions are discussed.

the boundary, respectively, whereas α, β, ΒR are known coefficients. Eq. (1) applies to a wide range of engineering problems that distinguish themselves through the type of the right-hand side function f. More precisely, Eq. (1) can be reduced to a standard Poisson equation or Helmholtz equation if f is a function of position vector x or a linear function of physical field u, respectively. More generally, f may be a nonlinear function of physical field u and its derivatives. The method that is presented below is able to handle these types of problems. 2.2. Moving Least Squares shape functions In the context of the meshless approximation/interpolation schemes, the MLS method [14] is widely used, since it can directly approximate the field variables in a local manner and, additionally, can be easily extended to n-dimensional problems. A brief summary of the MLS approximation scheme is given next. Within the MLS context, the approximation uh(x) of the unknown field function u(x) is expressed as: m

uh ðxÞ ¼ ∑ pi ðxÞαi ðxÞ ¼ pT ðxÞaðxÞ

ð2Þ

i¼1

where pT(x) is a polynomial basis in the space coordinates, that consists most often of monomials of the lowest order to ensure completeness, m is the total number of the terms in the basis and, α(x) is the vector of coefficients, given by αðxÞ ¼ ðα0 ðxÞ; α1 ðxÞ; :::; αm ðxÞÞT

ð3Þ

which is a function of x. Due to the local nature of the approximation, the polynomial basis can be written as pT ðx−xi Þ ¼ ½1; ðx−xi Þ; ðy−yi Þ; :::; ðx−xi Þðy−yi Þm−1 ; ðy−yi Þm 

2. Mathematical problem and approximation procedure

in 2D problems. Herein, a second order (m ¼ 6) polynomial basis is used,

2.1. Governing equations

pT ðxÞ ¼ ½1; ðx−xi Þ; ðy−yi Þ; ðx−xi Þ2 ; ðx−xi Þðy−yi Þ; ðy−yi Þ2 

Consider the following partial differential equation: ∂u ¼ Δu þ f ðx; u; u;x ; u;y ; u;xx ; u;xy ; u;yy Þ ∂t

ð1Þ

for the physical field u, subject to the following boundary and initial conditions:

 Dirichlet boundary condition for the unknown field: u ¼ uD on

∂ΩD

 Neumann boundary condition for the component of the fieldgradient normal to the boundary: q ¼ qN on

field-gradient normal to the boundary: ∂ΩR

 Initial condition at t¼0: u ¼ u0

on

n

n

i¼1

i¼1

JðxÞ ¼ ∑ wi ðx;xi Þ½uh ðxÞ−uðxÞ2 ¼ ∑ wi ðx;xi Þ½pT ðxi ÞaðxÞ−uðxi Þ2 ð6Þ

 Mixed (Robin) boundary condition for the component of the on

ð5Þ

The derivatives of the unknown field function are up to second order. Thus, in the context of strong form collocation the monomial basis has to be at least second order. For our computations we used the lowest (second) order of monomials that ensured both convergence and accuracy for the harmonic operator [5]. There exists a unique local approximation associated with each point in the domain. To obtain the local approximation of the function u(x) and determine the form of α(x), the difference between the local approximation uh(x) and the function u(x) must be minimized. Thus, a weighted discrete error norm,

∂ΩN

au þ βq ¼ BR

ð4Þ

Ω

where x∈ℝd is the position vector, d ¼2 is the dimension of domain Ω(x,y), which has a piecewise smooth boundary ∂Ω. Δ represents the Laplacian operator, and q stands for the boundary normal flux defined by q ¼ −∂u=∂n, with n the unit outward normal to the boundary ∂Ωð∂Ω ¼ ∂ΩD ∪∂ΩN ∪∂ΩR Þ, uD, qN are specified values on

is constructed and minimized with respect to the vector α(x) of coefficients. The weight function wi ðx;xi Þ≡wðx−xi Þ is usually built in such a way that it takes a unit value in the vicinity of node i, where the function and its derivatives are to be computed, and vanishes outside a region Ωi surrounding the point xi (called the support domain of node i), with n being the number of nodes in the domain. The choice of the weight function wðx−xi Þ affects the resulting approximation uh ðxi Þ significantly. In the present work, a Gaussian weight function is used [5], yet the support domain does not have a standard point density value. Instead, a constant number of nodes are used for the approximation of the field function. That is, 8 9 2 < e−ðDðxÞ=di Þ 0 ≤ DðxÞ ≤1 = di wi ðxÞ ¼ ; ð7Þ DðxÞ : 0 41 ; di

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

where D(x)¼||x−xi|| and di is the size of the support domain with i ¼ 1; 2; 3; :::; q are the nodes that produce the support domain of node xi. Eq. (2) can be rewritten as uh ðxÞ ¼ ðpT ðxÞA−1 ðxÞBðxÞÞU ¼ ΦðxÞU

ð8Þ

where matrices A, B and U are defined as n

AðxÞ ¼ ∑ wi ðx;xi Þpðxi ÞpT ðxi Þ

ð9Þ

i¼1

 BðxÞ ¼ w1 ðxÞpðx1 Þ

w2 ðxÞpðx2 Þ

:::

wn ðxÞpðxn Þ



U ¼ ½u1 ; u2 ; :::; un T

ð10Þ ð11Þ

Finally, the MLS approximant is obtained as n

uh ðxÞ ¼ ∑ Φi ui ¼ ΦðxÞU

ð12Þ

i¼1

where the meshless shape function Φi(x) is defined as m

Φi ðxÞ ¼ ∑ pj ½A−1 ðxÞBðxÞji

ð13Þ

j¼1

The partial derivatives of the shape functions Φi(x) are obtained as m

Φi;k ðxÞ ¼ ∑ ½pj;k ðA−1 BÞji þ pj ðA−1 B;k þ A−1 ;k BÞji  where ðÞ;k denotes ∂ðÞ=∂xk and inverse of the matrix A given by

ð15Þ

following the analysis in [9].

3. Numerical scheme A numerical scheme for the numerical solution of the nonlinear Poisson type equations is presented in this section. The time evolution in the governing equations is monitored at discrete time steps using the one-step θ method (see below for details), while the shape functions and their derivatives, up to second order, are employed using the MLS formulation. The constructed shape functions and their derivatives are used to calculate the numerical solution of the non-linear elliptic PDE for the field variables at discrete time instants, with the nonlinear terms treated iteratively within each time step [9]. The iterative process initiates from an initial guess u0 for the unknown field function within the entire spatial domain Ω(x,y). Subsequently, the right-hand side of Eq. (1) is calculated, based on the initial guess and, then, the resulting linear system is solved using a direct method. The updated un+1 values are then used to calculate the right hand side of Eq. (1), initiating the aforementioned iterative procedure. The relative error of the new and the old value is then calculated and the iterative procedure ends when the value of the relative error is less that a threshold value, typically, 10−10. The application of this collocation scheme to the transient Poisson type equations using the Moving Least Squares approximation is described next. The non-linear partial differential equations are given as: ∂uðx; tÞ þ Luðx; tÞ ¼ f ðu;x ;u;y ; x; tÞ; t 40 ∂t

uðnþ1Þ −uðnÞ ðnþθÞ ¼f −ðθLuðnþ1Þ þ ð1−θÞLuðnÞ Þ δt

ð17Þ

where δt is the time step, n is iteration step, defined as t(n) ¼ nδt, with the non-linear terms f(n+θ)calculated explicitly in the next time step. Thus, the parabolic PDEs are replaced with the semidiscrete PDEs of the elliptic type for the field variables u(n+1) assuming the fields u(n) to be known from the computation at the previous time step. A crucial point in the context of the meshless methods is the nodal configuration that can be uniform or irregular. It has been stated that when the positivity conditions, stated as Φi(xj), ∇2Φi(xj)≥0 for i≠j and ∇2Φi(xj)o0, Φi(xj) being the approximation function of a point i evaluated at a point j, are fulfilled [15], the linearized system is diagonally dominant and, hence, it can be solved with an iterative solver. Additionally, it has been proved [16] that the discrete harmonic operator converges when the Dirichlet boundary conditions are applied. Therefore, an embedded uniform nodal configuration was employed here to ensure both accuracy and convergence. The method proved to be robust and computationally efficient so that it can be used for the solution of problems that are encountered in numerous practical applications.

4. Numerical examples

represents the derivative of the

−1 −1 A−1 ;k ¼ −A A ;k A

time derivatives in the θ method is given as:

ð14Þ

j¼1

A−1 ;k

1119

ð16Þ

with given initial and boundary conditions (Dirichlet and/or Neumenn), where L is a linear differential operator with derivatives up to second order. The finite difference approximation of the

In order to demonstrate both the efficiency and the accuracy of the proposed method, representative benchmark numerical examples of non-linear transient and steady state Poisson type problems are considered and the numerical results are compared with those obtained using traditional and well established numerical methods, or analytical solutions, where available. Examples of nonlinear Poisson type equations are treated numerically here by the Meshless Point Collocation (MPC) method. In the case where the PDE is linear, the solution procedure is straightforward since the resulting system is linear and solvable in a direct manner. On the other hand, if the resulting system is nonlinear, then a linearization procedure must be used. Herein, we mainly adopted the lagging coefficients method, which is considered to be the simplest method of linearization. For the cases where the lagging coefficients method failed to converge and give inaccurate results, we used other efficient linearization methods, such as the Newton method. Solutions to representative examples of the generalized nonlinear Poisson-type equation, namely, the Liouville and Bratu equations, and the Poisson-Boltzmann equation are presented below. Transient problems are considered then, associated with convective diffusion and reaction, and with microwave heating. The extension of the present two-dimensional model to three dimensions is straightforward. The only additional effort that is needed is to recalculate the shape functions for three dimensional nodal distributions, which simply increases the computational cost.

4.1. Steady state problems

Example 1. Generalized non-linear Poisson equation Initially, for the sake of comparison with literature results [17], a non-linear Poisson problem is considered that is defined in a square domain of length L ¼1. The governing equation is [17]   ∂uðx; yÞ 2 ∇2 uðx; yÞ þ ¼ 2y þ x4 ð18Þ ∂y

1120

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

is assumed to be

subject to the following Dirichlet boundary conditions: uð0; yÞ ¼ 0

and

uðx; 0Þ ¼ ax

and

uð1; yÞ ¼ y þ a uðx; 1Þ ¼ xðx þ aÞ

ð19Þ

The analytical solution of Eq. (18) is u(x,y)¼x(xy+α). In the present computations, the constant α is selected to be zero for simplicity. A convergence analysis is undertaken in order to depict both the efficiency and the accuracy of the proposed meshless scheme. The convergence study was conducted using several grid configurations, namely, 4141, 8181, 101101 and 141141, of which the 101101 grid was eventually selected as the optimal one. The initial guess for the iteration procedure was set to zero, and the iterations were terminated when the convergent criterion 10−10 was satisfied, requiring about 12 iterations (see Fig. 1). Table 1 lists the maximum absolute error for each grid configuration used, while numerical results for the field function u at selected positions within the working domain are listed in Table 2. A concrete advantage of the meshless methods is their computational efficiency and, more precisely, their reduced computational cost. The only time-consuming part of the entire computational procedure is the construction of the shape functions and their derivatives. The required computational time for the shape function procedure is of order of seconds on a PC with a dual-core processor of 2.10 GHz and a 3 GB RAM, the precise time frame depending on the grid size, whereas the computational time for the solution of the linearized system is less than 1 s for each iteration of the simulation. More precisely, the CPU time was 46.46 s for a convergence error of 10−6 using the MATLAB built-in function cputime. Fig. 2 shows the contour plots of both the analytical and the numerical solution that, practically, coincide. Example 2. Liouville and Bratu equations

∇2 u ¼

λ2 −u 2 e ; λ ¼ 40 8

ð20Þ

The spatial domain, where the numerical solution is to be calculated, is a unit circle, on the periphery of which the Dirichlet boundary condition u ¼1 is applied. Since the geometry under consideration is irregular, a uniform grid embedded in the spatial domain has been used [8]. A convergence analysis took place starting with a rather coarse grid, namely, a total of 395 nodes, of which 90 nodes along the periphery and the rest in the interior domain. Then, three more grids were tested, as it can be seen in Table 3. For every grid size, the number of boundary nodes remained constant, that is, 90. The support domain for each node included a constant number of 20 nodes. For this example, the initial guess for the iteration procedure was set to zero, and the iterations were terminated when the convergence error dropped below 10−8, requiring about 20 and 40 iterations for the coarse and the dense grid, respectively. The values of the unknown field function for the various nodal configurations used are listed in Table 3. From the results therein we can see that for a grid larger Table 2 Numerical results for the field function u at selected positions using the MPC method and the analytical solution from [17]. x

y

u (MPC)

Exact

0.00 0.25 0.50 0.75 1.00

0.00 0.25 0.50 0.75 1.00

0.0000 0.0156 0.1250 0.4219 1.000

0.0000 0.0156 0.1250 0.4219 1.0000

The Liouville equation is considered next, which is widely used to describe discrete mass systems in the framework of classical mechanics. Again, for the sake of comparison with literature results [10,18], the governing mathematical differential equation

1 0.9 0.8 0.7

Y

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

X Fig. 1. Dependence of relative error on the number of iterations for the non-linear Poisson Eq. (18) on a square domain.

Fig. 2. Contour plot of the analytical (line) and the numerical (symbols) solution of the non-linear Poisson Eq. (18) on a square domain.

Table 1 Absolute deviation from analytical solution [17] for the generalized non-linear Poisson equation, for several grid configurations. Grid Absolute error

41  41 1.336334  10

81  81 −6

3.209042  10

101  101 −7

1.816792  10

141  141 −7

9.965875  10−8

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

Table 3 Numerical results for the Liouville equation in the unit circle, for various grid sizes. Comparison with numerical results from [10].

u at r¼ 0 u at r¼ 0.7 u at r¼ 0.9 u,r at r ¼1

21  21

41  41

81  81

161  161

201  201

Ref. [18]

0.1102 0.5960 0.8593 1.4328

0.1117 0.5960 0.8589 1.4328

0.1124 0.5963 0.8589 1.4329

0.1126 0.5964 0.8589 1.4332

0.1126 0.5964 0.8589 1.4332

0.1127 0.5964 0.8589 1.4332

0.8 0.6 0.4 0.2 0

1

-0.2

0.8

-0.4

0.6

-0.6

0.4

-0.8

0.2

Y

1121

0

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Fig. 4. Contour plot of the numerical solution for the Bratu problem in the unit circle.

-0.2 -0.4 -0.6 -0.8 -1 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

X

results are listed in Table 4. Similar to the previous example, a uniform grid, embedded in the spatial domain has been used. The Dirichlet boundary condition u ¼1 is also employed. The calculations were repeated for a number of grid sizes using the MLS basis functions. Similar conclusions to the aforementioned ones are derived with respect to the convergence and the speed of calculations. The comparison with literature results [10] was also excellent. Fig. 4 shows the contour plot of the unknown field function.

Fig. 3. Contour plot of the numerical solution for the Liouville problem in the unit circle.

Example 3. Poisson–Boltzmann equation

Table 4 Numerical results for the Bratu equation in the unit circle, for various grid sizes. Comparison with numerical results from [10].

The Poisson–Boltzmann (PB) equation is a typical example of non-linear Poisson-type equations. It describes a wide range of physical problems ranging from biochemical and bio-molecular processes to electrostatic interactions between colloidal particles and electro-osmotic flows in micro channels. The governing equation can be written as

u at r¼ 0 u at r¼ 0.7 u at r¼ 0.9 ux at r ¼ 1

21  21

41  41

81  81

161  161

201  201

Ref. [10]

1.3169 1.1553 1.0564 −0.5900

1.3167 1.1552 1.0564 −0.5900

1.3167 1.1552 1.0564 −0.5900

1.3176 1.1559 1.0568 −0.5900

1.3176 1.1559 1.0568 −0.5900

1.3176 1.1559 1.0568 −0.5900

than 121  121 the results converged to the solution, as evidenced in Table 3 upon comparison with numerical results from other methods [10]. The CPU times for convergence error 10−6 were 2.17, 3.58, 11.99, 61.29 and 106.92 s for grid configurations 21  21, 41  41, 81  81, 161  161 and 201  201, respectively. Fig. 3 shows the contour plot of the unknown field function. A similar problem is that of thermal explosion. The governing equation is the Bratu equation ∇2 u ¼ −λeu ; λ ¼

1 e

ð21Þ

The choice of this value for λ was made again for the sake of comparison with literature results [18]. This equation is quite often employed for the determination of regimes of safe operation during combustion and other exothermic processes. The solution domain and boundary condition are the same as in the Liouville example. A convergence analysis has been conducted using coarse and fine grids as in the previous example and the numerical

∇⋅ðε∇uÞ ¼ κ 2 sinhðuÞ þ f

ð22Þ

where u is the unknown field function to be calculated, ε and κ are some known field functions, and f is the source term. The Poisson– Boltzmann equation is a non-linear equation and a linearization procedure is necessary for a numerical solution. To this end, the Newton iteration procedure is used, that is, sinhðunþ1 Þ ¼ sinhðun Þ þ ðunþ1 −un Þcoshðun Þ;

ð23Þ

where un represents the solution at the nth iteration step. Using this expression, the Poisson–Boltzmann equation at the next iteration step (un+1) reads κ2 unþ1 coshðun Þ−∇⋅ðε∇unþ1 Þ ¼ −κ 2 sinhðun Þ þ κ2 un coshðun Þ−f ;

ð24Þ

The iterative procedure is concluded when the absolute value of the difference between un+1 and un is less than a chosen error tolerance. In order to demonstrate the efficiency and the accuracy of the proposed meshless scheme we consider a case where the analytical solution is known [19]. For the problem under consideration, the spatial domain is considered to be a rectangular domain [−1,1]  [−1,1], ε ¼κ ¼1 and the source term f is equal to f¼ 4−sinh (x2+y2+excosy). The exact solution is u ¼x2+y2+excosy and was also used to generate the Dirichlet boundary conditions. Several

1122

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

grids were used, namely 41  41, 81  81, 121  121 and 161  161 and the L2 error norm, defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u u ∑ ðuMPC −uexact Þ2 ui ¼ 1 L2 ¼ u u N t ∑ ðuexact Þ2

1 0.8

ð25Þ

0.6

i¼1

0.4

Table 5 Absolute deviation from analytical solution [19] for the generalized non-linear Poisson equation, for several grid configurations. Grid 41  41 L2

7.270  10

81  81 −10

4.276  10

121  121 −11

8.281  10

161  161 −12

2.595  10

201  201 −12

1.056  10−12

0.2

Y

has been used to study the accuracy and the convergence of the proposed scheme, where uMPC stands for the numerical solution of the present scheme, uexact is the exact solution of the unknown field function u, and N is the total number of the nodes used. The L2 norm errors for the nodal configurations used are listed in Table 5, whereas the plots of the L2 norm error versus the total number of nodes used and the contour lines of the field function u are shown in Fig. 5. Next, we calculate the distribution of electrostatic potential in a static ionic solution or the concentration in diffusion and reaction systems [18,20] using the non-linear Poisson–Boltzmann equation. For instance, in the electrostatic problem, u is the electrostatic potential to be calculated, κ is the inverse Debye–Huckel length, ε is the dielectric coefficient and f is the source term. The solution domain is the rectangle [0,5/8]  [0,5/12], with κ ¼79 and f ¼0 and, the boundary conditions are defined as u ¼8 at x ¼5/8, u ¼8 at ∂u y¼5/12, ∂u ∂x ¼ 0 at x ¼0, and ∂y ¼ 0 at y¼0. A total of 81  81 uniformly distributed nodes have been employed in the calculations. The initial guess for the iteration procedure was set to zero. The entire procedure terminated after 9 iterations when the convergent criterion 10−8 was satisfied. For this calculation, the CPU time was 13.33 s. Fig. 6 shows both the contour and the wireframe plots for the numerical solution of the electrostatic potential over the entire spatial domain. As it can be seen, at the top and right boundaries, the potential drops rapidly to zero. The results compare quite satisfactorily with the ones that were obtained using the Thin Plate Splines (TPS) Radial Basis [18]. Fig. 7 shows the results along the bottom (y¼0) and left (x ¼0) walls. We have to notice the accuracy of the proposed scheme in the case of the steep gradient of the potential near the boundary. There is no unwanted oscillation for the numerical solution, showing the capability of the method to capture the boundary layers. As a final example, we investigate the interaction of circular electrostatic particles, described by the nonlinear Poisson–Boltzmann equation (Eq. (22)), setting the source term f as f¼ sinh(u) and ε ¼κ ¼ 1, as in the first example. The spatial domain is depicted in Fig. 8 and involves two circles with radius r ¼ 0.4, which are centered at (x,y)¼ (−1,0) and (x,y) ¼(1,0). The boundary conditions for the potential u is defined as u ¼1 at the interior boundaries and u ¼2 at the outer boundaries. Since the geometry is irregular, a regular nodal distribution is used that is embedded in the spatial domain [8]. In this way, convergence of the harmonic operator and accuracy of the numerical solution are ensured. The initial guess u(0) ¼ 1 was used and the solution was found to converge to less than 10−10 in 6 iterations. Fig. 9a shows the absolute value of the residuals as a function of the number of iterations, while in Fig. 9b the contour plots of the solution obtained is plotted. Fig. 10 shows the results for the potential function u along the middle line x ¼0, in order to verify the numerical accuracy of the proposed scheme.

0 -0.2 -0.4 -0.6 -0.8 -1 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

X

Fig. 5. (a) Contour plot for the Poisson–Boltzmann problem and (b) L2 norm error versus the total number of nodes used.

The results are compared with the results obtained using the hybrid finite element formulation with F-Trefftz (HFS-FEM) method [19]. 4.2. Transient problems Transient non-linear Poisson-type problems are widely encountered in the modeling of physical phenomena. For example, transient heat conduction and mass diffusion with source terms arise in many different areas of physics and engineering. Representative problems include transient diffusion with chemical reaction in a catalyst pellet and microwave heating process. Efficient solution of non-linear Poisson-type problems is essential to solving numerically a variety of other complex problems, e.g., time-dependent multidimensional viscous and visco-elastic flows. In the present study we examine two representative physical examples of the application of the transient non-linear Poisson equation in order to depict the applicability of the meshless, strong form, collocation method.

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

1123

Fig. 6. (a) Contour and (b) wireframe plots for the Poisson–Boltzmann problem.

4.2.1. Convection–diffusion–reaction As a first case study, we consider the convective diffusion and reaction problem in a rectangular domain of length L¼ 6 and height H ¼0.7 [21]. The governing equation is ∂u ∂u ∂u þ vx þ vy ¼ De ∇2 u−RðuÞ ∂t ∂x ∂y

ð26Þ

where u is the unknown field function and vx and vy are the velocity components that define the movement of the function u. Flow is fully developed in the x-direction (plug flow), while the reaction kinetics is assumed linear, that is, R(u)¼ k1u. The boundary conditions are defined as u¼1 at x¼0 and ∂u=∂n ¼ 0 over the rest of the boundary, while the initial condition is given by u¼0. The numerical solution of the convection–diffusion–reaction problem is obtained using the methodology presented in Section 3, with the differential operator defined as L ¼ vx ð∂u=∂xÞ þ vy ð∂u=∂yÞ−De ∇2 u þ k1 u. We consider the case where the flow is in x-direction (vy ¼0). Emphasis is placed here on the accuracy of the method at high Peclet numbers. A uniform grid configuration is used and, a total amount of 6989 nodes has been used for the calculation of the concentration distribution. The CPU time was 16.25 s for convergence error 10−8. Fig. 11a and 11b show the results for vx ¼1 m/s and vx ¼6 m/s, and for k1 ¼0.278 s−1. The agreement with the operator splitting method is excellent [21].

Fig. 7. Profile of solution along (a) the bottom (y¼ 0) and (b) the left (x¼ 0) wall for the Poisson–Boltzmann problem using the MPC (lines) method.

y

1

2

1

x 2

4

4

Fig. 8. Geometry configuration for the electrostatic Poisson–Boltzmann problem around a pair of particles.

the forced heat equation 4.2.2. Microwave heating The problem of heating a square slab using microwave radiation energy is considered next. The specific problem is challenging from both the physical and the numerical point of view, since a thermal boundary layer is created. The process can be modeled by

 ∂T ¼ ∇2 T þ ηðTÞEj2 ∂t

ð27Þ

describing the absorption and conduction of heat, along with Maxwell equations that govern the propagation of microwave radiation through the material. In the above equation, η(T) is the

1124

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

2

Y

1

0

-1

-2

-4

-3

-2

-1

0

1

2

3

4

X Fig. 9. (a) Error convergence and (b) contour plots for the electrostatic problem around a pair of particles.

Fig. 11. Comparison of MPC (lines) and operator splitting (squares) solution of the convective diffusion and reaction problem for (a) ux ¼ 1 m/s and (b) ux ¼6 m/s for k1 ¼ 0.278 s−l.

where γ is the decay constant. Finally, the governing equation is written as ∂T ¼ ∇2 T þ βe−γx T n ∂t

ð29Þ

Typical values of n¼ 2 or 3 are used in practical modeling work. The non-linear equation was linearized based on the linearization scheme proposed in [22]. The aforementioned equation can then be written as [22]

Fig. 10. Variation of potential along the middle line x ¼0 for the electrostatic problem. Comparison with the hybrid finite element method.

thermal absorptivity, |E| is the amplitude of the electric field, assuming a constant thermal diffusivity normalized to unity, rendering t a scaled time. In the present study the material properties are set constant and the amplitude of the electric field is exponentially dependent on the spatial variables, that is, jEj ¼ e−ðγx=2Þ

ð28Þ

∂T ¼ ∇2 T þ βe−γx ð1−nÞT nk þ nβe−γx T n−1 k T f or t∈ðkδt; ðk þ 1ÞδtÞ ∂t

ð30Þ

where Tk is the solution at time t¼ kdt. The error tolerance was set to 10−6 and the initial condition was used as first iteration. Numerical results at selected positions within the working domain for the cases of n ¼2 and n ¼3 are listed in Table 6. This calculation required a CPU time of 20.24 s. The corresponding values were found to be in very good agreement with the numerical results that were obtained by different numerical methods, namely Dual Reciprocity Method (DRM) and Laplace Transform Dual Reciprocity Method (LTDRM) [22].

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

1125

Table 6 Comparison of numerical results of different methods for the forced heat equation at specific positions [22]. x

0.1428 0.2857 0.4285 0.7142

y

0.8571 0.7142 0.5714 0.4285

DRM

LTDRM

MPC

a

b

c

a

b

c

a

b

c

1.139 1.404 1.612 1.496

1.195 1.485 1.613 1.354

1.240 1.510 1.550 1.252

1.137 1.395 1.597 1.485

1.192 1.486 1.598 1.346

1.237 1.502 1.541 1.248

1.138 1.402 1.609 1.495

1.194 1.486 1.611 1.352

1.239 1.508 1.549 1.251

4.2.3. Reaction–diffusion Brusselator system Finally, we deal with the numerical solution of the coupled pair of nonlinear partial differential equations:  2  ∂u ∂ u ∂2 u ¼ D1 þ þ a1 u þ Aðu; vÞ þ f 1 ðx; y; tÞ ∂t ∂x2 ∂y2  2  ∂v ∂ v ∂2 v ¼ D2 þ ð31Þ þ a2 v þ βu þ Bðu; vÞ þ f 2 ðx; y; tÞ ∂t ∂x2 ∂y2 with given initial and Dirichlet and/or Neumann's boundary conditions in the two-dimensional region Ω, where D1, D2, α1, α2 and β are given constants, A and B are functions of the field variables u and v, and f1 and f2 are assumed to be prescribed sources. In our case, for the two-component reaction system, u(x,t) and v(x,t) stand for concentrations and D1, D2 for the diffusion coefficients of the chemical species. More specifically, the non-linear reaction–diffusion Brusselator system [23,24] was considered here. This particular system is used to describe chemical reaction and diffusion with non-linear oscillations [23] and references there in. The two-dimensional reaction–diffusion Brusselator system under consideration is  2  ∂u ∂ u ∂2 u ¼ D1 þ 2 −ðA þ 1Þu þ u2 v þ B 2 ∂t ∂x ∂y  2  ∂v ∂ v ∂2 v ¼ D2 −u2 v þ Au þ ð32Þ ∂t ∂x2 ∂y2 defined in the two-dimensional region Ω ¼[0,1]2, with D1 ¼ D2 ¼0.002, A ¼0.5 and B ¼1. The initial conditions are defined as 1 2 1 3 x − x 2 3 1 2 1 3 vðx; y; 0Þ ¼ y − y 2 3

uðx; y; 0Þ ¼

ð33Þ

while Neumann boundary conditions were considered: ∂uðx; y; tÞ  ∂uðx; y; tÞ  ∂uðx; y; tÞ  ∂uðx; y; tÞ  ¼ ¼ ¼ ¼0     ∂x ∂x ∂y ∂y x¼0 x¼1 y¼0 y¼1 ∂vðx; y; tÞ  ∂vðx; y; tÞ  ∂vðx; y; tÞ  ∂vðx; y; tÞ  ¼ ¼ ¼ ¼0     ∂x ∂x ∂y ∂y x¼0 x¼1 y¼0 y¼1 ð34Þ The non-linear problem is not known to have an analytical solution. It is stated that for small values of the diffusion coefficients D1, D2, and if 1−A+B2 40 then the steady state solution of the Brusselator system converges to equilibrium point (u,v)¼(B,A/ B) [24]. For the numerical solution of the problem under investigation, we initially discretize the time derivatives of the system using the first order forward difference formula and then apply the θweighted (0≤θ≤1) scheme to the space derivatives at two successive time levels k and k+1, that is, h i uðkþ1Þ −uðkÞ ¼ B þ θ ðu2 vÞðkþ1Þ −ðA þ 1Þuðkþ1Þ þ D1 ð∇2 uÞðkþ1Þ þ δt

h i ð1−θÞ ðu2 vÞðkÞ −ðA þ 1ÞuðkÞ þ D1 ð∇2 uÞðkÞ h i vðkþ1Þ −vðkÞ ¼ θ Auðkþ1Þ −ðu2 vÞðkþ1Þ þ D2 ð∇2 vÞðkþ1Þ δt þð1−θÞ½AuðkÞ −ðu2 vÞðkÞ þ D2 ð∇2 vÞðkÞ  (k)

(k)

(k)

(k)

(k+1)

ð35Þ

(k)

¼ t +δt, with δt being where u ¼u(x,y,t ), v ¼v(x,y,t ), t the time step. Eq. (35) can be linearized by approximating the nonlinear term (u2v) (k+1) with Taylors series formula (see also [24]) ðu2 vÞðkþ1Þ ¼ ðu2 ÞðkÞ vðkþ1Þ þ 2uðkÞ vðkÞ uðkþ1Þ −2ðu2 ÞðkÞ vðkÞ

ð36Þ

Eq. (35) using Eq. (36) become [24] −θδt½ðu2 ÞðkÞ vðkþ1Þ þ 2uðkÞ vðkÞ uðkþ1Þ −ðA þ 1Þuðkþ1Þ þ D1 ð∇2 uÞðkþ1Þ  n o ¼ uðkÞ þ δt½B þ ð1−3θÞðu2 ÞðkÞ vðkÞ þ ð1−θÞ −ðA þ 1ÞuðkÞ þ D1 ð∇2 uÞðkÞ 

ðkþ1Þ

u

vðkþ1Þ −θδt½Auðkþ1Þ −ðu2 ÞðkÞ vðkþ1Þ −2uðkÞ vðkÞ uðkþ1Þ þ D2 ð∇2 vÞðkþ1Þ  n o ¼ vðkÞ þ δt½ð3θ−1Þðu2 ÞðkÞ vðkÞ þ ð1−θÞ AuðkÞ þ D2 ð∇2 vÞðkÞ 

ð37Þ

We can define Φ, Φ,s, Φ,ss, with s ¼x,y, in the context of MLS, as the matrices that give the unknown field function approximation values and their spatial derivatives up to second order. These matrices can be split into two matrices Φd and Φb corresponding to Nd interior points and Nb boundary points " # Φd (N ¼ Nd +Nb) in the following form Φ ¼ ∈ℝNN , with N the Φb total number of nodes. The final solution w(k) can be written, using the field vectors u(k) ¼ Φu and v(k) ¼ Φv, as " # uðkÞ ð38Þ wðkÞ ¼ Μ ðkÞ ∈ℝ2N1 v  Φ 0 with Μ ¼ ∈ℝ2N2N . Using Eq. (38) the system in Eq. (37) 0 Φ can be written in matrix notation as Ηðkþ1Þ wðkþ1Þ ¼ ΗðkÞ wðkÞ þ R ðkþ1Þ " ðkþ1Þ # Ηu where Ηðkþ1Þ ¼ ∈ℝ2N2N , with Ηðvkþ1Þ

ð39Þ

Ηðkþ1Þ u 2

 Φd −2θδtðuðkÞ ⋅uðkÞ Þ∘Φd þ ðA þ 1ÞθδtΦd −θδtD1 ðΦd;xx þ Φd;yy Þ ¼4

 Gu Φb 2

 −θδtAΦd þ 2θδtðuðkÞ ⋅vðkÞ Þ∘Φd

Ηðkþ1Þ ¼4 v

f0g

" and ΗðkÞ ¼ 2

4 ΗðkÞ u ¼ 2

4 ΗðkÞ v ¼

n

n

−θδtðuðkÞ Þ2 ∘Φd f0g

o3 5

o3 Φd þ θδtðuðkÞ Þ2 ∘Φd −θδtD1 ðΦd;xx þ Φd;yy Þ 5

 G v Φb

# ðkÞ

Ηu

∈ℝ2N2N ,

ΗðvkÞ

Φd −ð1−θÞðΑ þ 1ÞδtΦd þ ð1−θÞδtD1 ðΦd;xx þ Φd;yy Þ f 0g ð1−θÞδtAΦd f 0g



n



n o3 ð1−3θÞδtðuðkÞ Þ2 ∘Φd 5; f 0g

Φd þ ð3θ−1ÞδtðuðkÞ Þ2 ∘Φd þ ð1−θÞδtD1 ðΦd;xx þ Φd;yy Þ f 0g

o3 5

1126

G.C. Bourantas, V.N. Burganos / Engineering Analysis with Boundary Elements 37 (2013) 1117–1126

domain. Additionally, the refinement procedure is also straightforward and easy to implement and, in combination with the linearization method, renders the proposed meshless method quite suitable for this type of problems. Finally, it must be noted that the computational cost of the method is sufficiently low to make the method attractive in a number of practical applications.

Table 7 Meshless point collocation solution at selected nodes (x,y). t

1 2 3 4 5 6 7 8 9 10 ↓ ∞

(0.2, 0.2)

(0.8, 0.9)

v

u

v

u

v

u

v

0.5329 0.7043 0.8192 0.9117 0.9727 1.0001 1.0063 1.0046 1.0020 1.0004 ↓ 1

0.1709 0.3730 0.4956 0.5363 0.5304 0.5142 0.5036 0.4995 0.4988 0.4992 ↓ 0.5

0.5515 0.7289 0.8470 0.9346 0.9855 1.0044 1.0066 1.0039 1.0014 1.0002 ↓ 1

0.2575 0.4306 0.5217 0.5402 0.5261 0.5101 0.5016 0.4989 0.4988 0.4993 ↓ 0.5

0.5558 0.7268 0.8428 0.9308 0.9832 1.0036 1.0065 1.0040 1.0015 1.0002 ↓ 1

0.2410 0.4201 0.5170 0.5394 0.5268 0.5109 0.5019 0.4990 0.4988 0.4993 ↓ 0.5

0.5819 0.7537 0.8702 0.9512 0.9934 1.0066 1.0064 1.0033 1.0010 1.0000 ↓ 1

0.3180 0.4683 0.5356 0.5400 0.5221 0.5073 0.5004 0.4987 0.4989 0.4994 ↓ 0.5

"

#

R ðkþ1Þ ¼ Ε

(0.5, 0.5)

u

and

ðkþ1Þ

(0.4, 0.2)

¼

Εðkþ1Þ

J ðkþ1Þ " # Bδt hi

∈ℝ2N1 ; with "

and J

ðkþ1Þ

¼

0

#

qi

Gu and Gv are the matrices of order nb  N defining the boundary conditions, Dirichlet (G≡Φ) or Neumann (ndΦ,j, n is the unit vector normal to the boundary and j¼ x,y) and ht and qt are scalar functions defining the boundary conditions. From Eq. (39) we can write for the solution wðkþ1Þ ¼ ðΗðkþ1Þ Þ−1 ðΗðkÞ wðkÞ þ R ðkþ1Þ Þ

ð40Þ

For comparison reasons, the numerical results of the above method were compared against the results of the meshless local integral equation [23]. For that reason a uniform grid of 21  21 nodes was used and the time step δt was set to δt ¼0.1, requiring a CPU time of 12.48 s.The results are listed in Table 7 for the meshless point collocation method. The agreement with the meshless local integral equation method is excellent. 5. Conclusions The meshless point collocation method has been successfully used for the numerical solution of non-linear transient and steady state Poisson-type equations. A robust and easy to implement numerical method is applied, yielding accurate numerical results. The great number of applications that can be modeled with the help of the non-linear Poisson equation renders the present method very useful in a variety of problems. In fact, the traditional numerical methods, such as FDM, used for solving this type of problems are computationally less efficient compared with the meshless methods. The advantages of the meshless methods rest mainly with the interpolation procedure. A set of nodes that are randomly or regularly distributed over the spatial domain can be used in order to compute the values of an unknown field function in a straightforward manner, even in a fully irregular spatial

References [1] Li Z, Pao CV, Qiao ZA. Finite difference method and analysis for 2D nonlinear Poisson–Boltzmann equations. J Sci Comput 2007;30:61–81. [2] Lin YP, Zhang T. Finite element methods for nonlinear Sobolev equations with nonlinear boundary conditions. J Math Anal Appl 1992;165:180–91. [3] Kasab JJ, Karur SR, Ramachandran PA. Quasilinear boundary element method for nonlinear Poisson type problems. Eng Anal Boundary Elem 1995;15:277–82. [4] Atluri SN, Shen SP. The meshless local Petrov–Galerkin (MLPG) method. Encino USA: Tech Science Press; 2002 440 p. [5] Liu GR. Mesh free methods, moving beyond the finite element method. Boca Raton, FL: CRC Press; 2002. [6] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. Int J Numer Meth Eng 1994;37:229–56. [7] Fasshauer. G. Meshfree approximation methods with MATLAB. Singapore: World Scientific; 2007. [8] Bourantas GC, Skouras ED, Nikiforidis GC. Adaptive support domain implementation on the moving least squares approximation for Mfree methods applied on elliptic and parabolic PDE problems using strong-form description. Comput Model Eng Sci 2009;43:1–25. [9] Bourantas GC, Korfiatis DP, Loukopoulos VC, Thoma K-ATh. Numerical simulation of the unsteady non-linear heat transfer problems. Application on nanosecond laser annealing of Si. Appl Surf Sci 2012;12:21–5. [10] Balakrishnan K, Ramachandran PA. Osculatory interpolation in the method of fundamental solution for nonlinear Poisson problems. J Comput Phys 2001;172:1–18. [11] Chen CS. The method of fundamental solutions for non-linear thermal explosions. Commun Numer Methods Eng 1995;11:675–81. [12] Zhu T, Zhang J, Atluri SN. A meshless local boundary integral equation (LBIE) method for solving nonlinear problems. Comput Mech 1998;22:174–86. [13] Balakrishnan K, Sureshkumar R, Ramachandran PA. An operator splittingradial basis function method for the solution of transient nonlinear Poisson problems. Comput Math Appl 2002;43:289–304. [14] Lancaster P, Salkauskas K. Surfaces generated by moving least squares method. Math Comput 1981;37:141–58. [15] Aluru NR. A point collocation method based on reproducing kernel approximations. Int J Num Methods Eng 2000;47:1083–121. [16] Kim DW, Liu WK. Maximum principle and convergence analysis for the meshfree point collocation method. SIAM J Numer Anal 2006;44:515–39. [17] Wang H, Qin QH. A meshless method for generalized linear or nonlinear Poisson-type problems. Eng Anal Boundary Elem 2006;30:515–21. [18] Liu X, Liu GR, Tai K, Lam KY. Radial point interpolation collocation method (RPICM) for the solution of nonlinear Poisson problems. Comput Mech 2005;36:298–306. [19] Wang H, Qin Q-H, Liang X-P. Solving the nonlinear Poisson-type problems with F-Trefftz hybrid finite element model. Eng Anal Boundary Elem 2012;36:39–46. [20] Yang C, Li D, Maliyah JH. Modeling forced liquid convection in rectangular microchannels with electrokinetic effects. Int J Heat Mass Transfer 1998;41:4229–49. [21] Balakrishnan K, Sureshkumar R, Ramachandran PA. An operator splittingradial basis function method for the solution of transient nonlinear Poisson problems. Comput Math Appl 2002;43:289–304. [22] Zhu S, Satravaha P. An efficient computational method for modeling transient heat conduction with nonlinear source terms. Appl Math Modeling 1996;20:513–22. [23] Shirzadi A, Sladek V, Sladek J. A local integral equation formulation to solve coupled nonlinear reaction–diffusion equations by using moving least square approximation. Eng Anal Boundary Elem 2013;37:8–14. [24] Islam S, Ali A, Haq S. A computational modeling of the behavior of the twodimensional reaction–diffusion Brusselator system. Appl Math Modelling 2010;34:3896–909.