Engineering Analysis with Boundary Elements 136 (2022) 77–92
A
0
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements
journal homepage: www.elsevier.com/locate/enganabound
A new radial basis integration method applied to the boundary element
analysis of 2D scalar wave equations
A. Narváez ∗, J. Useche
Universidad Tecnologica de Bolivar, Cartagena, Colombia
A R T I C L E I N F O A B S T R A C T
Keywords: A new integration method named the Radial Basis Integration Method (RBIM) that include the Kriging
Domain integration Integration Method (KIM) Narváez and Useche (2020) as a particular case and performs boundary only offline
Boundary element method precomputations for the creation of a meshless quadrature was developed for its application in boundary
Radial basis integration method elements. Herein, as in DR-BEM, the inertial term is approximated using radial basis functions, however, its
Dual reciprocity boundary element method
(DR-BEM) particular solution is not needed. The quadrature is now obtained in a simpler way than in KIM, because the
Scalar wave equation evaluations of domain integrals necessary to compute the weights of quadrature points, is done transforming
those to the boundary instead of using the Cartesian Transformation Method. Using RBIM, weakly singular
domain integrals may be computed with good accuracy over complex domains. The results obtained in some
scalar wave propagation problems using both Houbolt-𝛼 and Newmark-𝛼 time marching methods show that
this procedure can be even more accurate than other used in BEM analysis.1. Introduction dimensional inner integral. It has the advantage of being able to elim-
inate some types of singularities in the integrand [22], however, the
Scalar wave propagation problems are commonly encountered in method becomes computationally inefficient when the radial integral
various engineering disciplines. For the case of interior problems based cannot be computed analytically. LIM has the advantage that it only
on a direct BEM formulation in which a fundamental stationary solution needs to evaluate integrals on lines generated from integration points
is used in combination with direct time integration algorithms, several located in the boundary elements, although it must be combined with
strategies have emerged for the treatment of the inertial domain inte- Fast methods in order to be more efficient and approximate. DI-BEM
grals that appear there. The first method used to solve these integrals
is the Cell Integration Method (CIM), where a domain partition into unlike RI-BEM and LIM always allows the exact transformation to the
simpler regions called cells is made. Its use with boundary elements boundary of domain integrals using a primitive radial basis function
is known as classic BEM [1–3]. With the idea of obtaining integral (for which there are no restrictions on its choice) in combination with
formulations that only involve boundary integrals, new methods have Green’s theorem.
been developed, among these there are some with an approach similar On the other hand, within the category of methods based on
to CIM but with the difference that they do not require a mesh fitting meshless techniques, the following stand out: The Monte Carlo Method
task, rather they use a Cartesian grid for integration purposes [4–6]. In (MCM) [23,24], The Quasi-Monte Carlo Method [25],The Cartesian
the second instance, we can chronologically name those combined with Transformation Method (CTM) [26] and The Kriging Integration Method
BEM that convert domain integrals into boundary ones, among them we (KIM) [27]. All of them are based on a quadrature defined for inte-
have: The Dual Reciprocity Method (DR-BEM) [7–13] , the Triple Reci- gration points distributed within the domain, however, except in the
procity Method (TR-BEM) [14,15], the Multiple Reciprocity Method CMT such distributions can be random and allow solving problems with
(MR-BEM) [16–19], the Radial Integration Method (RI-BEM) [20], the singularities without any special treatment, although in KIM this only
Direct Interpolation with Boundary Elements (DI-BEM) [21] and more seems to be true for weak singularities. All these methods can evaluate
recently the Line Integration Method (LIM). The first three are based integrals on irregular domains with multiple boundaries. Moreover,
on the reciprocity theorem. In contrast to DR-BEM, MR-BEM only
needs boundary collocation points, however, its application is restricted regarding efficiency by taking into account the precision vs the number
to certain types of problems. RIM is characterized by converting the of integration points, KIM is as good as CTM even with fewer points and
domain integrals into boundary ones by previously carrying out an one therefore it saves computational time in the integral evaluation process.
∗ Corresponding author.
E-mail addresses: anarvaez@utb.edu.co (A. Narváez), juseche@utb.edu.co (J. Useche).https://doi.org/10.1016/j.enganabound.2021.12.005
Received 8 September 2021; Received in revised form 18 November 2021; Accepte
vailable online 10 January 2022
955-7997/© 2021 Elsevier Ltd. All rights reserved.d 6 December 2021
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
t
𝑢
𝑢
b
𝒖
A
g
In BEM the discretization process for the domain and/or the bound-
ary gives rise to a matrix equation of motion which is solved us-
ing several time marching stepping algorithms. When solving wave
propagation problems, Houbolt [28] and Newmark [29] are the most
common used methods. However, exist another stable and accurate
alternatives [2,30–33].
Making a review of the different existing domain integration meth-
ods in dynamic problems, it is evident that no extensive studies has
been carried out with RI-BEM to address scalar wave propagation
problems. Only has been reported a new integration procedure incor-
porated into RI-BEM (hereafter referred as MRI-BEM) to analyze some
problems in both convex and concave domains [34]. This approach
uses an auxiliary point, conveniently located so that the function to
be integrated in the radial direction is always evaluated over intervals
within the domain. Although what is observed is that there are more
application of RI-BEM related to acoustic eigenvalue and Helmholtz
problems [35–37]. Similarly, the Triple reciprocity method (TRM) and
DI-BEM has been applied to solve Helmholtz problems in [15,21],
respectively.
As it has been noted, there are few methods in BEM to compute
domain integrals that appears in applications concerning scalar wave
equation problems.
In this paper, a meshless-based integration method as in [26,27,38],
named the Radial Basis Integration Method (RBIM) has been developed
for the computation of inertial domain integrals that are found in the
BEM formulation of the 2D scalar wave equation. Using this quadrature,
it is possible to compute them directly, therefore, it is a potential
alternative in relation to the existent methods. The coupling between
RBIM and BEM will be named RBI-BEM.
The other sections of this article are organized as follows: In the
next section, the boundary integral representation of the scalar wave
equation is shown. In the third section, the basis of the Radial Basis
Integration Method is revised. In the fourth section, a review of the
modified version of Radial Integration Method is performed. In the
section fifth, in order to implement this method several subroutines
are given as possible options . In section six, RBI-BEM is applied to
solve some benchmark problems and comparisons with DR-BEM and
the MRI-BEM are made. Finally, in section seven, some conclusions are
presented.
2. Boundary integral representation of the scalar wave equation
The equation that governs transient scalar waves in the region 𝛺
discarding source terms is (Fig. 1):
2
∇2𝑢(𝒙 𝜕𝑢 (𝒙, 𝑡), 𝑡) = 1 (1)
𝑐2 𝜕𝑡2
where 𝑐 refers to the wave propagation velocity, 𝑢(𝒙, 𝑡) represent the
potential field variable , and 𝑡 is the time.
The time dependent boundary conditions are defined as follows:
𝑢(𝒙, 𝑡) = ?̄?(𝒙, 𝑡), on 𝛤𝑢 (2)
𝜕𝑢(𝒙, 𝑡)
= 𝑞(𝒙, 𝑡), on 𝛤𝑞 (3)𝜕𝑛
where 𝑛 stand for the boundary outward normal. Additionally we have
he following initial values:
(𝒙, 0) = 𝑢0(𝒙), on 𝛺 (4)
̇ (𝒙, 0) = ?̇?0(𝒙), on 𝛺 (5)
The boundary element formulation for Eq. (1) is defined as fol-
lows [39]:
𝑐(𝝃)𝑢(𝝃) + ∫ 𝑞
∗(𝝃,𝒙)𝑢(𝒙)d𝛤 = ∗∫ 𝑢 (𝝃,𝒙)𝑞(𝒙)d𝛤𝛤 𝛤
− 1 𝑢∗(𝝃,𝒙)?̈?(𝒙)d𝛺 (6)
𝑐2 ∫ 𝑯𝛺
78Fig. 1. Physical domain, its boundary and some point relationships.
where 𝒙 are field points, 𝝃 represent source points, while 𝑢∗ and
𝑞∗ represent the two well known fundamental solutions for the field
variable 𝑢 and its gradient, respectively [40].
After performing the spatial discretization process, the following
matrix equation is obtained:
𝑯𝒖𝒏+𝟏 = 𝑮𝒒
1
𝑛+1 − 𝑴?̈?𝒏+𝟏 (7)𝑐2
2.1. Time integration method (TIM)
The time integration methods used in this work are the proposed by
Mansur in [32]. These start redefining Eq. (7) as follows:
(1 + 𝛼)𝑯𝒖𝒏+𝟏 − 𝛼𝑯𝒖
1
𝒏 = 𝑮𝒒𝑛+1 − 𝑴?̈?𝑐2 𝒏+𝟏
(8)
The idea behind introducing the 𝛼 parameter in Eq. (8) is to get a
better control of the damping. Experience indicates that when assigning
positive values to such a parameter the damping is increased and
vice versa. The acceleration vector may be approximated with the
Newmark [29] or the Houbolt [28] algorithms.
In the Houbolt method, the following expressions are used for time
derivatives:
[ ]
?̇?𝒏+𝟏 =
1 11𝒖
6𝛥𝑡 𝒏+𝟏
− 18𝒖𝒏 + 9𝒖𝒏−𝟏 − 2𝒖𝒏−𝟐
[ ]
?̈? 1𝒏+𝟏 = 2𝒖𝒏+𝟏 − 5𝒖𝒏 + 4𝒖 − 𝒖𝛥𝑡2 𝒏−𝟏 𝒏−𝟐
For the case when 𝑛 = 0 the vectors 𝒖−𝟏 and 𝒖−𝟐 must be computed
based on a first order Euler scheme as:
𝒖−𝟏 = 𝒖𝟎 − ?̇?𝟎 × 𝛥𝑡
𝒖−𝟐 = 𝒖𝟎 − 2?̇?𝟎 × 𝛥𝑡
Besides, in the Newmark scheme, the following expressions are used
for time derivatives:
[ ]
𝛾 𝒖 − 𝒖
?̇? = 𝑛+1 𝑛 (𝛽 − 𝛾) (𝛾 − 2𝛽)𝑛+1 + ?̇?𝛽𝛥𝑡2 𝛽 𝒏
− ?̈?
2𝛽 𝑛
𝛥𝑡
[ ]
𝒖
?̈? 𝑛+1
− 𝒖𝑛 1 ?̇? (1 − 2𝛽)𝑛+1 = + 𝑛 − ?̈?𝛽𝛥𝑡2 𝛽𝛥𝑡 2𝛽 𝒏
At the beginning the accelerations in general are unknown and must
e computed as:
̈ = 2
[ ]
𝟎 𝒖𝟏 − 𝒖𝟎 − ?̇? 𝛥𝑡 (9)𝛥𝑡2 𝟎
fter having approximated the acceleration vector, the following
eneric expression is obtained:
̂ 𝒖 −𝑮𝒒 = 𝒈 (10)𝒏+𝟏 𝑛+1 𝒏
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92Table 1
Definitions for generic expression.
TIM ?̂? 𝒈𝒏
Houbolt-𝛼 𝑛 = 0 (1 + 𝛼)𝑯 + 22 2 𝑴 𝛼𝑯𝒖
2
𝑐 𝛥𝑡 𝟎 + 𝑐2 2 𝑴[𝒖𝟎 + ?̇?𝟎 × 𝛥𝑡]𝛥𝑡
Houbolt-𝛼 𝑛 > 0 (1 + 𝛼)𝑯 + 2 12 2 𝑴 𝛼𝑯𝒖𝒏 − 2 2 𝑴[−5𝒖 + 4𝒖𝑐 𝛥𝑡 𝑐 𝛥𝑡 𝒏 𝒏−𝟏 − 𝒖𝒏−𝟐]
Newmark-𝛼 𝑛 = 0 (1 + 𝛼)𝑯 + 2 22 2 𝑴 𝛼𝑯𝒖𝟎 + 2 2 𝑴[𝒖𝟎 + ?̇?𝟎 × 𝛥𝑡]𝑐 𝛥𝑡 𝑐 𝛥𝑡
Newmark-𝛼 𝑛 > 0 (1 + 𝛼)𝑯 + 12 2 𝑴 𝛼𝑯𝒖
1
𝒏 + 2 2 𝑴[𝒖𝒏 + ?̇?𝒏 × 𝛥𝑡 + ?̈?𝒏(0.5 − 𝛽) × 𝛥𝑡2]𝛽𝑐 𝛥𝑡 𝛽𝑐 𝛥𝑡Fig. 2. Integration based on whole domain.
Fig. 3. Integration based on auxiliary subdomains.
Depending on the method used to approximate the acceleration the
schemes may be called as Houbolt-𝛼 or Newmark-𝛼, respectively. Some
definitions for the Eq. (10) are listed in Table 1.
3. Radial basis integration method
The Radial Basis Integration Method (RBIM) uses a distributed set
of points inside the whole domain or inside auxiliary sub-domains as
integration points. In order to obtain greater precision it is recom-
mended that the points be distributed in a uniform or quasi-uniform
way, because the location of the singularity in the domain integrals
is variable on the domain. One possible way to obtain such points is
by using a bounding box and hexagonal patterns with a point at its
centroid. The points outside the domain must be rejected using an
appropriate algorithm (Figs. 2–3). Another alternative is to generate
quasi-random distributions as in mesh-free methods [25,41] .
The present method is based on our most recent work, i.e. the
Kriging Integration Method [27]. The main difference is that the new
approach contemplates the use of other radial basis functions different
to the Gaussian but without free parameters, and they are also comple-
mented with polynomial terms to ensure greater convergence. With the79absence of these parameters, a significant computation time saving can
also be achieved with similar results. On the other hand, the evaluation
of the domain integrals required to obtain the weights are now exactly
converted to boundary integrals instead of being evaluated using the
Cartesian Integration Method, which in general requires a significantly
higher number of integration points and more specialized programming
when controlling appropriate number of integration points in specific
areas in quite complex domains such as those containing hundreds of
arbitrary perforations of different sizes.
The theoretical foundations to get the weights of this quadrature is
very similar in both methods, however, for a better understanding it is
better to show it in detail as described below.
The weights of the quadrature points can be obtained by carrying
out a similar procedure as in [38,42], which is based on performing
certain matrix/vector operations taking into account the following
augmented radial basis approximation:
𝑛𝑟 𝑛𝑝
∑ ∑
𝑢(𝒙) ≈ 𝛼𝑖𝑢𝑖 + 𝛽𝑗𝑝𝑗 (𝒙) = 𝜱𝑇 ?̂? (11)
𝑖=1 𝑗=1
where 𝑛𝑟 is the amount of points employed with the RBF, 𝑛𝑝 refer to
the number of polynomial terms, 𝑢𝑖 are the values of the variable to be
approximated, 𝜱 = [𝛼1, 𝛼2,… ., 𝛼𝑛 , 𝛽1, 𝛽2,… ., 𝛽 ]𝑇𝑛 where 𝛼 and 𝛽 are𝑟 𝑝
unknown coefficients, while
?̂? = [𝑢1, 𝑢2, 𝑢3,… , 𝑢𝑛 , 0, 0, 0,… , 0]𝑇 whose size is𝑟
(𝑛𝑟 + 𝑛𝑝) × 1. The vector 𝜱 can be obtained as follows:
[ 𝑇 ]−1 [ ]
𝜱 = 𝐴−1𝒗 = 𝑪 𝑷 𝒗𝒓𝑷 𝟎 𝒗 (12)𝒑
where the elements of matrix 𝑪 are defined by 𝐶𝑖𝑗 = 𝐶(𝒙𝒊,𝒙𝒋) = 𝑓 (𝑟𝑖𝑗 )
in which 𝑓 is a function that belongs to a group of RBFs evaluated for
each point combination 𝑖 and 𝑗 (See Fig. 2).
𝒗𝒓 = [𝐶(𝒙𝟏,𝒙), 𝐶(𝒙𝟐,𝒙),… , 𝐶(𝒙𝒏 ,𝒙)]𝑇 , 𝒗𝒑 = 𝒗𝒑(𝒙) = [1, 𝑝1(𝒙), 𝑝2(𝒙),𝒓
… , 𝑝 (𝒙)]𝑇𝑛 and the matrix 𝑷 may be defined in compact form as𝑝
follows:
[ ]
𝑷 = 𝒗𝒑(𝒙𝟏), 𝒗𝒑(𝒙𝟐),… , 𝒗𝒑(𝒙𝒏 )𝒓 𝑛𝑝×𝑛𝑟
The required weights could be computed for the whole domain or
alternatively for auxiliary sub-domains. The integral of function 𝑢(𝒙)
over the 𝑘th auxiliary sub-domain using the present method may be
computed as follows:
𝑢(𝒙)𝑑𝛺 = 𝝓𝑇∫ 𝑘 ∫ ?̂?𝑑𝛺𝑘 = ∫ (𝐴
−1
𝑘 𝒗)
𝑇 ?̂?𝑑𝛺𝑘 (13)
𝛺𝑘 𝛺𝑘 𝛺𝑘
with the matrix 𝐴−1𝑘 defined as in Eq. (12) for the 𝑘th sub-domain .
However, the integral defined above can also be calculated according
to the following quadrature:
𝑛
⟨ ⟩ 𝑟
∑
∫ 𝑢(𝒙)𝑑𝛺 ≈ 𝒘𝜴 , ?̂? = ⟨𝒘𝒌, 𝒖⟩ = 𝑤
𝑖 𝑢
𝒌 𝑘 𝑖 (14)𝛺𝑘 𝑖=1
Therefore, it is inferred that
𝒘 = (𝐴−1𝜴 ∫ 𝑘 𝒗)𝑑𝛺
−1
𝑘 = 𝐴𝑘 ∫ 𝒗𝑑𝛺𝑘 (15)𝒌 𝛺𝑘 𝛺𝑘
or
[ ] [ 𝑻 ]𝒘 𝑪 (𝑷 ) −1
[ ]
𝒌 = 𝒌 𝒌 𝒗𝒓 𝑑𝛺𝑘𝒛𝒌 𝑷 𝒌 𝟎 ∫𝛺 𝒗𝑘 𝒑
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
w
u
∫
t
w
R
v
𝑭
𝐼
The domain integrals that appear in the Eq. (15) may be converted
into the boundary by using the Green Theorem. Therefore,
( 𝑥 )
∫ 𝐶𝑘(𝒙𝒏,𝒙)𝑑𝛺𝑘 = ∫ ∫ 𝐶𝑘(𝒙𝒏,𝒙)𝑑𝑥 𝑛𝑥𝑑𝛤𝑘𝛺𝑘 𝛤𝑘 𝑎 (16)
= ∫ 𝐷𝑘(𝒙𝒏,𝒙) 𝑛𝑥𝑑𝛤𝑘𝛤𝑘
( 𝑥 )
∫ 𝑃
𝑖 𝑖
𝑘(𝒙)𝑑𝛺𝑘 = ∫ ∫ 𝑃𝑘(𝒙)𝑑𝑥 𝑛𝑥𝑑𝛤𝑘𝛺𝑘 𝛤𝑘 𝑎 (17)
= 𝐸𝑖∫ 𝑘(𝒙) 𝑛𝑥𝑑𝛤𝑘𝛤𝑘
where 𝑎 is an arbitrary constant. As indicated in [26], only when the
interior integral is evaluated numerically is it necessary to set the value
of 𝑎 as follows:
∑𝑚
𝑗=1 𝑥𝑗𝑎 =
𝑚
here 𝑥𝑗 is the abscissa of the 𝑗th node of the discretized boundary
sing 𝑚 boundary element nodes.
Thus, the discretized form of Eqs. (16) and (17) are given by:
𝑙𝑘
∑
𝐶𝑘(𝒙𝒏,𝒙)𝑑𝛺𝑘 = ∫ 𝐷𝑘(𝒙𝒏,𝒙) 𝑑𝑦𝛺 𝑒𝑘 𝑒=1 𝛤𝑘 (18)
𝑙𝑘 +1 ( )
∑ 𝑑𝑦(𝜉)
= ∫ 𝐷
𝑘(𝒙𝒏,𝒙(𝜉)) 𝑑𝜉
𝑒=1 −1 𝑑𝜉
𝑙𝑘
∑
𝑃 𝑘(𝒙)𝑑𝛺 = 𝐸𝑘∫ 𝑖 𝑘 ∫ 𝑖 (𝒙) 𝑑𝑦𝛺 𝑒𝑘 𝑒=1 𝛤𝑘 (19)
𝑙𝑘 +1 ( )
∑ 𝑑𝑦(𝜉)
= 𝐸𝑘∫ 𝑖 (𝒙(𝜉)) 𝑑𝜉
𝑒=1 −1 𝑑𝜉
with ∑𝑚 𝑒, ∑𝑚 𝑒 and 𝑑𝑦(𝜉)
𝑑𝜙
∑
𝑥(𝜉) = 𝑗=1 𝜙𝑗𝑥𝑗 𝑦(𝜉) =
𝑚 𝑗 𝑒
𝑗=1 𝜙𝑗𝑦𝑗 = 𝑦 ,𝑑𝜉 𝑗=1 𝑑𝜉 𝑗
in which 𝜙𝑗 are shape functions.
Among the functions that we have been tested to get the weights of
he present method are :
𝑟𝑗
• 2The Gaussian 𝑓 = 𝑒−𝑐( 𝑟 )𝑗 𝑚𝑎𝑥 , with the parameters 𝑐 and 𝑟𝑚𝑎𝑥
defined as in [43].
• The thin plate spline 𝑓𝑗 = 𝑟2𝑗 𝑙𝑛(𝑟𝑗 ), which it was recommended
in [44]. It can also be expressed as 𝑟𝑓𝑗 = 𝑟𝑗 𝑙𝑛(𝑟
𝑗
𝑗 ) to avoid
numerical problems at 𝑟𝑗 = 0.
• The polyharmonic spline 𝑓 = 𝑟4𝑗 𝑗 𝑙𝑛(𝑟𝑗 ), which can also be ex-
pressed as 𝑟𝑓𝑗 = 𝑟3𝑗 𝑙𝑛(𝑟
𝑗
𝑗 ) to avoid numerical problems at 𝑟𝑗 =
0.
where subscript 𝑗 indicates that the function is centered at point 𝑗.
Actually the problem of evaluating the last two functions when 𝑟𝑗
tends to zero occurs only when the direct substitution 𝑟𝑗 = 0 is made,
however, lim (𝑟4𝑟 →0 𝑗 𝑙𝑛(𝑟𝑗 )) = 0, so for this case 𝑓𝑗 = 0. If the suggested𝑗
equivalent expression 𝑟𝑓 3 𝑗𝑗 = 𝑟𝑗 𝑙𝑛(𝑟𝑗 ) is not used, alternatively a very
small value such as 𝜖 = 2.2204 × 10−16 can be added to the logarithmic
term, i.e. 𝑓𝑗 = 𝑟4𝑗 𝑙𝑛(𝑟𝑗 + 𝜖) and the same value given by the limit would
be obtained.
It should be noted that in the present method the functions used to
calculate the weights of any distribution of points as part of a quadra-
ture that estimates the value of the integral of an unknown function
must have the property of being interpolant (which possesses the Kro-
necker delta property) since such weights are associated with the values
for which both functions coincide when evaluated at those locations.
The RBFs that are good candidates to get the best integration weights
are those for which the inverse of the interpolating matrix 𝐴𝑘 can be
calculated without inconvenience, that is, those that have a relatively
low condition number compared to the number of points used. It is also
important that such radial basis functions have adequate local behavior
80that guarantees convergence. In this sense, positive-defined compact
radial basis functions, including those based on semivariograms such
as the Gaussian and also conditionally positive-defined functions such
as the Thin Plate Spline (TPS) and the Polyharmonic Spline augmented
with polynomial terms, may be suitable for obtaining the weights.
On the other hand, to approximate the source term in the RBI-BEM
formulation, both interpolant or approximant meshless functions could
be used. If the source term has a known mathematical expression or
some values are known at some points (not necessarily at the poles),
then approximant functions such as the based on Moving Least Squares
(MLS), or any radial basis interpolant function can be used. However,
sometimes the source term is a function of the unknown variable of
the problem and even its derivatives and the meshless functions must
pass through the values of the unknown function at such collocation
points, and it is also required to impose boundary conditions, so in
this particular case they must be interpolant. In this sense, taking into
account that the radial basis functions can be used in both situations
and their derivatives are easier to calculate, they were chosen in this
work to approximate the source term as is done in DR-BEM. Unlike
the radial basis functions used to calculate the quadrature weights, the
functions used to approximate the source term do not necessarily have
to be positive-defined or with compact support.
In this work we used the Polyharmonic Spline 3 𝑟𝑓𝑗 = 𝑟𝑗 𝑙𝑛(𝑟
𝑗
𝑗 ) sup-
plemented with the monomial family {1, 𝑥, 𝑦} because with this, lower
error percentages were obtained against integration obtained analyti-
cally or by other methods than when using the Thin Plate Spline in
some tests. For the Polyharmonic Spline the antiderivative 𝐷𝑘(𝒙𝒏,𝒙)
can be integrated analytically after performing the variable changes
?̃? = 𝑥 − 𝑥𝑛 and ?̃? = 𝑦 − 𝑦𝑛. We can choose for convenience 𝑎 = 𝑥𝑛 in
order to get a shorter expression, so the new integration limits become 0
and ?̃?. Such integral was calculated with the symbolic algebra program
Wx-Maxima, obtaining the following simplified expression:
𝒙 𝒙 (225?̃??̃?
4 + 150?̃?3?̃?2 + 45?̃?5) ln (?̃?2 + ?̃?2)
𝐷𝑘( 𝒏, ) = 450
240?̃?5 arctan (?̃?∕?̃?) − 240?̃??̃?4 − 70?̃?3?̃?2 − 18?̃?5
+
450
On the other hand, each monomial could be transformed to the bound-
ary using a different value of 𝑎 because the integration is calculated
analytically and therefore its value can be arbitrary. With the idea of
getting the shortest expression of the integrand, 𝑎 = 0 was chosen for
such integrals, so we got the following terms:
2
𝐸1𝑘(𝒙) = 𝑥, 𝐸
2
𝑘(𝒙) =
𝑥 , 𝐸3
2 𝑘
(𝒙) = 𝑥𝑦
Then, using the previous expressions, the integrals in Eqs. (18) and
(19) can be calculated numerically to obtain the integration weights
from the matrix definition of Eq. (15).
All singular domain integrals present in boundary elements can be
expressed in discretized form as follows:
𝑁𝑆
∑
𝐼𝑖 = ∫ 𝑏(𝒙)𝑢
∗(𝝃𝒊,𝒙)d𝛺 ≡ ∫ 𝑏(𝒙)𝑢
∗(𝝃𝒊,𝒙)d𝛺𝑘 (20)
𝛺 𝑘=1 𝛺𝑘
where 𝑖 = 1, 2,… ,𝑀 . Here 𝑏(𝒙), 𝑢∗, 𝒙 and 𝝃𝒊 stand for the source
term, the fundamental solution, the field points and source points,
respectively. The source term may be approximated as in the Dual
reciprocity method as follows:
𝑚
∑
𝑏(𝒙) = 𝛼𝑗𝑓𝑗 (𝒙) = 𝒇𝑻 (𝒙)𝜶 (21)
𝑗=1
here 𝑚 represent the sum of boundary and domain nodes , 𝑓𝑗 are the
BFs for each application point and 𝛼𝑗 are unknown coefficients . In
ector notation the following definition are given:
𝒇𝑻 (𝒙) = [𝑓1(𝒙), 𝑓2(𝒙),… , 𝑓𝑚(𝒙)] and 𝜶 = 𝑭−𝟏𝒃 where
[ ]
= 𝒇 (𝒙𝟏),𝒇 (𝒙𝟐),… ,𝒇 (𝒙𝒎) . Therefore, Eq. (20) can be rewritten as:
𝑁𝑆
∑
= ∗𝑖 ∫ 𝑢 (𝝃
𝒊,𝒙)𝒇𝑻 (𝒙)d𝛺 𝑭−𝟏𝑘 𝒃𝑘=1 𝛺𝑘
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92Applying the present method we have:
[ ]
𝐼𝑖 = 𝑺𝒊𝟏 ,𝑺𝒊𝟐 ,… , 𝑺 −𝟏𝒊𝒋 ,… , 𝑺𝒊𝒎 𝑭 𝒃
where
𝑁𝑆
∑
𝑺 = 𝑓 ∗ 𝒊𝒊𝒋 ∫ 𝑗 (𝒙)𝑢 (𝝃 ,𝒙)d𝛺𝑘
𝑘=1 𝛺𝑘 (22)
𝑁𝑆 𝑁𝑆
∑ ∑
≈ 𝒘𝑇𝒌𝒇 𝒋(𝑸𝒌) ∙ 𝒖
∗(𝝃𝐢,𝐐𝐤) ≡ 𝒘𝑇𝒌𝒉𝒋(𝝃
𝐢,𝐐𝐤)
𝑘=1 𝑘=1
Here 𝒘𝒌 refer to the vector of the integration weights and the
symbol ∙ represents an element-wise product. The vector 𝒉𝒋 contains
the function values for each quadrature point inside the 𝑘th subdomain.
Unlike RI-BEM in RBI-BEM by applying Green’s theorem the inner
integral is computed with respect the abscissa coordinate, therefore 𝑅𝑗
it is not necessary to be defined as a function of 𝑟𝑖 (Fig. 1). In Eq. (22)
only if 𝑸𝒎 ≠ 𝝃𝐢𝒌 there is a contribution to the integral, being 𝑸
𝒎
𝒌 the 𝑚th
integration point in the 𝑘th subdomain.
The column vector of domain integrals may be expressed in compact
way as follows:
𝐼 = 𝑺𝑭−𝟏𝒃 (23)
where
⎡𝑺𝟏𝟏 𝑺𝟏𝟐 ⋯ 𝑺𝟏𝒎 ⎤
⎢ ⎥
𝑺 = 𝑺𝟐𝟏 𝑺𝟐𝟐 ⋯ 𝑺𝟐𝒎⎢ ⎥
⎢ ⋮ ⋮ ⋱ ⋮ ⎥
⎢𝑺 ⎥
⎣ 𝒎𝟏 𝑺𝒎𝟐 ⋯ 𝑺𝒎𝒎⎦
Following a similar procedure for the approximation of ?̈?(𝑡) in the
scalar wave equation we have:
𝐼 = 𝑺𝑭−𝟏?̈? =𝑴?̈? (24)
It must be highlighted that if the domain is partitioned into subdo-
mains we get the following advantages:
• The size of matrix 𝑨−𝟏 is reduced, thus, the calculation of the
weights is done more quickly.
• Different types of singularities could be treated by only applying
MRI-BEM or LIM over the auxiliary subdomains that represent the
near field of the singular point.
• It could be possible to accelerate RBI-BEM with ACA or Fast
Multipole Methods (FMMs).
3.1. Weakly singular domain integrals
In RBI-BEM, in the same way as in the Cartesian Integration Method
(CTM) [26] and the combined approach called DIBEM-RIM [37], reg-
ularization techniques can be used to treat weakly singular integrals,
however, experience has shown that the results obtained using regu-
larization can be even less precise than the achieved by the untreated
ones. Additionally, it was inferred from the mathematical development
carried out in [27] that the presence of logarithmic singularities does
not significantly affect the direct evaluation of domain integrals, there-
fore, in this work these integrals are calculated using a quadrature
obtained from the Radial Basis Integration Method and discarding the
possible integration point where the radial distance between this and
the source point becomes zero.
4. Review of the modified radial integration method
This is a modified version of the radial integration method. It was
specially developed to convert domain integrals into boundary ones
over non-convex domains. In general, the source term that appears in
domain integrals is approximated by RBFs, therefore, these integrals
are defined as follows:
𝑚
∑
𝐼 = 𝑏(𝒙)𝑢∗(𝝃𝑖𝑖 ∫ ,𝒙)d𝛺 = 𝛾𝑗𝐼𝑗 (𝝃
𝑖) (25)𝛺 𝑗=1
81Fig. 4. Physical domain and relationship among variables.
where 𝐼 (𝝃𝑖) = ∫ 𝝓 (𝑅 )𝑢∗𝑗 𝛺 𝑗 𝑗 (𝝃
𝑖,𝒙)𝑑𝛺, which is transformed by using
the Modified Radial Integration Method to the following equivalent
boundary integral [45]:
𝝃𝑖 𝜕𝜌
𝐹𝑗 (𝝃𝑖)𝐼𝑗 ( ) = ∫ 𝑑𝛤 (𝒙) (26)𝛤 𝜕𝑛 𝜌
in which
𝜚
𝐹𝑗 (𝝃𝑖) = ∫ 𝝓𝑗 (𝑅𝑗 )𝑢
∗(𝑟𝑖)𝜌𝑑𝜌 (27)
0
The distances 𝑟𝑖 and 𝑅𝑗 are expressed as function of the new variable
𝜚 from an auxiliary point as shown in Fig. 4. These distances are
calculated as follows:
√
𝑟𝑖 = 𝜚2 + 𝑙2 − 2𝜚𝑙𝑐𝑜𝑠(𝜑) (28)
√
𝑅 2 2𝑗 = 𝜚 + 𝑙 − 2𝜚𝑙𝑐𝑜𝑠(𝜃) (29)
The right side of Eq. (25) may alternatively be written in matrix
notation as:
[ ]
[ ]
𝐼 𝒊 𝒊 𝒊 −𝟏 𝒃𝑩𝑖 = 𝐼1(𝝃 ) , 𝐼2(𝝃 ) , … , 𝐼𝑚(𝝃 ) 𝑭 𝒃𝑫
5. Implementation
There are several possibilities to implement the present method,
among them we have:
5.1. Subroutine 1
1: Discretize the boundary of the whole domain.
2: Generate points inside a bounding box that circumscribes the do-
main, then search and store the 𝑁𝑝 integration points inside the do-
main. Alternatively, any algorithm/software as in meshfree meth-
ods to get the points may be used.
3: if 𝑁𝑝 > 2 then
4: Perform the Radial Basis Integration Method
5: Store the respective coordinates and their weights in a list to be
used in the RBI-BEM formulation.
6: end if
5.2. Subroutine 2
1: Perform a partition of the domain into auxiliary sub-domains by
intercepting background cells with the domain and discretizes its
boundaries.
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
5
5
o
a
a
6
m
t
t
s
u
g
A
o
e
T
b
t
n
t
p
s
i
2: Generate integration points inside a bounding box that circum-
scribes the whole domain.
3: for 𝑖 ← 1, 𝑁𝑆 do
4: Search and store the 𝑁𝑝 integration points inside the 𝑖th sub-
domain.
5: if 𝑁𝑝 > 2 then
6: Perform the Radial Basis Integration Method in the 𝑖th aux-
iliary sub-domain.
7: Store the respective coordinates and their weights in a multi-
dimensional list to be used in the RBI-BEM formulation.
8: end if
9: end for
.3. Subroutine 3
1: Perform a partition of the domain into auxiliary sub-domains by
intercepting background cells with the domain and discretizes its
boundaries.
2: for 𝑖 ← 1, 𝑁𝑆 do
3: Generate the integration points inside the 𝑖th subdomain using
any algorithm/software as in mesh free methods.
4: Perform the Radial Basis Integration Method in the 𝑖th auxiliary
sub-domain.
5: Store the respective coordinates and their weights in a multi-
dimensional list to be used in the RBI-BEM formulation.
6: end for
.4. Subroutine 4
1: Perform a partition of the domain by using a Centroidal Voronoi
Tessellation
2: for 𝑖 ← 1, 𝑁𝑆 do
3: Generate the integration points inside the 𝑖th subdomain using
any algorithm/software as in mesh free methods.
4: Perform the Radial Basis Integration Method in the 𝑖th auxiliary
sub-domain.
5: Store the respective coordinates and their weights in a multi-
dimensional list to be used in the RBI-BEM formulation.
6: end for
It is important to highlight that the quadrature is generated offline
nd used by the main code similarly as with any standard quadrature.
. Numerical examples
In the present work four examples were selected from some bench-
ark problems. For the generation of the integration quadrature with
he Radial Basis Integration Method, subroutine 1 was consider because
he examples addressed here do not belong to the category of large-
cale problems and relatively few elements and internal points were
sed. In order to get better results in the MRI-BEM four and six inte-
ration points were used for boundary and radial integrals respectively.
lthough other types of radial basis functions could be compared,
nly the well known linear and cubic RBFs, were tested in several
xamples using both domain and boundary nodes for its construction.
he singular boundary integrals were computed with the usual rigid
ody motion approach.
Here as in [26,27] the ratio ℜ = 𝑑∕𝑑𝑖𝑛𝑡 was used as a measure of
he appropriate quantity of Radial Basis Integration points. Where the
umerator represent the mean distance among all internal points, while
he denominator is the mean distance of the Radial Basis Integration
oints.
The average relative error 𝑒𝑟 between the present method and others
olutions, taking into account all time steps within the used time
nterval is computed based on the norm 𝐿2.
82Fig. 5. Square domain.
Table 2
Results at point 𝐴 for example 1.
Approach 𝑓𝑗 Houbolt-𝛼 Newmark-𝛼
𝑒𝑟 t(s) 𝑒𝑟 t(s)
DR-BEM 0.043 7.88 0.044 9.27
RBI-BEM 1 + 𝑟𝑗 0.057 9.57 0.054 9.64
MRI-BEM 0.081 845 0.059 812
In Tables 2–5 the term 𝑡(𝑠) stand for the execution time in seconds
f the respective BEM approaches codes.
For comparison purposes, we tried to maintain the same parameters
nd 𝛥𝑡 for the same time marching method in the different BEM
approaches, however, due to the instability in some methods some
modifications were necessary. Which is discussed at the end of the
paper.
6.1. Example 1
A unitary square region representing a clamped membrane with the
below prescribed initial conditions is considered (Fig. 5).
𝑢(𝑥, 𝑦, 0) = 𝑥(𝑥 − 1)𝑦(𝑦 − 1) (𝑥, 𝑦) ∈ 𝛺 (30)
?̇?(𝑥, 𝑦, 0) = 0 (𝑥, 𝑦) ∈ 𝛺 (31)
The respective theoretical solution may be found in [3]. Here the
wave propagation velocity was taken as 1 m∕s. The variation of 𝑢 at
membrane’s center was gotten using 80 constant boundary elements
and 81 domain points in combination with Houbolt-𝛼 and Newmark-𝛼
for a normalized time inside the interval [0, 5]. For the present method
the integration points were generated using the NodeLab Matlab’s
package considering ℜ = 1.98. Comparison of potential obtained from
the present method, DR-BEM, MRI-BEM and the theoretical solution
using the linear function 𝑓𝑗 = 1+ 𝑟𝑗 are provided in Figs. 6–7, in which
the results for all compared methods are shown to be in close agreement
with the obtained theoretically. According to Table 2, DR-BEM was
slightly more accurate than RBI-BEM, while that with MRI-BEM the
less approximate solution was obtained. The achieved results for RBI-
BEM and MRI-BEM with the Newmark-𝛼 method were better than with
the Houbolt-𝛼 method, however for both time marching schemes used
with DR-BEM the solutions were practically the same. Additionally,
no significant difference is observed in the execution times between
DR-BEM and RBI-BEM, but the same cannot be said for the difference
between these two and MRI-BEM.
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92Fig. 6. Displacement at point (0.5, 0.5) with Houbolt-𝛼 method.Fig. 7. Displacement at point (0.5, 0.5) with Newmark-𝛼 method.Fig. 8. Square domain.83Table 3
Results at point 𝐴 for example 2.
Approach 𝑓𝑗 Houbolt-𝛼 Newmark-𝛼
𝑒𝑟 t(s) 𝑒𝑟 t(s)
DR-BEM 0.110 26.86 0.132 10.33
RBI-BEM 1 + 𝑟𝑗 0.108 112.42 0.116 156.3
MRI-BEM 0.84 25995.3 0.854 28594
DR-BEM 0.804 27.07 0.88 13.16
RBI-BEM 𝑟3𝑗 0.214 294.27 0.14 273.3
MRI-BEM 0.932 15640 0.94 21788
6.2. Example 2
Here as in the before example both the same geometry with clamped
sides and propagation wave velocity were considered but in this case
under others initial conditions which are given as (Fig. 8):
𝑢(𝑥, 𝑦, 0) = 0 (𝑥, 𝑦) ∈ 𝛺 (32)
?̇?(𝑥, 𝑦, 0) = 𝑉0 = 1 (𝑥, 𝑦) ∈ 𝛺0 (33)
where 𝛺0 belongs to the domain center. The analytical solution is
shown in [2].
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
Fig. 9. Displacement at point (0.5, 0.5) with 𝑓𝑗 = 1 + 𝑟𝑗 and Houbolt-𝛼 method.
Fig. 10. Displacement at point (0.5, 0.5) with 𝑓𝑗 = 1 + 𝑟𝑗 and Newmark-𝛼 method.
Fig. 11. Displacement at point (0.5, 0.5) with 𝑓𝑗 = 𝑟3𝑗 and Houbolt-𝛼 method.
84
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92Fig. 12. Displacement at point (0.5, 0.5) with 𝑓 3𝑗 = 𝑟𝑗 and Newmark-𝛼 method.Fig. 13. Geometry and boundary condition data for example 3.
For this problem the initial non-zero velocity conditions were im-
posed directly on the poles belonging to 𝛺0. The variation of 𝑢 at the
center of membrane was obtained using 80 constant boundary elements
plus 625 domain nodes. This problem was solved with CIM in [2]
and using a mesh that almost triples the amount of nodes used here
, however, this nodal discretization becomes relatively prohibitive for
MRI-BEM due to the high computation time required to calculate the
integrals. For the present method the integration points were generated
using the NodeLab Matlab’s package for which the ratio ℜ = 1.92
was obtained. The displacement solution at the center of membrane
was obtained using Houbolt-𝛼 and Newmark-𝛼 methods. Comparison
of potential obtained from the present method, DR-BEM, MRI-BEM and
the theoretical solution are provided in Figs. 9–10 and Figs. 11–12
respectively. The results showed for DR-BEM and RBI-BEM in Fig. 9–10
are very similar to those obtained by the theoretical solution. However,
as shown in Figs. 11–12 bad results were obtained when using 𝑓𝑗 = 𝑟3𝑗
in DR-BEM which is related to instability and convergence problems.
It is evident in Table 3 that RBI-BEM is the most accurate of all. DR-
BEM was slightly less accurate than RBI-BEM for 𝑓𝑗 = 1+ 𝑟𝑗 using both
time marching schemes, while the performance of MRI-BEM was very85Table 4
Results at point 𝐴 for example 3.
Approach 𝑓𝑗 Houbolt-𝛼 Newmark-𝛼
𝑒𝑟 t(s) 𝑒𝑟 t(s)
DR-BEM 0.47 2.34 0.359 3.25
RBI-BEM 1 + 𝑟𝑗 0.273 6.48 0.309 8.23
MRI-BEM 0.217 390 0.215 192
DR-BEM 0.82 2.39 0.785 2.62
RBI-BEM 𝑟3𝑗 0.153 6.48 0.167 8.27
MRI-BEM 0.237 393.4 0.287 298.35
unsatisfactory in all cases, since a high level of average relative error
in the results and a very high execution time were achieved.
6.3. Example 3
In this example a concave L-shaped domain under mixed boundary
conditions is considered (Fig. 13). Here c=1 m/s was adopted.
In this problem we used 56 constant boundary elements plus 50
domain nodes. It was also solved by MRI-BEM in [34] using a similar
discretization. The solution for 𝑢 at point A(0.05, 0.25) in the interval
0 < 𝑡 ≤ 15 𝑠 was obtained using the Newmark-𝛼 and the Houbolt-𝛼
methods. Comparisons among present method, DR-BEM, MRI-BEM and
Comsol using a relative fine mesh are provided in Figs. 14–17. In which
it is possible to see that only the curves of RBI-BEM and MRI-BEM are
in close agreement with those obtained by Comsol.
From Table 4 it is possible to make the following inferences:
• MRI-BEM performed better than RBI-BEM when using 𝑓𝑗 = 1+ 𝑟𝑗 .
• In RBI-BEM the standard Houbolt (𝛼 = 0) and the standard
Newmark (𝛼 = 0) without numerical damping (𝛾 = 1∕2, 𝛽 = 1∕4)
were used, while that combining MRI-BEM with 𝑓𝑗 = 𝑟3𝑗 and the
Newmark-𝛼 different parameters had to be set in order to get
stability and accuracy.
• The performance of RBI-BEM was notably superior to the MRI-
BEM when 𝑓 3𝑗 = 𝑟𝑗 was employed.
• The performance of DR-BEM was unsatisfactory in all tested cases.
6.4. Example 4
Here a domain with a more complicated geometry whose boundary
is clamped was selected taking as reference the problem analyzed
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
Fig. 14. L-shaped domain under mixed boundary condition: Potential at point 𝐴(0.05, 0.25) with 𝑓𝑗 = 1 + 𝑟𝑗 and the Houbolt-𝛼 method.
Fig. 15. L-shaped domain under mixed boundary condition: Potential at point 𝐴(0.05, 0.25) with 𝑓𝑗 = 𝑟3𝑗 and the Houbolt-𝛼 method.
Fig. 16. L-shaped domain under mixed boundary condition: Potential at point 𝐴(0.05, 0.25) with 𝑓𝑗 = 1 + 𝑟𝑗 and the Newmark-𝛼 method.
86
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
𝑢
Fig. 17. L-shaped domain under mixed boundary condition: Potential at point 𝐴(0.05, 0.25) with 𝑓 3𝑗 = 𝑟𝑗 and the Newmark-𝛼 method.Table 5
Results at point 𝐴 for example 4.
Approach 𝑓𝑗 Houbolt-𝛼 Newmark-𝛼
𝑒𝑟 t(s) 𝑒𝑟 t(s)
DR-BEM 0.655 43.5 0.788 10.31
RBI-BEM 1 + 𝑟𝑗 0.625 61.7 0.753 45.12
MRI-BEM 0.537 6275.87 0.593 9720
DR-BEM 0.97 9.21 1.03 13.61
RBI-BEM 𝑟3𝑗 0.245 121.2 0.252 104.65
MRI-BEM 0.702 6851 0.738 4888
in [46] (Fig. 18). The closed region 𝛺 is bounded by the parametric
curve:
𝑥 = 𝑟0 cos 𝜃 , 𝑦 = 𝑟0 sin 𝜃 , 𝜃 ∈ [0, 2𝜋) (34)
where 𝑟0 is defined as:
[ √ ]1∕3
𝑟0 = cos(4𝜃) + 3.6 − sin
2(4𝜃) (35)
The initial values are defined using the following expressions:
𝑢(𝑥, 𝑦, 0) = 𝑒
−0.16(𝑥2+𝑦2)
(𝑥, 𝑦) ∈ 𝛺 (36)
10
̇ (𝑥, 𝑦, 0) = 0 (𝑥, 𝑦) ∈ 𝛺 (37)
The variation of 𝑢 at the center of membrane was obtained using 60
constant boundary elements plus 405 domain nodes. A similar scaled
problem using 60 boundary points and 1205 interior nodes was ana-
lyzed employing a Meshless approach in [46], however, for comparison
purposes, only 405 internal nodes were used in our analysis, taking into
account that MRI-BEM requires high computation time. The NodeLab
Matlab’s package was used in the present method for the generation of
integration points setting a ratio ℜ = 2. As in the previous examples
the same time marching stepping schemes were used. Due to the ana-
lytical solution is unavailable, all the BEM approaches were compared
with results obtained with Comsol, using 17092 triangular quadratics
elements and 34483 nodal points. Results are depicted in Figs. 19–22.
From Table 5 it is possible to make the following inferences:
• MRI-BEM performed better than RBI-BEM for 𝑓𝑗 = 1 + 𝑟𝑗 using
both time marching schemes.
• The performance of RBI-BEM was notably superior to the MRI-
BEM when 𝑓 = 𝑟3 was employed, although with the standard𝑗 𝑗
87Fig. 18. Irregular shape with a parametric based boundary.
Table 6
Minimum time step required for stability for all examples using Houbolt-𝛼.
𝑓𝑗 Example Approach
DR-BEM MRI-BEM RBI-BEM
1 0.0005 s 0.0005 s 0.0005 s
1 + 𝑟𝑗 2 0.0005 s 0.001 s 0.0005 s
3 0.15 s 0.0005 s 0.0005 s
4 0.07 s 0.01 s 0.2 s
1 0.0005 s 0.0005 s 0.0005 s
𝑟3𝑗 2 0.07 s 0.01 s 0.0005 s
3 0.85 s 0.0005 s 0.0005 s
4 1.0 0.03 s 0.007 s
Newmark method (𝛼 = 0), the parameters 𝛾 = 0.53 and 𝛽 = 0.3
had to be set in order to get more accuracy .
• The performance of DR-BEM was unsatisfactory in all tested cases,
and several set of parameters had to be set with the standard
Newmark (𝛼 = 0) in order to achieve stability and better results.
• MRI-BEM is very time consuming compared with DR-BEM and
RBI-BEM.
It can be highlighted that the present method achieved a similar
accuracy than Comsol using only few collocation points.
Although it seems unfair to compare several cases of the different
approaches under different time steps, the reasons are actually due
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
Fig. 19. Irregular domain with a smooth edge: Potential at point 𝐴(0, 0) with 𝑓𝑗 = 1 + 𝑟𝑗 and the Houbolt-𝛼 method.
Fig. 20. Irregular domain with a smooth edge: Potential at point 𝐴(0, 0) with 𝑓 3𝑗 = 𝑟𝑗 and the Houbolt-𝛼 method.
Fig. 21. Irregular domain with a smooth edge: Potential at point 𝐴(0, 0) with 𝑓𝑗 = 1 + 𝑟𝑗 and the Newmark-𝛼 method.
88
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
Fig. 22. Irregular domain with a smooth edge: Potential at point 𝐴(0, 0) with 𝑓 3𝑗 = 𝑟𝑗 and the Newmark-𝛼 method.
Fig. 23. Variation of condition number for matrix 𝐾 vs time step using the Houbolt-𝛼 method in example 1.
Fig. 24. Variation of condition number for matrix 𝐾 vs time step using the Houbolt-𝛼 method in example 2.
89
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92Fig. 25. Variation of condition number for matrix 𝐾 vs time step using the Houbolt-𝛼 method in example 3.Fig. 26. Variation of condition number for matrix 𝐾 vs time step using the Houbolt-𝛼 method in example 4.Table 7
Minimum time step required for stability for all examples using Newmark-𝛼.
𝑓𝑗 Example Approach
DR-BEM MRI-BEM RBI-BEM
1 0.0005 s 0.0005 s 0.0005 s
1 + 𝑟𝑗 2 0.001 s 0.005 s 0.0005 s
3 0.09 s 0.0005 s 0.0005 s
4 0.1 s 0.04 s 0.1 s
1 0.0005 s 0.0005 s 0.0005 s
𝑟3𝑗 2 0.07 s 0.014 s 0.0005 s
3 0.3 s 0.0008 s 0.0005 s
4 0.8 0.07 s 0.02 s
to instability issues (see Tables 6–7), since it was not possible to
find a more precise or convergent solution by decreasing the time
step, maintaining the same parameters of the analyzed time marching
integration methods, nor making a variation of them. For this reason,
only the best results were shown in this work.90It is possible to find justifications for the poor performance of DR-
BEM or MRI-BEM in combination with the radial basis functions used
for certain examples by reviewing the following: After substituting the
boundary conditions, the resultant matrix of the system of equations is
obtained, henceforth called K, which also depends on values calculated
from the inverse of the interpolation matrix F. It has been observed
that if the condition number of F is relatively high (which increases
with the number of points used to interpolate, and it was also reported
in [35]), the condition number of matrix K is also affected. In order to
find out the difference between the curves obtained for the condition
numbers vs the time step for the different formulations, Figs. 23–26
were plotted. There it is observed that the condition number of the
matrix K varies with the value of 𝛥𝑡 and its tendency is generally
decreasing for higher values of 𝛥𝑡. The values of the condition numbers
are similar between the different methods, except when DR-BEM and
the cubic function are combined. The same behavior was observed
when using Newmark-𝛼 instead. Furthermore, it is observed that the
values of 𝛥𝑡 used to obtain the lowest possible average relative error
in the results for the different examples generally corresponded to the
highest condition number values.
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92
7
b
b
a
b
a
f
c
i
i
In the Tables 6–7 it can be seen that both time marching schemes
can be unstable for certain time steps, possibly influenced by the ill-
conditioning of the matrices (which worsens when the number of poles
is significantly increased) in combination with certain inadequate radial
basis function for the problem under analysis. Regarding the size of the
intervals of 𝛥𝑡 where the solution was stable, there it is also observed
that with the RBI-BEM formulation the largest ranges were obtained on
average, while with DR-BEM the smallest ones.
In the analysis of the different problems, there were two situations
where it was not possible to achieve good results for the different
approaches: The first one was when using the shortest time step within
the range where there was stability. This is because that value does
not necessarily coincide with the appropriate 𝛥𝑡 according to the wave
propagation speed, and therefore, the precision of the time marching
method is reduced and deviated results are obtained concerning refer-
ence values. The second case occurred even when using an appropriate
time step, as in the second example when using MRI-BEM. In this
case, we can affirm that it may be due to the inability of such radial
basis functions to adequately represent the variation of the variable
to be approximated under the particular conditions of that problem
(Discontinuous initial conditions for the velocity), without neglecting
the effect that it can cause matrix ill-conditioning. The above can be
supported by what was found in the following works:
• In [47] the Thin Plate Spline radial basis function was used in all
examples because it was the only one that produced good results
with the DR-BEM. Here the following is also stated : ‘‘An excessive
number of poles can disturb the results, since the DR-BEM is not a
simple interpolation technique and conditioning matrix problems
and inadequacy of certain radial functions are related in many
computational tests’’.
• In [48] the function 𝑟3𝑗 was proved to be unpredictable in the
DR-BEM. Similarly in [49] it was found that not all RBFs are
appropriate to solve a certain problem with the DR-BEM.
• In [34] it is not explained for what reasons the function 𝑟3𝑗 was
only used for all the examples and therefore it is not known if
there are problems with the use of other functions in MRI-BEM.
On the other hand, the analyzes on the influence of the number of
internal points on the quality of the solution that has been made
in works where MRI-BEM has been implemented, only a maxi-
mum of 66 internal points have been considered, not relatively
large quantities such as the used in the second and fourth example
of our work, which were 625 and 405 respectively. Although
the results of other simulations were not included in the present
paper, it is good to highlight that for the second example no
better results were obtained with MRI-BEM when using a smaller
number of internal points (125 or 441) nor when using a larger
quantity (1089).
• According to [50] the DI-BEM performance was found to be
sensitive to the kind of RBF used to interpolate the source term.
They got the worst results using the function 1 + 𝑟3𝑗 with DI-BEM
to solve a certain 3D model , while its use with DR-BEM had been
a relative success in a reference paper.
. Conclusions
In the present article, a promising domain integration quadrature-
ased method was used to address wave propagation problems with the
oundary elements. In general, the condition number of the matrices
lone is not decisive to know if good results will be obtained, since
oth the type of BEM formulation (DR-BEM, MRI-BEM, RBI-BEM, etc.)
nd the specific conditions of the problem can influence. The RBI-BEM
ormulation appears to be more stable when using both the linear and
ubic RBFs and a relatively high quantity of internal points in compar-
son with the DR-BEM and the MRI-BEM ones. Other simulations not
2ncluded in this work carried out with functions 𝑓𝑗 = 𝑟𝑗 𝑙𝑜𝑔(𝑟𝑗 ) and 𝑓𝑗 =
91𝑟4𝑗 𝑙𝑜𝑔(𝑟𝑗 ) also produced very good results with RBI-BEM, however, the
performance of the present approach could be sensitive to some types of
radial basis function and therefore, it is necessary to review the stability
and quality of the solutions when using other global functions or those
with compact support. Which will be addressed in future research.
Although the results for the function 𝑟3𝑗 using both time marching
schemes in the first example were not shown, all the formulations
were in good agreement with the analytical solution, therefore, it is
confirmed that the type of problem also influences the performance of
the method used to transform the domain integral. On the other hand,
by using the present approach, the weakly singular domain integrals
appearing in BEM formulations may be evaluated directly without
using any treatment for the integrand. The present approach lacks some
limitations found in DR-BEM in relation to the choice of approximating
RBFs and is less time-consuming than the MRI-BEM. It can be used
to solve problems with simple or multi-connected convex and concave
geometries. Both Houbolt-𝛼 and Newmark-𝛼 methods are suitable to
solve this kind of problem, however, Newmark-𝛼 has a greater number
of control parameters that could be helpful in certain cases. The present
method may be implemented in BEM for solving several linear and
nonlinear engineering problems.
Acknowledgment
The authors gratefully acknowledges the Technological University
of Bolívar for supporting this work.
References
[1] Dallner R, Kuhn G. Efficient evaluation of volume integrals in the boundary
element method. Comput Methods Appl Mech Engrg 1993;109(1):95–109. http:
//dx.doi.org/10.1016/0045-7825(93)90226-N, URL https://www.sciencedirect.
com/science/article/pii/004578259390226N.
[2] Carrer JAM, Mansur WJ, Vanzuit RJ. Scalar wave equation by the boundary
element method: a D-BEM approach with non-homogeneous initial conditions.
Comput. Mech 2009;44(1):31. http://dx.doi.org/10.1007/s00466-008-0353-4.
[3] Carrer JAM, Mansur WJ. Scalar wave equation by the boundary element method:
a D-BEM approach with constant time-weighting functions.. Internat J Numer
Methods Engrg 2010;81(10):1281–97. http://dx.doi.org/10.1002/nme.2732.
[4] Ding J, Ye W. A grid-based integral approach for quasilinear problems. Comput
Mech 2006;38(2):113–8.
[5] Koehler M, Yang R, Gray LJ. Cell-based volume integration for boundary integral
analysis. Internat J Numer Methods Engrg 2012;90(7):915–27.
[6] Wang Q, Zhou W, Cheng Y, Ma G, Chang X, Huang Q. An adaptive cell-based
domain integration method for treatment of domain integrals in 3D boundary
element method for potential and elasticity problems. Acta Mech Solida Sin
2017;30(1):99–111. http://dx.doi.org/10.1016/j.camss.2016.08.002, URL https:
//www.sciencedirect.com/science/article/pii/S0894916616300751.
[7] Nardini D, Brebbia C. A new approach to free vibration analysis using boundary
elements. Appl Math Model 1983;7:157–62.
[8] Partridge PW, Brebbia CA. Computer implementation of the BEM dual reciprocity
method for the solution of general field equations. Commun Appl Numer Methods
1990;6(2):83–92. http://dx.doi.org/10.1002/cnm.1630060204.
[9] Agnantiaris J, Polyzos D, Beskos D. Some studies on dual reciprocity BEM for
elastodynamic analysis. Comput Mech 1996;17(4):270–7.
[10] Partridge PW, Brebbia CA, et al. Dual reciprocity boundary element method.
Springer Science & Business Media; 2012.
[11] Useche J, Narvaez A. Vibrational response of elastic membranes coupled to
acoustic fluids using a BEM–BEM formulation. In: De Clerck J, editor. Topics
in modal analysis I, vol. 7. Cham: Springer International Publishing; 2014, p.
333–40.
[12] Useche J. Vibration analysis of shear deformable shallow shells using the
boundary element method. Eng Struct 2014;62–63:65–74. http://dx.doi.org/10.
1016/j.engstruct.2014.01.010.
[13] Useche J, Alvarez H. Elastodynamic analysis of thick multilayer compos-
ite plates by the boundary element method. CMES-Comput Model Eng Sci
2015;107(4):277–96.
[14] Ochiai Y, Sladek V, Sladek J. Transient heat conduction analysis by triple-
reciprocity boundary element method. Eng Anal Bound Elem 2006;30:194–204.
http://dx.doi.org/10.1016/j.enganabound.2005.07.010.
[15] Guo S, Wu Q, Li H-G, Kuidong G. Triple reciprocity method for unknown
function’s domain integral in boundary integral equation. Eng Anal Bound Elem
2020;113:170–80. http://dx.doi.org/10.1016/j.enganabound.2019.12.014.
A. Narváez and J. Useche Engineering Analysis with Boundary Elements 136 (2022) 77–92[16] Nowak A, Brebbia C. The multiple-reciprocity method. a new approach for
transforming bem domain integrals to the boundary. Eng Anal Bound Elem
1989;6(3):164–7.
[17] Nowak A. The multiple reciprocity boundary element method. Southampton:
Computational Mechanics Publication; 1994, p. 1–41.
[18] Ochiai Y, Sekiya T. Steady thermal stress analysis by improved multiple-
reciprocity boundary element method. J. Therm Stresses 1995;18(6):603–20.
[19] Ochiai Y, Sekiya T. Steady heat conduction analysis by improved multiple-
reciprocity boundary element method. Eng Anal Bound Elem 1996;18(2):111–7.
[20] Gao X-W. The radial integration method for evaluation of domain integrals with
boundary-only discretization. Eng Anal Bound Elem 2002;26(10):905–16.
[21] Loeffler C, Cruz A, Bulcão A. Direct use of radial basis interpolation functions
for modelling source terms with the boundary element method. Eng Anal Bound
Elem 2015;50:97–108. http://dx.doi.org/10.1016/j.enganabound.2014.07.007.
[22] Gao X-W. Evaluation of regular and singular domain integrals with
boundary-only discretization—theory and fortran code. J Comput Appl Math
2005;175(2):265–90.
[23] Gipson G. The coupling of Monte Carlo integration with boundary integral
techniques to solve Poisson-type problems. Eng Anal 1985;2(3):138–45.
[24] Kagawa Y, Murai T. Use of Monte Carlo method for singular integral evaluation
in boundary elements. COMPEL 1991.
[25] Rosca VE, ao VML. Quasi-Monte Carlo mesh-free integration for meshless weak
formulations. Eng Anal Bound Elem 2008;32(6):471–9. http://dx.doi.org/10.
1016/j.enganabound.2007.10.015, Meshless Methods.
[26] Hematiyan M. A general method for evaluation of 2D and 3D domain inte-
grals without domain discretization and its application in BEM. Comput Mech
2007;39(4):509–20.
[27] Narváez A, Useche J. The kriging integration method applied to the boundary
element analysis of Poisson problems. Eng Anal Bound Elem 2020;121:1–20.
http://dx.doi.org/10.1016/j.enganabound.2020.09.001.
[28] Houbolt JC. A recurrence matrix solution for the dynamic response of elastic
aircraft. J Aeronaut Sci 1950;17(9):540–50.
[29] Newmark NM. A method of computation for structural dynamics. J Eng Mech
Div 1959;85(3):67–94.
[30] Yu G, Mansur WJ, Carrer JAM, Gong L. A linear 𝜃 method applied to
2D time-domain BEM. Commun Numer Methods Eng 1998;14(12):1171–
9. http://dx.doi.org/10.1002/(SICI)1099-0887(199812)14:12<1171::AID-
CNM217>3.0.CO;2-G, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/.
[31] Chen W, Tanaka M. A study on time schemes for DRBEM analysis of elastic
impact wave. Comput Mech 2002;28:331–8. http://dx.doi.org/10.1007/s00466-
001-0297-4.
[32] Carrer J, Mansur W. Alternative time-marching schemes for elastodynamic
analysis with the domain boundary element method formulation. Comput Mech
2004;34:387–99. http://dx.doi.org/10.1007/s00466-004-0582-0.
[33] Oyarzún P, Loureiro F, Carrer J, Mansur W. A time-stepping scheme based
on numerical green’s functions for the domain boundary element method: The
exga-DBEM newmark approach. Eng. Anal. Bound. Elem. 2011;35:533–42. http:
//dx.doi.org/10.1016/j.enganabound.2010.08.015.
[34] Najarzadeh L, Movahedian B, Azhari M. Numerical solution of scalar
wave equation by the modified radial integration boundary element
method. Eng Anal Bound Elem 2019;105:267–78. http://dx.doi.org/10.1016/j.
enganabound.2019.04.027, URL https://www.sciencedirect.com/science/article/
pii/S095579971830701X.
[35] Qu S, Li S, Chen H-R, Qu Z. Radial integration boundary element method
for acoustic eigenvalue problems. Eng Anal Bound Elem 2013;37(7):1043–51.
http://dx.doi.org/10.1016/j.enganabound.2013.03.016.92[36] AL-Jawary MA, Wrobel LC. Numerical solution of the two-dimensional
Helmholtz equation with variable coefficients by the radial integration bound-
ary integral and integro-differential equation methods. Int J Comput Math
2012;89(11):1463–87. http://dx.doi.org/10.1080/00207160.2012.667087.
[37] Campos LS, Loeffler CF, Netto FO, de Jesus dos Santos A. Testing the ac-
complishment of the radial integration method with the direct interpolation
boundary element technique for solving Helmholtz problems. Eng Anal Bound
Elem 2020;110:16–23. http://dx.doi.org/10.1016/j.enganabound.2019.09.022.
[38] Gandomkar M, Dibajian S, Farzin M, Hashemolhoseini S. An integration pro-
cedure for meshless methods using kriging interpolations. Indian J Sci Technol
2013;6(1):3859–67.
[39] Wrobel LC. Boundary element method, volume 1: Applications in thermo-fluids
and acoustics, vol. 1. John Wiley & Sons; 2002.
[40] Katsikadelis JT. The boundary element method for engineers and scientists. 2nd
Ed.. Oxford: Academic Press; 2016.
[41] Slak J, Kosec G. On generation of node distributions for meshless PDE
discretizations. SIAM J Sci Comput 2019;41. http://dx.doi.org/10.1137/
18M1231456.
[42] Gu Y, Wang Q, Lam K. A meshless local kriging method for large deformation
analyses. Comput Methods Appl Mech Engrg 2007;196(9–12):1673–84.
[43] Wang J, Liu G. On the optimal shape parameters of radial basis func-
tions used for 2-D meshless methods. Comput Methods Appl Mech Engrg
2002;191(23–24):2611–30.
[44] Punzi A, Sommariva A, Vianello M. Meshless cubature over the disk us-
ing thin-plate splines. J Comput Appl Math 2008;221(2):430–6. http://dx.doi.
org/10.1016/j.cam.2007.10.023, URL https://www.sciencedirect.com/science/
article/pii/S0377042707005742, Special Issue: Recent Progress in Spline and
Wavelet Approximation.
[45] Najarzadeh L, Movahedian B, Azhari M. Stability analysis of the thin plates
with arbitrary shapes subjected to non-uniform stress fields using boundary
element and radial integration methods. Eng Anal Bound Elem 2018;87:111–21.
http://dx.doi.org/10.1016/j.enganabound.2017.11.010.
[46] Young D, Gu M, Fan C. The time-marching method of fundamental solutions
for wave equations. Eng Anal Bound Elem 2009;33(12):1411–25. http://dx.
doi.org/10.1016/j.enganabound.2009.05.008, URL https://www.sciencedirect.
com/science/article/pii/S0955799709001374, Special Issue on the Method of
Fundamental Solutions in honour of Professor Michael Golberg.
[47] Loeffler C, Mansur W, Barcelos H, Bulcão A. Solving Helmholtz problems with the
boundary element method using direct radial basis function interpolation. Eng
Anal Bound Elem 2015;61:218–25. http://dx.doi.org/10.1016/j.enganabound.
2015.07.013.
[48] Uğurlu B. A dual reciprocity boundary element solution method for the
free vibration analysis of fluid-coupled kirchhoff plates. J Sound Vib
2015;340:190–210. http://dx.doi.org/10.1016/j.jsv.2014.12.011.
[49] Chanthawara K, Kaennakham S, Toutip W. The numerical study and comparison
of radial basis functions in applications of the dual reciprocity boundary element
method to convection-diffusion problems. 1705, 2016, 020029. http://dx.doi.
org/10.1063/1.4940277,
[50] Barbosa J, Loeffler C, Lara L. The direct interpolation boundary element
technique applied to three-dimensional scalar free vibration problems. Eng Anal
Bound Elem 2019;108:295–300. http://dx.doi.org/10.1016/j.enganabound.2019.
09.002.