THE USE OF THE MESH FREE METHODS (RADIAL BASIS FUNCTIONS) IN THE MODELING OF RADIONUCLIDE MIGRATION AND MOVING BOUNDARY VALUE PROBLEMS
LEOPOLD VRANKAR, FRANC RUNOVC and GORAN TURK
About the authors
Leopold Vrankar
Slovenian Nuclear Safety Administration Železna cesta 16, 1001 Ljubljana, Slovenia E-mail: leopold.vrankar@uni-mb.si
Franc Runovc
University of Ljubljana,
Faculty of Natural Sciences and Engineering
Aškerčeva cesta 12, 1000 Ljubljana, Slovenia
E-mail: franc.runovc@ntf.uni-lj.si
Goran Turk
University of Ljubljana,
Faculty of Civil and Geodetic Engineering
Jamova cesta 2, 1000 Ljubljana, Slovenia
E-mail: gturk@fgg.uni-lj.si
Abstract
Recently, the mesh free methods (radial basis functions-RBFs) have emerged as a novel computing method in the scientific and engineering computing community. The numerical solution of partial differential equations (PDEs) has been usually obtained by finite difference methods (FDM), finite element methods (FEM) and boundary elements methods (BEM). These conventional numerical methods still have some drawbacks. For example, the construction of the mesh in two or more dimensions is a nontrivial problem. Solving PDEs using radial basis function (RBF) collocations is an attractive alternative to these traditional methods because no tedious mesh generation is required. We compare the mesh free method, which uses radial basis functions, with the traditional finite difference scheme and analytical solutions. We will present some examples of using RBFs in geosta-tistical analysis of radionuclide migration modeling. The advection-dispersion equation will be used in the Eulerian and Lagrangian forms. Stefan's or moving boundary value problems will also be presented. The position of the moving boundary will be simulated by the moving data centers method and level set method.
Keywords
mesh free methods, radial basis functions, finite difference methods, finite element methods, boundary elements methods, geostatistics, Eulerian method, Lagrangian method, level set method
1 introduction
In recent years, the mesh free methods have emerged as a novel computing method in the scientific and engineering computing community. Traditionally, the most popular methods have been the finite element methods (FEM), the finite difference methods (FDM), and the boundary element method (BEM). In spite of their great success in solving scientific and engineering problems over the past four decades, these conventional numerical methods still have some drawbacks that impair their computational efficiency and even limit their applicability to more practical problems, particularly in a three-dimensional space.
The term mesh free method refers to the ability of the method to solve the given differential equations from a set of unstructured nodes, i.e. without any pre-defined connection or relationship among the nodes. Instead of generating a mesh, mesh free methods use scattered nodes, which can be randomly distributed or patterned through the computational domain.
During the past decade, increasing attention has been given to the development of mesh free methods using RBFs for the numerical solution of PDEs. There are two major developments in this direction. The first is the method of fundamental solutions (MFS) coupled with the dual reciprocity method (DRM), which evolved from the dual reciprocity boundary element method (DRBEM). More details about MFS can be found in review papers [1] and [2]. RBFs played a key role in the theoretical establishment and applications in the development of the DRM. With the combined features of the MFS and DRM, a mesh free numerical scheme for solving PDEs has been achieved.
ACTA GeOTeCHNICfl SLOVENICA, 2007/i 43.
L. VRflNKflR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONSl IN THE MODELING OF RflDIONUCLIDe MIGRATION AND
The second mesh free method using RBFs is the Kansa method [3-4], where the RBFs are directly implemented for the approximation of the solution of PDEs introducing the concept of solving PDEs, using radial basic functions for hyperbolic, parabolic and elliptic PDEs. A key feature of the RBF method is that it does not need a grid. In contrast to the MFS-DRM boundary method, the Kansa method is considered to be a domain type method, which has many features similar to the finite element method.
In our case, we use the Kansa method which uses radial basis functions. We compare the results with a traditional finite difference scheme and analytical solutions.
We will present two examples: some previous work results of using RBFs in geostatistical analysis of transport modeling of the radionuclide migration and moving boundary value problems.
2 radial basis functions method
The base of this approach is its employment of highorder interpolating functions to approximate solutions of differential equations. All RBFs possess the property that their values are only determined by distance and have nothing to do with directions. Kansa [4] introduced multiquadric functions to solve hyperbolic, parabolic and elliptic differential equations with collocation methods. This method is an asymmetric collocation set-up in which boundary conditions are treated separately from the interior problem. One of the most powerful RBF methods is based on multiquadric basis functions (MQ), first used by R. L. Hardy [5]. It is important to mention that the MQ was until now efficiently used in transport modeling [6].
A radial basis function is a function <(x) = <(||x — x. ||), which depends only on the distance between x e Rd and a fixed point x. e Rd. Here, q> is continuous and bounded on any bounded sub-domain QCRd. Let r denote the Euclidean distance between any pair of points in the domain O. The commonly used radial basis functions are: linear (<(r) = r), cubic (<(r) = r3), thin-plate spline (<(r) = r2 log r), Gaussian (<(r) = e~ar) and multiquadric (MQ) (<(r) = (r2 + c2))3). Commonly used values for ft are -1/2 and 1/2. The parameter c > 0 is a shape parameter controlling the fitting of a smoothing surface to the data. We usually use MQ or inverse MQ RBFs.
To introduce RBF collocation methods, we consider a PDE in the form of
Lu(x) = f (x), in He Rd, (1)
Bu(x) = g (x), on dH, (2)
where u is the concentration, d denotes the dimension, 3O is the boundary of the domain O, L is the differential operator on the interior, and B is the operator that specifies the boundary conditions of the Dirichlet, Neumann or mixed type.
Using Kansa's asymmetric multiquadric collocation method, the unknown PDE solution u is approximated by RBFs in the form:
N	M
u » U(x) = £a< (x) + £ Ylp1 (x), (3)
j=i	1=1
where q> can be any of above mentioned radial basis function, pi,...,pM e ndm , is a polynomial of degree m or less. Let (x.)N=1 be the N collocation points in H U dH . We assume the collocation points are arranged in such a way that the first N, points are in O, whereas the last Nb points are on 3O. To evaluate N+M unknown coefficients, N+M linearly independent equations are needed. Ensuring that U(x) satisfies (1) and (2) at the collocation points results in a good approximation of the solution u. The first N equations are given by
N	M
£a.L<j (x,) + £ylLpi (x) = f (x,), for i = 1,...,N,, (4)
j=1	1=1
N	M
£a.B< (xi) + £ YlBpi (x) = g(xi), for i = N, +1,...,N, + Nb .(5)
j=1 1-1
The last M equations could be obtained by imposing some extra condition on v(.):
£ ap (x.) = 0, k = 1,..., M. (6)
j=1
The choice of a basis function is another flexible feature of RBF methods. RBFs can be globally supported, infinitely differentiable, and contain a free parameter, c. This leads to a full coefficient matrix or a dense interpolation matrix. The shape parameter c affects both the accuracy of the approximation and the conditioning of the interpolation matrix. The optimal shape parameter c is still an open question. In our case we used an iterative mode by monitoring the spatial distribution of the residual errors
2. ACTA GEOTECHNICA SLOVENICfl, 2007/l
L. VRANKAR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONS) IN THE MODELING OF RADIONUCLIDE MIGRATION AND
in O and 30 as a function of c. The iterations are terminated when errors are smaller then a specified bound. This map is then used to guide the search of the optimal shape parameter c that gives the best approximation of the solution.
3 the modeling of
radionuclide migration
Assessment of the release and the transport of long-lived radioactive nuclides from a repository to the biological environment is an important part of the safety analysis of repository concepts. In this assessment, mathematical models describing the mechanisms involved in the nuclide transport from the repository to the biosphere are essential tools. For example, the groundwater models are mathematical representations of the flow of water and the transport of solutes in the subsurface. Models are used to compute the hydraulic head, velocity, concentration, etc., from hydrologic and mass inputs, hydro geologic and mass-transfer parameters, and conditions at the boundary of the domain.
3-1 GEOSTATISTICS
Many processes are inherently uncertain, and this uncertainty is handled through the use of stochastic realizations. The goal of stochastic simulation is to reproduce geological texture in a set of equiprobable simulated realizations. In mathematical terms, the most convenient method for simulation is sequential Gaussian simulation, because all successive conditional distributions from which simulated values are drawn are Gaussian with parameters determined by the solution of a simple kriging system.
3.1.1 simulation of hydraulic conductivity
The hydraulic conductivity was generated at different points based on different input data. The hydraulic conductivity at 8 different points is given (values are: 66.00, 71.00, 73.00, 75.00, 76.52, 77.02, 79.74, 83.41 [m/y].
The distribution of hydraulic conductivity for one specific simulation is shown in Fig. 1 (the numbers on the axes in Fig. 1 are distances in metres).
The coordinates of these values are also presented in Fig. 1 and marked with +. The following variogram parameters are chosen: positive variance contribution
or sill is equal to 0.8 and the nugget effect is 0.2. Simple kriging is chosen as the type of kriging. A spherical model is chosen as a type of a variogram structure.
Figure 1. Distribution of hydraulic conductivity.
3.2 SIMULATION OF DARCY VeLOCITY
For the calculation of velocity in principal directions, Darcy's law is used. The velocities were determined from the pressure of the fluid p by solving the Laplace differential equation. Homogeneous and anisotropic porous media and incompressible fluid are assumed. The equation has the following form:
+ 0, c
where Kx and Ky are the components of a hydraulic conductivity tensor. The corresponding Neumann and Dirichlet boundary conditions, defined along the boundary, are presented in Fig. 2.
n-VpM)
Figure 2. Boundary conditions.
For the calculation of velocity in principal directions we use Darcy's law:
K^ dp,
upa dx
Ky dp
v = — ,
y upa dy
(8)
ACTA GEOTECHNICA SLOVENICA, 2007/i 45.
L. VRflNKflR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONSl IN THE MODELING OF RflDIONUCLIDe MIGRATION AND
where p is the density of the fluid and a is the gravitational acceleration. The Laplace equation was solved by using direct collocation [7]. Fig. 3 presents the velocity vector. In our case we considered anisotropic porous media, thus we obtained the velocity vectors of different size which depends on porosity and permeability.
Figure 3. Calculated Darcy's velocity.
3.3 SIMULATION OF THE CONTAMINANT CONCENTRATIONS WITH ADVEC-TION-DISPERSION EQUATION (AD)
For the parabolic problem, we consider the implicit scheme [8] and radial basis functions formulation is presented in [9].
The simulation was implemented for a rectangular area which was 600 m long and 300 m wide. The source (initial condition) was Thorium (Th-230) with the activity of 1 MBq. The location of the radioactive source is presented in Fig. 3 with the symbol ♦.
The longitudinal dispersivity ax is 200 m and the transversal dispersivity ay is 20 m. For the porosity w we used values between 0.25 and 0.26. The retardation constant R is 800. The average contaminant concentrations distribution for the radial basis function method is shown in Fig. 4 (the axes in Fig. 4 present distances in metres).
3.3.1 The Eulerian form of the ad
Figure 4. Distribution of average of contaminant concentrations (Radial basis function).
The velocities obtained from Laplace equation are used in the advection-dispersion equation. The AD equation for transport through the saturated porous media zone at a macroscopic level with retardation and decay is
„du Dx d2u Dy d2u du R— = —^—- +--^—- - vx--R u,
dt w dx w dy dx
(x, y) en, 0 < t < T, = g (x, y, t )
(9)
l(x, y )eaa
u |t=0 = h(x,y),
0 < t < T, (x,y) en,
where x is the Eulerian groundwater flow axis and y is the Eulerian transverse axis in the 2D problem, u is the concentration of contaminant in the groundwater [Bqm3], Dx and Dy are the components of dispersion tensor [m2/y] in saturated zone, w is porosity of the saturated zone [-], R is the retardation factor in the saturated zone [-] and X=ln2/t0 is the radioactive decay constant [1/y]. t0 is a half-life. In these cases "y" means years.
The traditional finite difference scheme was used for solving the Laplace and advection-dispersion equation. For the approximation of the second derivatives we used the central difference with respect to x and y. The average contaminant concentrations distribution for the finite difference method is shown in Fig. 5 (the axes in Fig. 5 present distances in metres).
Figure 5. Distribution of average of contaminant concentrations (FDM).
46. ACTA GEOTECHNICA SLOVENICfl, 2007/l
L. VRANKAR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONS) IN THE MODELING OF RADIONUCLIDE MIGRATION AND
3.3.2 The Lagrangian form of the ad
In this case the time-derivative term and the advection term of Eq. (9) are expressed as a material derivative:
du du „ -=--h v • Vu.
dt dt
(10)
After including the material derivative into the advec-tion-dispersion equation (9), we have:
du dt
+X u =
d^ du+D. du
Ru dx2 Ru dy2
(11)
The solution of non-homogeneous ordinary differential equations (ODEs) can be found as a superposition of homogeneous and particular solutions. The particular solution of ODEs was found by the method of constant modification. In our case it was assumed that X is constant in each time step [9].
The calculation of radioactive concentrations in partly heterogeneous porous media was also carried out. The results of radioactive concentrations for one particular simulation in partly heterogeneous porous media for the Eulerian and Lagrangian method are presented in Figs. 6 and 7. The partly heterogeneous porous media means in our case the area which was split into three regions. In the central region we prescribed the conductivity of porous media with the value of 320 [m/y], and elsewhere the prescribed conductivity was 50 [m/y].
Figure 6. Concentrations and conductivity in partly heterogeneous porous media size 50 and 320 [m/y] (Eulerian method).
Figure 7. Concentrations and conductivity in partly heterogeneous porous media size 50 and 320 [m/y] (Lagrangian method).
4 moving boundary value problems
4.1 PROBLEM DEFINITION AND MOVING DATA CENTERS METHOD
Many physical processes involve heat conduction and materials undergoing a change of phase. Examples include the safety studies of nuclear reactors (the molten corium concrete interaction), casting of metals, geophysics and industrial applications involving metals, oil, and plastics. The molten core discharged to the containment cavity will interact with the concrete basement if it is not, or cannot be, cooled below the solid us temperature of the concrete. The molten core concrete interaction results in decomposition and melting of very large quantities of carbon dioxide and steam.
Several numerical methods have been developed to solve various Stefan's problems. Crank [10] provides a good introduction to the Stefan's problems and presents an elaborate collection of numerical methods for these problems. According to Crank the numerical methods for moving boundary problems can be classified in three categories: front-tracking methods, front-capturing methods and hybrid methods. We follow front-tracking methods (moving data centers method) which use an explicit representation of the interface, given by a set of points lying on the interface location, which must be updated at each time step.
Heat treatment of metals is often used to optimize mechanical properties. During heat treatment, the metallurgical state of the alloy changes. This change can involve the phase present at a given location or the morphology of the various phases. One of these processes, which is both of large industrial and scientific interest and amenable to modeling, is the dissolution of the second-phase particles in a matrix with a uniform initial composition. The position of the moving boundary will be simulated by moving data centers method.
4.1.1 The physical model
We consider the solid state phase transformation problem in binary metallic alloys which is described in [11]. In this problem a volume of constant composition is surrounded by a diffusive phase. In the interface between the particle and the diffusive phase a constant concentration is assumed, and the gradient of the concentration causes the movement of the interface.
ACTA GEOTECHNICA SLOVENICA, 2007/i 47.
L. VRflNKflR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONSl IN THE MODELING OF RflDIONUCLIDe MIGRATION AND
4.1.2 The mathematical model
4.2 THE LEVEL SET METHOD
We studied the domain O containing a diffusive phase Oäp and the part where the material characteristic Opart remains of constant composition cpart . The particle dissolves due to Fickian diffusion in the diffusive phase. The concentration at the interface r, separating Opart and Oäp , is assumed to be given by the constant value co.. The concentration gradient on the side of Oäp at r causes its displacement. The governing equations and boundary conditions of this problem are:
du
—(x,t) = DAu(x,t), x eüäp(t),t >0, (12)
dt
u(x, t ) = upart, x enpart (t),t > 0, (13)
u(x,t) = uso1, x er(t),t > 0,
. du
(14)
(upart - usd)vn(x,t) = D—(x,t), x er(t),t > 0, (15) dn
The level set method has gained much popularity for solving moving boundary problems. It was firstly introduced by Osher and Sethian [13]. The level set function captures the interface position as its zero level set, and it is advected by introduction of a hyperbolic equation into the governing set of equations.
4.2.1 The level set formulation
In the level set formulation of the moving interface, the interfaces, denoted by r, are represented implicitly through a level set function 0(x, t), where x is a position of the interface, t is a point of time. Usually, the 0 is defined as a signed distance function to the interface. The moving interface is then captured in all times by locating the set of r(t) for which 0 vanishes. The level set function is advected with time by a transport equation which is known as the level set equation:
^ + Vn IWI = 0, dt
¿(x,0) = ¿.(x) (17)
du dn
(x,t) = 0, xedüäp(t)\r(t),t>0, (16)
where x is the coordinate vector of a point in O, D means the diffusivity constant, n is the unit normal vector on the interface pointing outward with respect to Opart(t) and vn is the normal component of the velocity of the interface. The initial concentration u(x, 0) inside the diffusive phase is given.
4.1.3 The Numerical solution Method
Our interest is to give an accurate discretization of the moving boundary conditions. Here we present an inter-polative moving data centers method by which the data centers are computed for each time step and the solution is interpolated from the old data centers to the new ones. The equations are solved with the collocation method using RBFs. An outline of the algorithm is:
I.	Compute the concentrations profiles solving Eqs. (12-14) and (16).
II.	Predict the position of the boundary s1 at the new time-step: s1 (t + At) using Eq. (15).
III.	Once the boundary is moved, the concentration u
can be computed in the new region using Eq. (12). The solution is interpolated from the old point loca-
tion to the new one.
where 0o(x) embeds the initial position of the interface and vn is the normal component of the velocity of the interface:
v„ = v
'IWI '
(18)
where	is the unit normal to the surface N.
If we take into account a continuous extension of the interface velocity v, the evolution of the level set function can be done by the hyperbolic equation for the level set equation:
dt
+ v = 0, ¿(x,0) = ¿0(x) (19)
In our case the continuous extension of the velocity v is taken as the (steady) solution of the following evolution equation [14]:
ËVL ±ËVL. n = 0,
dr dx
(20)
where r denotes a fictitious time step not related to the main time step and the sign is determined from the normal direction of the level set function.
The RBFs are incorporated into level set methods to construct a more efficient approach. At the initial time,
48. ACTA GEOTECHNICA SLOVENICfl, 2007/l
L. VRANKAR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONS) IN THE MODELING OF RADIONUCLIDE MIGRATION AND
all the time dependent variables should be specified over entire domain. The initial value problem (17) can be considered equivalent to an interpolation problem, and hence the starting point of the use of RBFs to solve partial differential equations is the interpolation problem. Further, the spatial portion is approximated by the RBFs and the temporal variations are approximated by the time dependent expansion coefficients.
4.2.2 rbf implicit modeling of the Level set function
Interpolation of the level set function
In the present implicit modeling, the MQ RBFs is used to interpolate the scalar implicit level set functions 0(x) with N points by using N MQs centered at these points. The resulting RBF interpolant of the implicit function can be written as
N	M
) = £ a^j (x) + ^ 7l p} (x), (21)
j=i	1=1
Because of the introduction of this polynomial, the RBF interpolant of 0(x) in Eq. (21) must be subject to the side constraints (6).
If the interpolation data values , fN G^ at the point locations x1,—,xN	are given, the RBF
interpolant of 0(x) in Eq. (21) can be obtained by solving the system of N+2 linear equations for N+2 unknown generalized expansion coefficients:
4>(xi ) = f, i = I,---, N,
N	N	N
^= 0, ^aixi = 0, ^aiyi = 0,
which can be re-written in a matrix form as Ha= f, (23)
where
(22)
H=
A P PT 0
G^(n+3)>
(N+3)
(24)
A=
^1(xi) ••• ^N (xi) ^1(x N ) ••• ^N (x N )
G^NxN , (25)
P=
1 xi yi
1 xn yN
G^Nx3, (26)
a = [a • aN po P1 P2 ] G^Nx3,(27) f = ff • fN 0 0 0 f G^NX3, (28)
The generalized expansion coefficients a can be obtained by
a = H-1f, (29)
The resulting RBF interpolant of the implicit function in Eq. (21) can be re-written compactly as
4>(x) = Y T (x) a, (30)
where
Y(x) = f^>(x) - ^n (x) 1 x y f G^(NX3)X1, (31) Equation of motion
Since the Hamilton-Jacobi PDE (17) is time dependent, it is further assumed that all knots are fixed in space and the space and time are separable, and therefore the RBF interpolant of the implicit function in Eq. (30) becomes time dependent as
4>(x ,t ) = Y T (x )a(t), (32) Substituting Eq. (32) in (17) yields
Y T - a+ vn |(VY)Ta| = 0, (33)
where
G,1 =
G
-t
|(VY)r a = [(G»2 +(GT2a)2 j1'2, (34) " 010 g^(n+3)X1,
	
dx	dx
	
dy	dy
(35)
T
ACTA GEOTECHNICA SLOVENICA, 2007/i 49.
L. VRflNKflR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONSl IN THE MODËLING OF RflDIONUCLIDË MIGRATION AND
The initial value problem can be considered equivalent to the interpolation problem since the expansion coefficients at the initial time are found as a solution of the interpolation problem [15]. Therefore the preliminary starting point of the use of RBFs to solve PDEs is the interpolation problem that is equivalent to solving the initial value problem. The original equation (17) is thus converted into a time-dependent interpolation problem for the initial values of expansion coefficients and the propagation of the front is governed by the time dependent equation (33).
For time advance the initial values of a in Eq. (33) we used a collocation formulation of the method of lines. The governing equation of motion of the front (33) is extended to the whole domain O and the normal velocities vn at the front are thus replaced by the extension velocities in O. All nodes of domain are taken as fixed nodes of RBF interpolation. We also take into consideration constraints which must be introduced to guarantee that the generalized coefficients a can be solved.
Using the present collocation method for N points and above mentioned constraints, a set of resulting ODEs can be compactly written as:
H - a + B(a) = 0, (36) dt
where
vn (x1)|(v^r K))«!
= v'n(Xi)|(Wr(Xi))a eKNx3)xi, (37) 0 0 0
b) the exact explicit time integration. The first order homogeneous ODE of the form d a
B
The set of ODEs can be solved by several ODE solvers such as the first-order forward Euler's method and higher-order Runge-Kutta, Runge-Kutta-Fehlberg, Adams-Bashforth, or Adams-Moulten methods [16].
We used
a) the first-order forward Euler's method, an approximate solution to Eq. (36) is the following:
dt
+ Ea = 0, (39)
where
E = H-1(ViGi + v2G>2), (40) has the solution:
a(t + dt) = expm(-Edt)a(t - dt). (41)
The expression expm (-Edt) is a MATLAB exponential matrix function represents the series expansion or a rational fraction:
expm(-Edt) = I - Edt + (dt2 /2!)E2 - (dt3 /3!)E3 t • • •. (42)
E is the coefficient matrix resulting from the application of the Eq. (35).
4.2.3 Numerical example
For the simulation we used data from [11]: the concentration inside the part where the material characteristics remain constant upart = 0.53, the concentration on the interface usd = 0 , the initial concentration of the diffusive phase u° = 0.1, the diffusivity constant D = 1, the domain length l = 1 and the initial position of the interface s0 = 0.2. Let N be total number of points, r the number of points that lie inside a constant composition and N - r the number of points that lie inside the diffusive phase. Due to the movement of the interface, the point locations are adapted at each time step. The MQ exponent fi had the values 0.5 and 1.5. In Fig. 8 the movement of the interface positions calculated with a different MQ exponent, fi is presented. In the numerical experiments we compared our numerical solutions with the analytical solutions that exist for the problem presented in Chapter 4. (See [12]). The results are presented in Fig. 8 (see next page).
The next example is the rotation of the solid body. Consider the rotation of a circular bubble of the radius r =0.25 centered at (0.5, 0.15) in a vortex flow with the velocity field (v1, v2) = (-y, x). A half cycle of rotation is presented in Fig. 9.
a(tn+1) = a(tn ) - dtH-1B(a(tn ), where dt is the time step, and
(38)
50. ACTA GEOTECHNICA SLOVENICA, 2007/l
L. VRANKAR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONS) IN THE MODELING OF RADIONUCLIDE MIGRATION AND
Movement of the interface
0.45
-r
m 0.35 -
0.25 -
-analytic
— MQ(beta=0.5) O MQ(beta=1.5) LS M
0.4 0.5 Time
Figure 8. Interface position vs. time.
Figure 9. Zero contour of the level set function at different points and time during the rotation of a circle.
5 discussion and conclusion
In the case of the radionuclide migration (see also [6], [9] and [12]), two evaluation steps were performed. In the first step the velocities in principal directions were determined from the pressure of the fluid obtained from the Laplace differential equation. In the second step the advection-dispersion equation was solved to find the concentration of the contaminant. In this case the method of evaluation was verified by comparing results with those obtained from the finite difference method.
The traditional finite difference scheme was also used for solving the Laplace and advection-dispersion equation. For the approximation of the first derivative second-order central difference or one-sided difference were used. But for the approximation of the second derivatives we used the second-order central second difference. The time dependent part was implemented with the implicit scheme. The discretization grid has actually 12x12 points.
ACTA GEOTECHNICA SLOVENICA, 2007/i 51.
L. VRflNKflR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONSl IN THE MODELING OF RflDIONUCLIDe MIGRATION AND
Te results show that differences exist between both numerical schemes (Figs. 4 and 5). Different sets of input data yield differences between the schemes. In the simulation, a very large scatter caused by different given values of hydraulic conductivity was observed. Another reason for differences could result from the kriging method applied in the sequential Gaussian methods.
The comparison of concentrations calculated with the Eulerian and Lagrangian methods in partly heterogeneous porous (Figs. 6 and 7) media shows that the Lagrangian methods provide a wider concentration cloud in the area of high conductivity. It seems that it shows the influence of non-smooth change between low and high conductivity.
In general, the Eulerian approach is more convenient and more frequently used. But if it is important to study sharp changes (in our case between areas of low and high conductivity) of the solutions where important chemistry and physics take place, it is better to use the Lagrangian RBF scheme.
The comparison of the moving boundary positions calculated with the moving data centers method and MQ (0=0.5) and MQ (0=1.5) (Fig. 8) shows that MQ (0=1.5) determines the position of the interfaces much more accurately than MQ (0=0.5). The simulations have also shown that the value of the shape parameter c which was computed by the residual error procedure was in the range between 0.01 and 0.09. This confirms the fact that for a fixed number of centers N, smaller shape parameters produce more accurate approximations.
The comparison of the moving boundary positions calculated with the moving data centers method (MQ (0=1.5)) and the level set method (Fig. 8) also shows that the moving data centers method gives better results in this case. To achieve better accuracy, the resultant system of the RBF-PDE problem usually becomes badly conditioned. Several different strategies [17] have been somewhat successful in reducing the ill-conditioning problem when using RBF methods in PDE problems. These strategies include: variable shape parameters, domain decomposition, preconditioning of the interpolation matrix, and optimizing the center locations.
From chapter 4.2 we can see that RBFs (MQ) can be easily included in the level set formulation. Fig. 9 shows that we can get logical results using MQ in a 2-dimen-sional example.
We can conclude that the Kansa method is a valid alternative to the FDM because of its simple implementation and its easy use in the level set formulations. The only
geometric properties that are used in the RBF approximation are the pair-wise distances between points. Figs. 4 and 5 show that the RBF solution has far less diffusion than the finite difference method with the included upwinding.
In the future work we will use the Gershgorin circle theorem that could be useful a tool for choosing appropriate RBFs. For each value of the shape parameter, eigenvalues and their distribution can be studied to obtain knowledge concerning properties of an approximation matrix and their role being played in finding a better approximation of computed data to the equation solution. The solutions can be improved by using an affine space decomposition that decouples the influence between the interior and boundary collocations.
references
[1]	G. Fairweather, A. and Karageorghis (1998). The method of fundamental solutions for elliptic boundary value problems. Adv. Comput. Math., 9(1/2), 69-95.
[2]	Golberg, M.A. and Chen, C. S. (1998). The method of fundamental solutions for potential, Helmholtz and diffusion problems, in Boundary Integral Methods-Numerical and Mathematical Aspects, M. A. Golberg (Ed.). Computational Mechanics Publications, 103-176.
[3]	Kansa, E. J. (1990). Multiquadrics - A scattered data approximation scheme with applications to computational fluid dynamics - I. Surface approximations and partial derivative estimates. Computers Math. Applic. 19 (8/9), 127-145.
[4]	Kansa, E. J. (1990). Multiquadrics - A Scattered data approximation scheme with applications to computational fluid dynamics - II. Solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers Math. Applic. 19 (8/9), 147-161.
[5]	Hardy, R. L. (1971). Multiquadric equation of topography and other irregular surfaces. J. Geophys. Res. Vol. 176, 1905-1915.
[6]	Vrankar, L., Turk, G. and Runovc, F. (2004). Modelling of radionuclide migration through the geosphere with radial basis function method and geostatistics. Journal of the Chinese Institute of Engineers, Vol. 27, No. 4, 455-462.
[7]	Fedoseyev, A. I, Friedman, M. J. and Kansa, E. J. (2002). Improved multiquadric method for elliptic
52. ACTA GEOTECHNICA SLOVENICfl, 2007/l
L. VRANKAR & ET AL.: THE USE OF THE MESH FREE METHODS (RADIAL BASIC FUNCTIONS) IN THE MODELING OF RADIONUCLIDE MIGRATION AND
partial differential equations via PDE collocation on the boundary. Computers Math. Applic. Vol. 43, 439-455.
[8]	Smith, G. D. (1978). Numerical solution of partial differential equations: Finite difference methods. Oxford Applied Mathematics and Computing Science Series, Oxford University Press.
[9]	Vrankar, L., Turk, G. and Runovc, F. (2004). Combining the radial basis function Eulerian and Lagrangian schemes with geostatistics for modelling of radionuclide migration through the Geosphere. Computers. Math. Applic, Vol. 48, No. 10-11, 1517-1529.
[10]	Crank, J. (1984). Free and Moving Boundary Problems, Clarendon Press, Oxford.
[11]	Javierre, E., Vuik, C., Vermolen, F. J. and Van der Zwaag, S. (2005). A comparison of numerical models for one-dimensional Stefan problems, Reports of the Delft Institute of Applied Mathematics, Netherlands.
[12]	Vrankar, L., Turk, G., Runovc, F. and Kansa, E. J. (2006). Solving one-dimensional phase change problems with moving grid method and mesh free radial basis functions. International Conference Nuclear Energy for New Europe, Portorož. 1005.1.1005.8.
[13]	Osher, S. and Sethian, J. A. (1988). Front propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics. Vol. 78, 12--49.
[14]	Li, Z. and Ito, K. (2006). The Immersed Interface Method - Numerical Solutions of PDEs Involving Interfaces and Irregular Domains. Society forIin-dustrial and Applied Mathematics, Philadelphia, USA.
[15]	Wang, S. Y., Lim, K. M., Khoo, B. C. and Wang, M. Y. (2007). An extended level set method for shape and topology optimization. Journal of Computational Physics. Vol. 221, 395-421.
[16]	Greenberg, M. D. (1988). Advanced Engineering Mathematics, 2nd ed., Prentice-Hall International Inc., Upper Saddle River, NJ, USA.
[17]	Kansa, E. J. and Hon, Y. C. (2000). Circumventing the Ill - Conditioning problem with multiquadric radial basis functions: Applications to elliptic partial differential equations, Computers Math. Applic, Vol. 39, No. 7-8, 123-137.
ACTA GEOTECHNICA SLOVENICA, 2007/i 53.