- Email: [email protected]
Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method
Computers and Mathematics with Applications (
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method Yan Gu a, *, Xiaoqiao He b , Wen Chen c , Chuanzeng Zhang d a
College of Mathematics and Statistics, Qingdao University, Qingdao 266071, PR China Department of Civil and Architectural Engineering, City University of Hong Kong, Hong Kong c Department of Engineering Mechanics, College of Mechanics and Materials, Hohai University, Nanjing 210098, PR China d Department of Civil Engineering, University of Siegen, Paul-Bonatz-Str. 9-11, D-57076 Siegen, Germany b
article
info
Article history: Received 6 November 2016 Received in revised form 15 August 2017 Accepted 19 August 2017 Available online xxxx Keywords: Boundary element method Thin-walled structures Nearly singular integral Three-dimensional problems Anisotropic media
a b s t r a c t In this paper, an advanced boundary element method (BEM) is developed for solving threedimensional (3D) anisotropic heat conduction problems in thin-walled structures. The troublesome nearly singular integrals, which are crucial in the applications of the BEM to thin structures, are calculated efficiently by using a nonlinear coordinate transformation method. For the test problems studied, promising BEM results with only a small number of boundary elements have been obtained when the thickness of the structure is in the orders of micro-scales (10−6 ), which is sufficient for modeling most thin-walled structures as used in, for example, smart materials and thin layered coating systems. The advantages, disadvantages as well as potential applications of the proposed method, as compared with the finite element method (FEM), are also discussed. © 2017 Published by Elsevier Ltd.
1. Introduction The study of boundary value problems for thin structural problems, such as various thin films in electronic devices, sensors and actuators in smart materials, has received considerable attention in recent years [1–3]. This interest is partly related to the extensive use of smart materials and micro-electro-mechanical systems (MEMS) in various engineering applications. The popular finite element method (FEM) has long been a dominant numerical technique in the simulation of various engineering applications [4,5]. However, as illustrated in Refs. [1,6,7], the aspect ratio issues associated with the FEM limit its application to thin-structural problems. To maintain element aspect ratio, a great amount of finite elements should be discretized as the thickness of the structures decreases, which can be arduous, time-consuming and computationally expensive for certain problems. The issue is confirmed indirectly in FEM bibliography, where the relative thickness of 3D thin-walled structures is usually up to the order of 10−3 , see Refs. [8,9], for example. During the past decades, the boundary element method (BEM) has rapidly improved and is nowadays considered as a competing numerical method to the FEM [10–14], thanks to its advantages of semi-analytical nature and boundaryonly discretizations. One of the difficulties in the applications of the BEM to structures with thin shapes is the accurate calculation of the nearly singular integrals [6,15]. Generally speaking, the integrals are nearly singular in the sense that the calculation point is getting towards, but not on, the boundary element. In such case, the integrand, instead of remaining smooth, will develop a sharp peak as the calculation point moves closer to the integration element, thus leading to difficulties
*
Corresponding author. E-mail address: [email protected] (Y. Gu).
http://dx.doi.org/10.1016/j.camwa.2017.08.030 0898-1221/© 2017 Published by Elsevier Ltd.
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
2
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
in accurate evaluation of such integrals [16–18]. Over the past two decades, some considerable effort has been devoted to proposing novel computational techniques to overcome such issues associated with the boundary element analysis. The methods developed so far include, but are not limited to, element subdivision methods [19–21], analytical and semianalytical methods [22–24] and various nonlinear transformations [1,25–31]. Comprehensive reviews of the earlier work in this regard can be found in Refs. [26,32]. Although impressive results have been obtained, a number of drawbacks still remain and mainly include the fact that some techniques failed to provide accurate results when the thickness of the structure is very small [23,26] while others are tailored to some certain kind of integrals [1,33]. In addition, to our knowledge, very few studies on the calculation of nearly singular integrals arising in general anisotropic BEM formulations have been reported in the BEM community [34]. In this study, we focus on a nonlinear coordinate transformation method, named as the ‘sinh’ transformation. The main idea of the method is to smooth out the severe oscillation of the nearly singular integrals using a sinh function before the traditional Gaussian quadrature is applied. It is shown that, as illustrated in Refs. [28,35], this approach is accurate, stable and superior to most of existing methods in terms of overall accuracy and stability. The basis of the sinh transformation was published in the early 2000s by Johnston and Elliott [27,36] and were later essentially improved and extended by many other authors. Gu et al. [37] and Lv et al. [38] extended the method to 2D and 3D nearly singular integrals defined on highorder geometry elements (curved surface elements), respectively. In 2013, the sinh transformed BEM developed in [37] was extended to 2D thin structures and applied to the interfacial stress analysis of multi-coating systems. In a more recent study, Gu et al. [39] extended the transformation to nearly singular integrals arising in general 3D anisotropic boundary element analysis. In this paper, the sinh transformed BEM formulation developed in [39] is extended to 3D thin structural problems associated with anisotropic heat conduction. It is shown that the original 3D anisotropic nearly singular integrals can be transformed into an element-by-element sum of regular integrals, each one expressed in terms of intrinsic (local) coordinates. As a consequence, the actual computation can be performed by using the standard n × n Gaussian quadrature. For the test problems studied, very promising results have been obtained when the thickness of the structures is in the micro-scales. A brief outline of the rest of the paper is as follows. The BEM formulations for 3D anisotropic potential problems are described in Section 2. The nearly singular integrals in three dimensions and their numerical implementation are introduced in Section 3. Followed in Section 4, three benchmark numerical examples with thin shapes are well-studied to demonstrate the accuracy and efficiency of the proposed methodology. Finally, the conclusions and remarks are provided in Section 5. 2. BEM formulations for 3D anisotropic heat conduction problems In this study, we consider anisotropic steady heat conduction applications in the absence of inner heat sources. Hence the function u (x), which stands for the temperature distribution in the computational domain Ω , satisfies the following equation [39] k11
∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) + 2k13 + 2k23 = 0, + k22 + k33 + 2k12 2 2 2 ∂ x1 ∂ x2 ∂ x1 ∂ x3 ∂ x2 ∂ x3 ∂ x1 ∂ x2 ∂ x3
(1)
subject to the following boundary conditions u (x) = u (x) ,
x ∈ ΓD (Dirichlet boundary condition) ,
q (x) = ∇ u (x) · m = q (x) ,
x ∈ ΓN (Neumann boundary condition) ,
(2) (3)
where m = (k11 n1 + k12 n2 + k13 n3 ) i + (k12 n1 + k22 n2 + k23 n3 ) j + (k13 n1 + k23 n2 + k33 n3 ) k ,
(4)
is a vector in the direction of the conormal to the boundary, the boundary of Ω is ∂ Ω ≡ Γ = ΓD ∪ ΓN which we shall assume ( ) to be piecewise smooth, q (x) denotes the heat flux at point x, kij i,j=1,2,3 are thermal conductivity coefficients which are assumed to be symmetric and positive-definite so that the partial differential equation (1) is elliptic, the barred quantities u (x) and q (x) denote the given values of temperature and heat flux specified along the boundary and n = (n1 , n2 , n3 ) indicates the unit outward normal on the boundary points. It is interesting to note that when k11 = k22 = k33 = 1 and k12 = k13 = k23 = 0, the boundary value problem (1)–(3) reduces to the potential problems in isotropic bodies. The fundamental solutions of Eq. (1), which is a basic component of all the subsequent analyses, can be written as [38,39] u∗ (x, y ) =
1
1
√⏐ ⏐ , r (x, y ) 4π ⏐kij ⏐
q∗ (x, y ) = ∇ u∗ (x, y ) · m,
(5)
where r (x, y ) =
√
tij (xi − yi ) xj − yj ,
(
)
(i, j = 1, 2, 3) ,
(6)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
3
the field point y = (y1 , y2 , y3 ). In Eq. (6), the values of is [ the ] distance function between the source point [ x] = (x1 , x2 ,[x3 )]−and 1 tij i,j=1,2,3 can be computed using the relation tij i,j=1,2,3 = kij i,j=1,2,3 , see Refs. [40,41], for example. In the boundary element analysis, the solutions of Eq. (1) can be expressed by using the following boundary integral equations [39,42] c (y ) u (y ) =
∫
[ Γ
q (x) u∗ (x, y ) − u (x) q∗ (x, y ) dΓ (x) ,
y ∈ Γ,
]
(7)
where u (y ) is the temperature field at a typical field point y, c (y ), usually taken to be c (y ) = 0.5 for a smooth boundary, is a coefficient which depends only on the local geometry of Γ . The key problem is the solution of Eq. (1) using the specified boundary conditions (2) and (3) to calculate the remaining boundary values of u (x) and q (x). The quintessence of the BEM is to discretize the boundary into a finite number of segments, not necessarily equal, which are called boundary elements. Two approximations are made over each of these elements. One is about the geometry of the boundary, while the other has to do with the variation of the unknown boundary quantity over the element. To represent the geometry of the boundary, we introduce here, as usual in the BEM, a piecewise eight-node curved quadrilateral element, in terms of local shape functions xj = N k (x) xkj ,
j = 1, 2, 3,
(8)
where k = 1, 2, . . . , 8 and N (x) denote shape functions used to interpolate the geometry, see Ref. [39], of eight nodal points. Similarly, we can express the boundary quantities over each element as k
u (x) = N k (x) uk , k
q (x) = N k (x) qk ,
xkj
are coordinates
(9)
k
where u and q are nodal values of u and q, respectively. After discretizing the boundary into N curved quadrilateral elements, the boundary integral equation (7) may be written as c (y ) u (y ) =
N ∫ ∑ [ j=1
q (x) u∗ (x, y ) − u (x) q∗ (x, y ) dΓ (x) .
Γj
]
(10)
By taking y to all M boundary nodes, say {xm }M m=1 , and taking into account the variation of the boundary quantities on the elements, we can then rewrite Eq. (10) as m m
c u =
N ∑
]
{[∫ u
∗
Γj
j=1
(
x, x
m
)
N (x) dΓ (x) k
qkj
]
[∫
∗
q
− Γj
(
x, x
m
)
N (x) dΓ (x) k
} ukj
.
(11)
Let, for brevity, k Umj =
∫
u∗ x, xm N k (x) dΓ (x) ,
(
Γj
)
k Qmj =
∫
q∗ x, xm N k (x) dΓ (x) ,
(
Γj
)
(12)
and then the boundary integral equation (11) can be rewritten as c m um =
N ∑ (
k k k k Umj qj − Qmj uj ,
)
(13)
j=1 k k where um is M × 1 vector of M nodal values of u, qkj and ukj are M × 1 vectors of nodal values of q and u, Umj and Qmj are M × M matrices of coefficients. It is well-known that the efficiency of the BEM highly relies on the accurate evaluation of various integrals arising in Eq. (12). These integrals become singular or nearly singular when the calculation point locates on or gets close to the boundary surface. Here in this paper, we do not address the singularity issues since many direct and indirect techniques for calculating singular integrals have been developed and used successfully. Readers interested can refer to Refs. [43,44]. In this paper we employ the methods proposed by Guiggiani et al. [45] to calculate the singular integrals arising in Eq. (12). k We can now absorb the coefficients of um on the left-hand side of Eq. (13) with the leading diagonals of the matrix Qmj to k
form a new matrix Q mj , one has N ( ∑
k
k k Umj qj − Q mj ukj
)
= 0.
(14)
j=1
Eq. (14) may be written in matrix form as Q {u} = [U] {q} .
[ ]
(15)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
4
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 1. Quadrangular element splits into (a) four, (b) three and (c) two sub-triangles at the nearly singular point. ‘- - - - -’ denotes the splitting line of the element.
It is noted that for each node only one unknown quantity is admissible (either u or q). Therefore, we should rearrange the unknowns on the basis of the boundary conditions (2) and (3). If the value of u is unknown at a node, then the respective column remains at the left-hand side of Eq. (15), otherwise this column is multiplied by the known values of u, its sign is switched and is shifted to the right-hand side of the equation. Similarly, if the value of q is unknown, then the respective column of [U] is shifted with opposite sign to the left-hand side of the equation. After completing this process, the right-hand side of equation contains only known quantities and, thus, the matrix multiplication results in a single vector. The final system equation can then be expressed as the following linear algebraic system Aα = b,
(16)
where A is the coefficient matrix, α the unknown values (u and/or q) on the boundary, and b the right-hand side vector. Once all the values of α are determined by Eq. (16), the temperature and its derivatives at any point inside the domain or on the boundary can be easily evaluated. 3. Accurate calculation of 3D nearly singular integrals in anisotropic media 3.1. Nearly singular integrals The integrals arising in general 3D boundary element analysis can be expressed in the following form
∫
f
I= Γ
r
dΓ = 2α
1
∫
1
∫
−1
−1
f (s, t ) r 2α (s, t )
|J | dsdt ,
(17)
where Γ is, generally, a curved surface, s and t denote the local intrinsic coordinates, α = 1/2, 3/2 or 5/2 is a constant which indicates the orders of near singularities, |J | is the Jacobian of the transformation from the boundary surface to square [−1, 1] × [−1, 1], and the function f (s, t ) is a polynomial which usually consists of the shape functions used to interpolate the solution and the terms from taking the derivative of the fundamental solutions. To represent the geometry of the boundary, we introduce here, as usual in the BEM applications, an eight-node curved quadrilateral element, in terms of local shape functions xk (s, t ) =
8 ∑
j
Nj (s, t ) xk ,
k = 1, 2, 3,
(18)
j=1
(
j
j
j
)
where Ni (s, t ), as defined in Ref. [39], are the shape functions used to interpolate the geometry, while xj = x1 , x2 , x3 are eight nodal points of the quadrilateral element. Suppose xc is the projection of the near-boundary point y in the local intrinsic coordinates (see Fig. 1(a)), it is convenient to introduce, in each parameter space [−1, 1] × [−1, 1], the following polar coordinate system (ρ, θ) s = s0 + ρ cos (θ) ,
t = t0 + ρ sin (θ) ,
(19)
where (s0 , t0 ) is the local intrinsic coordinates of the projection point x . It is easy to obtain the relation dsdt = ρ dρ dθ . Generally, the integration domain should be split into a set of sub-triangles at the nearly singular point xc . When the nearly singular point is located inside the element, we split the square into four sub-triangles, as shown in Fig. 1(a). Three or two triangles can be obtained when the nearly singular point is located on one side or in a vertex of the square, as shown in Fig. 1(b) and (c). The nearly singular integral (17) can now be expressed as the following element-by-element summation c
I=
n ∫ ∑ m=1
θ2m θ1m
Rm (θ )
∫ 0
f (ρ, θ) r 2α (ρ, θ)
|J | ρ dρ dθ,
(20)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
5
where n is the total number of sub-triangles, θ1m , θ2 and [0, Rm (θ)] indicate the angular and radial span for each subtriangle, respectively (see Fig. 1). It should be noted that the explicit knowledge of s0 and t0 , the local intrinsic coordinates of the projection point xc , plays a key role in the present analysis. Their explicit expressions can be found in Ref. [39]. By using Taylor expansion in the neighborhood of the projection point xc , it is easy to establish the following formulae [38,45]
[
] m
xk − yk = xk − xck + xck − yk ⏐ ⏐ ) ∂ xk ⏐ ( c ( ) ∂ xk ⏐⏐ ⏐ = xk − yk + (s − s0 ) + (t − t0 ) + O ρ 2 ∂ s ⏐ s=s0 ∂ t ⏐ s=s0
t =t 0 ⎡t =t0 ⎤ ⏐ ⏐ ( ) ∂ xk ⏐ ∂ xk ⏐ ⏐ ⏐ = r0 nk (s0 , t0 ) + ρ ⎣ cos (θ) + sin (θ)⎦ + O ρ 2 ∂ s ⏐ s=s0 ∂ t ⏐ s=s0 t =t0
(21)
t =t 0
where k = 1, 2, 3, −1 ≤ s0 , t0 ≤ 1 indicate the local intrinsic coordinates of the projection point, nk denotes the component of the outward normal of the surface point, y = (y1 , y2 , y3 ) and x = (x1 , x2 , x3 ) denote the calculation and integration points, respectively. Let, for brevity,
⎡
⎤ ⏐ ⏐ ∂ xk ⏐ ∂ xk ⏐ ⏐ ⏐ Ak (θ) = ⎣ cos (θ) + sin (θ)⎦ , ∂ s ⏐ s=s0 ∂ t ⏐ s=s0 t =t0
k = 1, 2, or 3,
(22)
t =t 0
and A (θ) =
tij Ai (θ) Aj (θ),
√
(23)
the following equation can be easily obtained r 2 (ρ, θ) = tij (xi − yi ) xj − yj = r02 + ρ 2 A2 (θ) + O ρ 3 .
)
(
( )
(24)
From Eqs. (20) and (24), one has I=
n ∫ ∑ m=1
θ2m
Rm (θ )
∫
θ1m
f (ρ, θ)
( )]α |J | ρ dρ dθ.
(25)
r02 + ρ 2 A2 (θ) + O ρ 3
[
0
It can be seen that, if the value of r0 is very small, the integral (25) will present various orders of near singularities in the radial direction. The key to achieving high accuracy is then to find an algorithm to accurately calculate integral (25) for a small value of r0 , which will be analyzed in Section 3.2. 3.2. The sinh transformation for 3D nearly singular integrals As shown in Eq. (25), the nearly singular integrals under consideration all contain the term r02 + ρ 2 A2 (θ) + O ρ 3 . When faced with evaluating such integral, Johnston and Elliott [46] and Gu and Chen [37,39] suggested the following transformation with the form
( )
ρ = r0 sinh (k1 ξ1 − k2 ) ,
(26)
where k1 and k2 are chosen such that the transformation maps ρ ∈ [0, R (θ)] into ξ1 ∈ [−1, 1] so that the Gaussian quadrature can be applied. For this, the values of k1 and k2 can be calculated by m
⎛ k1 =
1 2
log ⎝
Rm (θ) r0
√
(
+
Rm (θ)
)2
r0
⎞ + 1⎠ and k2 = −k1 .
(27)
The Jacobian of transformation (26) is dρ = r0 k1 cosh (k1 ξ1 − k2 ) dξ1 .
(28)
Substituting (26) into the integral (25) yields I=
∫ n ∑ θ2m − θ1m k1 m=1
2
r02α−2
1
∫
−1
1
−1
f (ξ1 , ξ2 ) sinh (k1 ξ1 − k2 ) cosh (k1 ξ1 − k2 )
]α |J | dξ1 dξ2 ,
1 + sinh2 (k1 ξ1 − k2 ) A2 (ξ2 ) + O (ρ)
[
(29)
where the transformation in the angular direction
θ=
θ2m − θ1m
ξ2 +
θ2m + θ1m
2 2 has been already performed.
,
(30)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
6
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 2. Steady state heat conduction in thin hollow ring torus.
Using the procedure described above, the kernels in the transformed integral (29) are now fully regular regardless of the value of r0 . Thus, the transformed integral can now be evaluated accurately by using the standard Gaussian-quadrature. It is interesting to notice that the Taylor series used in the above analysis is just to perform the asymptotic analysis of the kernel function and to show, in the sense of mathematics, the feasibility of the proposed regularization method. In the actual computation, the proposed transformation can be directly substituted into the nearly singular integral (20), without using the Taylor expansion and the above asymptotic analysis. 4. Numerical results and discussions Three benchmark heat conduction problems in 3D anisotropic media with thin shapes are presented to verify the methodology developed above. The effect of the proposed method as well as the stability, and the convergence with respect to the number of boundary elements are carefully investigated. Unless otherwise specified, the thermal conductivity tensor kij are chosen as kij = 12 + 12 δij , where δij is the Kronecker delta. For the ease of comparison and validation of the numerical results, the analytical solution is taken to be 1
1
3
1
1
x2 x3 + 2. (31) 2 2 2 2 2 The numerical values obtained will be compared with analytical solutions in terms of the relative error defined by u (x) =
x21 +
x22 −
[ Relative Error =
x23 −
M 1 ∑
M
k=1
(
x1 x2 + x1 x3 +
k k Inumerical − Iexact k Iexact
)2 ]1/2
,
(32)
k k where Inumerical and Iexact are the numerical and analytical solutions at the kth calculation point, respectively.
4.1. Test problem 1: heat conduction in a thin hollow ring torus Firstly, we consider the steady state heat conduction in a thin hollow ring torus, as shown in Fig. 2. In mathematics, a torus can be defined parametrically by x (φ, θ) = (f (φ) cos θ, f (φ) sin θ, r sin φ) ,
(33)
where f (φ) = 4 + cos (φ), and φ and θ vary in the interval [0, 2π ). The inner radius of the ring torus is a = 1 which is constant in this study, while the outer radius b changes from 1.1 to 1.000000001. It can be seen that the thickness of the structure (δ ) changes in the ranges of 10−1 ≤ δ ≤ 10−9 . In this 3D model, Dirichlet boundary conditions are specified on the outer surface of the ring torus, while the Neumann boundary conditions are specified at the inner surface of the ring torus. The actual values of the Dirichlet and Neumann boundary conditions can be easily obtained using the corresponding analytical solutions given by Eq. (31). Here, a total number of N = 300 evenly distributed boundary elements, in terms of the angle parameterization, are chosen on the whole surfaces, regardless of the thickness of the structure. Figs. 3 and 4 illustrate the relative error curves of u (x) and ∂ u (x) /∂ x1 calculated at the interior point (4 + (a + b) /2, 0, 0), as functions of different thickness-to-length ratio. Both the untransformed BEM and the sinh-transformed BEM proposed in this paper are employed for the purpose of comparisons. As can be seen in Figs. 3 and 4, when the thickness of the structure is not very small, both the BEM schemes are effective and can yield accurate numerical results. However, the untransformed BEM performs less satisfactorily as the thickness decreases, i.e., after the thickness is less than 1E−1. In contrast, the results obtained using the proposed method are excellently consistent with the exact solutions, even when the thickness is in the order of 10−9 . In addition, the number of BEM elements does not change across the entire range of thickness, and therefore, the solution time and memory requirements are quite modest. Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
7
Fig. 3. Relative error curves of temperatures as functions of different thicknesses of the structure.
Fig. 4. Relative error curves of fluxes as functions of different thicknesses of the structure.
4.2. Test problem 2: heat conduction in a thin hollow peanut-shaped domain Next, we consider steady state heat conduction in a thin hollow peanut-shaped domain, as shown in Fig. 5. The geometry of the outer and inner surfaces of the computational domain is defined parametrically as x (φ, θ) = ((f (φ) + r ) sin φ cos θ, (f (φ) + r ) sin φ sin θ, (f (φ) + r ) cos φ) ,
(34)
x (φ, θ) = (f (φ) sin φ cos θ, f (φ) sin φ sin θ, f (φ) cos φ) ,
(35)
√
√
respectively, where φ ∈ [0, π) , θ ∈ [0, 2π), f (φ) = cos (2φ) + 1.1 − sin2 (2φ) and r = 10−6 denoting the thickness of the structure. Here, Neumann boundary conditions are specified on the outer surface of the structure and Dirichlet boundary conditions on the inner surface of the structure. To investigate the convergence of method, Figs. 6 and 7 display the relative error surfaces of u (x) and ∂ u (x) /∂ x1 using N = 400 and N = 800 evenly distributed boundary elements (in terms of the angle parameterization), respectively. In the case of N = 400 boundary elements, there are 25 elements discretized along the θ -direction (along the peanut-shaped Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
8
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 5. Steady state heat conduction in a peanut-shaped domain.
a
b
10-4
0.012
8
0.01 Relative errors
Relative errors
6 0.008 0.006 0.004 0.002 30
0 30
4
2
0 30
20 20
10
10 0
30
20 20 10
0
10 0
0
Fig. 6. Relative error surfaces of temperatures, where the number of BEM elements is 400 (a) and 800 (b), respectively.
domain) and 16 elements along the φ -direction (around the peanut-shaped domain). In the case of N = 800 boundary elements, there are 40 elements discretized along the θ -direction and 20 elements along the φ -direction. The error surfaces were yielded at M = 30 × 30 calculation points uniformly-spaced over the surface x (φ, θ) = ((f (φ) + r /2) sin φ cos θ, (f (φ) + r /2) sin φ sin θ, (f (φ) + r /2) cos φ) .
(36) −7
It is noted that the distance between the calculation points to the boundary is in the orders of 10 . As can be seen from Figs. 6 and 7, the present method is stable, accurate, and rapidly convergent as the number of boundary elements increases. It can be also observed that the BEM results with only 400 boundary elements are quite accurate for this 3D case, so that the size of resulting system of linear algebraic equations is quite small. 4.3. Test problem 3: heat conduction in a thin hollow hyperboloid domain Finally, we consider the steady state heat conduction in a thin hollow hyperboloid domain, as shown in Fig. 8. The geometry of the outer surface is defined as x (θ , t ) =
(√
1 + t 2 cos θ,
√
)
1 + t 2 sin θ, t ,
(37)
where the two parameters θ and t vary between 0 ∼ 2π and −2 ∼ 2, respectively. The inner boundary is taken to be a curve similar to the outer surface. Here the distance between the inner and outer surfaces changes from 10−1 to 10−8 . In this problem, Dirichlet boundary conditions are prescribed on the outer surface, while Neumann boundary conditions are specified at the remaining surface. To show the numerical accuracy, M = 100 interior points uniformly-spaced over the surface between the inner and outer boundaries are chosen. Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
a
)
–
9
b 10-3
0.015
1.5
0.01 Relative errors
Relative errors
2
0.005
0 30
1 0.5 0 30
30
20 20 10
10 0
0
30
20 20 10
10 0
0
Fig. 7. Relative error surfaces of fluxes, where the number of BEM elements is 400 (a) and 800 (b), respectively.
Fig. 8. Steady state heat conduction in a thin hollow hyperboloid domain.
Here, both the BEM and the FEM are employed for the purpose of comparison. For the BEM model, the total number of N = 600 boundary elements is chosen, regardless of the thickness of the structure. The FEM discretization varies according to the structure thickness. In general, more finite elements are required as the thickness of the structure decreases, in order to maintain reasonable element aspect ratio. Fig. 9 illustrates the relative error curves of the u (x) using both the BEM and the FEM, as functions of different values of thicknesses. The number of BEM and FEM elements employed in the calculation is provided in Fig. 10. As can be seen from Fig. 9, the BEM results are excellently consistent with the analytical solutions, even when the thickness of the structure is as small as 10−8 . And the number of BEM elements does not change across the entire range of thickness. The FEM solution, however, shows a very different behavior. When the thickness of the structure is less than 0.01, the problem is easily solved using the FEM. However, the FEM solution requires significantly more effort when the thickness of the structure further decreases due to aspect ratio limitations. Hence we can conclude that, in comparison with the FEM, the proposed BEM could be considered a competitive alternative for solving thin structural problems. 5. Conclusions This paper investigates the applicability of the BEM for three-dimensional anisotropic heat conduction problems on thin domains. It is shown that the BEM can deal with ultra-thin structures very efficiently and accurately, as long as the Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
10
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 9. Temperature relative errors by using the BEM and FEM for different thicknesses of the structure.
Fig. 10. Numbers of elements of the BEM and FEM required for different thicknesses of the structure.
nearly singular integrals existing in the BEM formulations are handled correctly. It can be seen that the BEM does not suffer from thickness or aspect ratio issues associated with FEM. Therefore, the number of BEM elements can be held constant as described in examples. The method proposed in this paper can be applied to other 3D anisotropic problems on thin domains in the analogous way. It must be pointed out that the proposed BEM approach to thin-bodies should be considered as a complement to the FEM. For structures that can be identified clearly as shells or bulky solids, the FEM should be used because of the efficiency FEM delivers. For thin structures where the solid models are difficult to obtain, the BEM can be used as an alternative, especially when high accuracy is desired, as in benchmarks. Acknowledgments The work described in this paper was supported by the National Basic Research Program of China (973 Project No. 2010CB832702), the National Natural Science Foundation of China (Nos. 11402075, 71571108), Projects of International (Regional) Cooperation and Exchanges of NSFC (Nos. 71611530712, 61661136002), the China Postdoctoral Science Foundation (Nos. 2015M570572, 2016T90608), and the Qingdao Postdoctoral Application Research Project (No. 15-9-1-49-jc). Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
11
References [1] J.F. Luo, Y.J. Liu, E.J. Berger, Analysis of two-dimensional thin structures (from micro- to nano-scales) using the boundary element method, Comput. Mech. 22 (1998) 404–412. [2] B.P. Lamichhane, S.G. Roberts, L. Stals, A mixed finite element discretisation of thin plate splines based on biorthogonal systems, J. Sci. Comput. 67 (2016) 20–42. [3] L. Marin, A meshless method for solving the Cauchy problem in three-dimensional elastostatics, Comput. Math. Appl. 50 (2005) 73–92. [4] A.H.D. Cheng, D.T. Cheng, Heritage and early history of the boundary element method, Eng. Anal. Bound. Elem. 29 (2005) 268–302. [5] M.G. Knight, L.C. Wrobel, J.L. Henshall, Micromechanical response of fibre-reinforced materials using the boundary element technique, Compos. Struct. 62 (2003) 341–352. [6] Y. Gu, W. Chen, C. Zhang, Stress analysis for thin multilayered coating systems using a sinh transformed boundary element method, Internat. J. Solids Struct. 50 (2013) 3460–3471. [7] R.E. Caflisch, Growth, Structure and pattern formation for thin films, J. Sci. Comput. 37 (2008) 3–17. [8] B.A. Szabó, G.J. Sahrmann, Hierarchic plate and shell models based on p-extension, Int. J. Numer. Meth. Eng. 26 (1988) 1855–1881. [9] G. Zboinski, Adaptive hpq finite element methods for the analysis of 3D-based models of complex structures. Part 1. Hierarchical modeling and approximations, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2913–2940. [10] M. Guiggiani, Computing principal-value integrals in 3-D BEM for time-harmonic elastodynamics—A direct approach, Commun. Appl. Numer. Methods 8 (1992) 141–149. [11] Y. Gu, W. Chen, C.Z. Zhang, Singular boundary method for solving plane strain elastostatic problems, Internat. J. Solids Struct. 48 (2011) 2549–2556. [12] J.T. Chen, Y.L. Chang, S.K. Kao, J. Jian, Revisit of the indirect boundary element method: Necessary and sufficient formulation, J. Sci. Comput. 65 (2015) 467–485. [13] L. Demkowicz, Asymptotic convergence in finite and boundary element methods: part 1: theoretical results, Comput. Math. Appl. 27 (1994) 69–84. [14] L. Demkowicz, Asymptotic convergence in finite and boundary element methods: part 2: The LBB constant for rigid and elastic scattering problems, Comput. Math. Appl. 28 (1994) 93–109. [15] V. Sladek, J. Sladek, M. Tanaka, Nonsingular BEM formulations for thin-walled structures and elastostatic crack problems, Acta Mech. 99 (1993) 173–190. [16] Y.J. Liu, Analysis of shell-like structures by the Boundary Element Method based on 3-D elasticity: formulation and verification, Int. J. Numer. Meth. Eng. 41 (1998) 541–558. [17] A. Nagarajan, S. Mukherjee, A mapping method for numerical evaluation of two-dimensional integrals with 1/r singularity, Comput. Mech. 12 (1993) 19–26. [18] M. Guiggiani, The evaluation of Cauchy principal value integrals in the boundary element method—a review, Math. Comput. Modelling 15 (1991) 175–184. [19] J.C. Lachat, J.O. Watson, Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics, Int. J. Numer. Meth. Eng. 10 (1976) 991–1005. [20] X.W. Gao, K. Yang, J. Wang, An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals, Eng. Anal. Bound. Elem. 32 (2008) 692–696. [21] X.W. Gao, Numerical evaluation of two-dimensional singular boundary integrals—Theory and Fortran code, J. Comput. Appl. Math. 188 (2006) 44–64. [22] G.S. Padhi, R.A. Shenoi, S.S.J. Moy, M.A. McCarthy, Analytic integration of kernel shape function product integrals in the boundary element method, Comput. Struct. 79 (2001) 1325–1333. [23] H.L. Zhou, Z.R. Niu, C.Z. Cheng, Z.W. Guan, Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for anisotropic potential problems, Comput. Struct. 86 (2008) 1656–1671. [24] Z.R. Niu, C.Z. Cheng, H.L. Zhou, Z.J. Hu, Analytic formulations for calculating nearly singular integrals in two-dimensional BEM, Eng. Anal. Bound. Elem. 31 (2007) 949–964. [25] J.C.F. Telles, A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals, Int. J. Numer. Meth. Eng. 24 (1987) 959–973. [26] Q. Huang, T.A. Cruse, Some notes on singular integral techniques in boundary element analysis, Int. J. Numer. Meth. Eng. 36 (1993) 2643–2659. [27] P.R. Johnston, D. Elliott, A sinh transformation for evaluating nearly singular boundary element integrals, Int. J. Numer. Meth. Eng. 62 (2005) 564–578. [28] P.R. Johnston, B.M. Johnston, D. Elliott, Using the iterated sinh transformation to evaluate two dimensional nearly singular boundary element integrals, Eng. Anal. Bound. Elem. 37 (2013) 708–718. [29] V. Sladek, J. Sladek, M. Tanaka, Optimal transformations of the integration variables in computation of singular integrals in BEM, Int. J. Numer. Meth. Eng. 47 (2000) 1263–1283. [30] Y.M. Zhang, Y. Gu, J.T. Chen, Boundary element analysis of the thermal behaviour in thin-coated cutting tools, Eng. Anal. Bound. Elem. 34 (2010) 775–784. [31] G. Xie, J. Zhang, X. Qin, G. Li, New variable transformations for evaluating nearly singular integrals in 2D boundary element method, Eng. Anal. Bound. Elem. 35 (2011) 811–817. [32] Y. Gu, Q. Hua, W. Chen, C. Zhang, Numerical evaluation of nearly hyper-singular integrals in the boundary element analysis, Comput. Struct. 167 (2016) 15–23. [33] H. Ma, N. Kamiya, A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two- and threedimensional elasticity, Comput. Mech. 29 (2002) 277–288. [34] Y. Gu, W. Chen, B. Zhang, W. Qu, Two general algorithms for nearly singular integrals in two dimensional anisotropic boundary element method, Comput. Mech. 53 (2014) 1223–1234. [35] B.M. Johnston, P.R. Johnston, D. Elliott, A new method for the numerical evaluation of nearly singular integrals on triangular elements in the 3D boundary element method, J. Comput. Appl. Math. 245 (2013) 148–161. [36] P.R. Johnston, Application of sigmoidal transformations to weakly singular and near-singular boundary element integrals, Int. J. Numer. Meth. Eng. 45 (1999) 1333–1348. [37] Y. Gu, W. Chen, C. Zhang, The sinh transformation for evaluating nearly singular boundary element integrals over high-order geometry elements, Eng. Anal. Bound. Elem. 37 (2013) 301–308. [38] J. Lv, Y. Miao, H. Zhu, The distance sinh transformation for the numerical evaluation of nearly singular integrals over curved surface elements, Comput. Mech. 53 (2014) 359–367. [39] Y. Gu, H. Gao, W. Chen, C. Zhang, A general algorithm for evaluating nearly singular integrals in anisotropic three-dimensional boundary element analysis, Comput. Methods Appl. Mech. Engrg. 308 (2016) 483–498. [40] Y. Gu, W. Chen, X.Q. He, Singular boundary method for steady-state heat conduction in three dimensional general anisotropic media, Int. J. Heat. Mass. Transfer 55 (2012) 4837–4848.
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
12
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
[41] Y.M. Zhang, Z.Y. Liu, J.T. Chen, Y. Gu, A novel boundary element approach for solving the anisotropic potential problems, Eng. Anal. Bound. Elem. 35 (2011) 1181–1189. [42] P.K. Banerjee, The Boundary Element Methods in Engineering, McGRAW-HILL Book Company Europe, 1994. [43] M. Guiggiani, P. Casalini, Direct computation of Cauchy principal value integrals in advanced boundary elements, Int. J. Numer. Meth. Eng. 24 (1987) 1711–1720. [44] V. Sladek, J. Sladek, Singular integrals and boundary elements, Comput. Methods Appl. Mech. Engrg. 157 (1998) 251–266. [45] M. Guiggiani, A. Gigante, A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method, J. Appl. Mech. 57 (1990) 906–915. [46] B.M. Johnston, P.R. Johnston, D. Elliott, A sinh transformation for evaluating two-dimensional nearly singular boundary element integrals, Int. J. Numer. Meth. Eng 69 (2007) 1460–1479.
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method Yan Gu a, *, Xiaoqiao He b , Wen Chen c , Chuanzeng Zhang d a
College of Mathematics and Statistics, Qingdao University, Qingdao 266071, PR China Department of Civil and Architectural Engineering, City University of Hong Kong, Hong Kong c Department of Engineering Mechanics, College of Mechanics and Materials, Hohai University, Nanjing 210098, PR China d Department of Civil Engineering, University of Siegen, Paul-Bonatz-Str. 9-11, D-57076 Siegen, Germany b
article
info
Article history: Received 6 November 2016 Received in revised form 15 August 2017 Accepted 19 August 2017 Available online xxxx Keywords: Boundary element method Thin-walled structures Nearly singular integral Three-dimensional problems Anisotropic media
a b s t r a c t In this paper, an advanced boundary element method (BEM) is developed for solving threedimensional (3D) anisotropic heat conduction problems in thin-walled structures. The troublesome nearly singular integrals, which are crucial in the applications of the BEM to thin structures, are calculated efficiently by using a nonlinear coordinate transformation method. For the test problems studied, promising BEM results with only a small number of boundary elements have been obtained when the thickness of the structure is in the orders of micro-scales (10−6 ), which is sufficient for modeling most thin-walled structures as used in, for example, smart materials and thin layered coating systems. The advantages, disadvantages as well as potential applications of the proposed method, as compared with the finite element method (FEM), are also discussed. © 2017 Published by Elsevier Ltd.
1. Introduction The study of boundary value problems for thin structural problems, such as various thin films in electronic devices, sensors and actuators in smart materials, has received considerable attention in recent years [1–3]. This interest is partly related to the extensive use of smart materials and micro-electro-mechanical systems (MEMS) in various engineering applications. The popular finite element method (FEM) has long been a dominant numerical technique in the simulation of various engineering applications [4,5]. However, as illustrated in Refs. [1,6,7], the aspect ratio issues associated with the FEM limit its application to thin-structural problems. To maintain element aspect ratio, a great amount of finite elements should be discretized as the thickness of the structures decreases, which can be arduous, time-consuming and computationally expensive for certain problems. The issue is confirmed indirectly in FEM bibliography, where the relative thickness of 3D thin-walled structures is usually up to the order of 10−3 , see Refs. [8,9], for example. During the past decades, the boundary element method (BEM) has rapidly improved and is nowadays considered as a competing numerical method to the FEM [10–14], thanks to its advantages of semi-analytical nature and boundaryonly discretizations. One of the difficulties in the applications of the BEM to structures with thin shapes is the accurate calculation of the nearly singular integrals [6,15]. Generally speaking, the integrals are nearly singular in the sense that the calculation point is getting towards, but not on, the boundary element. In such case, the integrand, instead of remaining smooth, will develop a sharp peak as the calculation point moves closer to the integration element, thus leading to difficulties
*
Corresponding author. E-mail address: [email protected] (Y. Gu).
http://dx.doi.org/10.1016/j.camwa.2017.08.030 0898-1221/© 2017 Published by Elsevier Ltd.
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
2
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
in accurate evaluation of such integrals [16–18]. Over the past two decades, some considerable effort has been devoted to proposing novel computational techniques to overcome such issues associated with the boundary element analysis. The methods developed so far include, but are not limited to, element subdivision methods [19–21], analytical and semianalytical methods [22–24] and various nonlinear transformations [1,25–31]. Comprehensive reviews of the earlier work in this regard can be found in Refs. [26,32]. Although impressive results have been obtained, a number of drawbacks still remain and mainly include the fact that some techniques failed to provide accurate results when the thickness of the structure is very small [23,26] while others are tailored to some certain kind of integrals [1,33]. In addition, to our knowledge, very few studies on the calculation of nearly singular integrals arising in general anisotropic BEM formulations have been reported in the BEM community [34]. In this study, we focus on a nonlinear coordinate transformation method, named as the ‘sinh’ transformation. The main idea of the method is to smooth out the severe oscillation of the nearly singular integrals using a sinh function before the traditional Gaussian quadrature is applied. It is shown that, as illustrated in Refs. [28,35], this approach is accurate, stable and superior to most of existing methods in terms of overall accuracy and stability. The basis of the sinh transformation was published in the early 2000s by Johnston and Elliott [27,36] and were later essentially improved and extended by many other authors. Gu et al. [37] and Lv et al. [38] extended the method to 2D and 3D nearly singular integrals defined on highorder geometry elements (curved surface elements), respectively. In 2013, the sinh transformed BEM developed in [37] was extended to 2D thin structures and applied to the interfacial stress analysis of multi-coating systems. In a more recent study, Gu et al. [39] extended the transformation to nearly singular integrals arising in general 3D anisotropic boundary element analysis. In this paper, the sinh transformed BEM formulation developed in [39] is extended to 3D thin structural problems associated with anisotropic heat conduction. It is shown that the original 3D anisotropic nearly singular integrals can be transformed into an element-by-element sum of regular integrals, each one expressed in terms of intrinsic (local) coordinates. As a consequence, the actual computation can be performed by using the standard n × n Gaussian quadrature. For the test problems studied, very promising results have been obtained when the thickness of the structures is in the micro-scales. A brief outline of the rest of the paper is as follows. The BEM formulations for 3D anisotropic potential problems are described in Section 2. The nearly singular integrals in three dimensions and their numerical implementation are introduced in Section 3. Followed in Section 4, three benchmark numerical examples with thin shapes are well-studied to demonstrate the accuracy and efficiency of the proposed methodology. Finally, the conclusions and remarks are provided in Section 5. 2. BEM formulations for 3D anisotropic heat conduction problems In this study, we consider anisotropic steady heat conduction applications in the absence of inner heat sources. Hence the function u (x), which stands for the temperature distribution in the computational domain Ω , satisfies the following equation [39] k11
∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) ∂ 2 u (x) + 2k13 + 2k23 = 0, + k22 + k33 + 2k12 2 2 2 ∂ x1 ∂ x2 ∂ x1 ∂ x3 ∂ x2 ∂ x3 ∂ x1 ∂ x2 ∂ x3
(1)
subject to the following boundary conditions u (x) = u (x) ,
x ∈ ΓD (Dirichlet boundary condition) ,
q (x) = ∇ u (x) · m = q (x) ,
x ∈ ΓN (Neumann boundary condition) ,
(2) (3)
where m = (k11 n1 + k12 n2 + k13 n3 ) i + (k12 n1 + k22 n2 + k23 n3 ) j + (k13 n1 + k23 n2 + k33 n3 ) k ,
(4)
is a vector in the direction of the conormal to the boundary, the boundary of Ω is ∂ Ω ≡ Γ = ΓD ∪ ΓN which we shall assume ( ) to be piecewise smooth, q (x) denotes the heat flux at point x, kij i,j=1,2,3 are thermal conductivity coefficients which are assumed to be symmetric and positive-definite so that the partial differential equation (1) is elliptic, the barred quantities u (x) and q (x) denote the given values of temperature and heat flux specified along the boundary and n = (n1 , n2 , n3 ) indicates the unit outward normal on the boundary points. It is interesting to note that when k11 = k22 = k33 = 1 and k12 = k13 = k23 = 0, the boundary value problem (1)–(3) reduces to the potential problems in isotropic bodies. The fundamental solutions of Eq. (1), which is a basic component of all the subsequent analyses, can be written as [38,39] u∗ (x, y ) =
1
1
√⏐ ⏐ , r (x, y ) 4π ⏐kij ⏐
q∗ (x, y ) = ∇ u∗ (x, y ) · m,
(5)
where r (x, y ) =
√
tij (xi − yi ) xj − yj ,
(
)
(i, j = 1, 2, 3) ,
(6)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
3
the field point y = (y1 , y2 , y3 ). In Eq. (6), the values of is [ the ] distance function between the source point [ x] = (x1 , x2 ,[x3 )]−and 1 tij i,j=1,2,3 can be computed using the relation tij i,j=1,2,3 = kij i,j=1,2,3 , see Refs. [40,41], for example. In the boundary element analysis, the solutions of Eq. (1) can be expressed by using the following boundary integral equations [39,42] c (y ) u (y ) =
∫
[ Γ
q (x) u∗ (x, y ) − u (x) q∗ (x, y ) dΓ (x) ,
y ∈ Γ,
]
(7)
where u (y ) is the temperature field at a typical field point y, c (y ), usually taken to be c (y ) = 0.5 for a smooth boundary, is a coefficient which depends only on the local geometry of Γ . The key problem is the solution of Eq. (1) using the specified boundary conditions (2) and (3) to calculate the remaining boundary values of u (x) and q (x). The quintessence of the BEM is to discretize the boundary into a finite number of segments, not necessarily equal, which are called boundary elements. Two approximations are made over each of these elements. One is about the geometry of the boundary, while the other has to do with the variation of the unknown boundary quantity over the element. To represent the geometry of the boundary, we introduce here, as usual in the BEM, a piecewise eight-node curved quadrilateral element, in terms of local shape functions xj = N k (x) xkj ,
j = 1, 2, 3,
(8)
where k = 1, 2, . . . , 8 and N (x) denote shape functions used to interpolate the geometry, see Ref. [39], of eight nodal points. Similarly, we can express the boundary quantities over each element as k
u (x) = N k (x) uk , k
q (x) = N k (x) qk ,
xkj
are coordinates
(9)
k
where u and q are nodal values of u and q, respectively. After discretizing the boundary into N curved quadrilateral elements, the boundary integral equation (7) may be written as c (y ) u (y ) =
N ∫ ∑ [ j=1
q (x) u∗ (x, y ) − u (x) q∗ (x, y ) dΓ (x) .
Γj
]
(10)
By taking y to all M boundary nodes, say {xm }M m=1 , and taking into account the variation of the boundary quantities on the elements, we can then rewrite Eq. (10) as m m
c u =
N ∑
]
{[∫ u
∗
Γj
j=1
(
x, x
m
)
N (x) dΓ (x) k
qkj
]
[∫
∗
q
− Γj
(
x, x
m
)
N (x) dΓ (x) k
} ukj
.
(11)
Let, for brevity, k Umj =
∫
u∗ x, xm N k (x) dΓ (x) ,
(
Γj
)
k Qmj =
∫
q∗ x, xm N k (x) dΓ (x) ,
(
Γj
)
(12)
and then the boundary integral equation (11) can be rewritten as c m um =
N ∑ (
k k k k Umj qj − Qmj uj ,
)
(13)
j=1 k k where um is M × 1 vector of M nodal values of u, qkj and ukj are M × 1 vectors of nodal values of q and u, Umj and Qmj are M × M matrices of coefficients. It is well-known that the efficiency of the BEM highly relies on the accurate evaluation of various integrals arising in Eq. (12). These integrals become singular or nearly singular when the calculation point locates on or gets close to the boundary surface. Here in this paper, we do not address the singularity issues since many direct and indirect techniques for calculating singular integrals have been developed and used successfully. Readers interested can refer to Refs. [43,44]. In this paper we employ the methods proposed by Guiggiani et al. [45] to calculate the singular integrals arising in Eq. (12). k We can now absorb the coefficients of um on the left-hand side of Eq. (13) with the leading diagonals of the matrix Qmj to k
form a new matrix Q mj , one has N ( ∑
k
k k Umj qj − Q mj ukj
)
= 0.
(14)
j=1
Eq. (14) may be written in matrix form as Q {u} = [U] {q} .
[ ]
(15)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
4
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 1. Quadrangular element splits into (a) four, (b) three and (c) two sub-triangles at the nearly singular point. ‘- - - - -’ denotes the splitting line of the element.
It is noted that for each node only one unknown quantity is admissible (either u or q). Therefore, we should rearrange the unknowns on the basis of the boundary conditions (2) and (3). If the value of u is unknown at a node, then the respective column remains at the left-hand side of Eq. (15), otherwise this column is multiplied by the known values of u, its sign is switched and is shifted to the right-hand side of the equation. Similarly, if the value of q is unknown, then the respective column of [U] is shifted with opposite sign to the left-hand side of the equation. After completing this process, the right-hand side of equation contains only known quantities and, thus, the matrix multiplication results in a single vector. The final system equation can then be expressed as the following linear algebraic system Aα = b,
(16)
where A is the coefficient matrix, α the unknown values (u and/or q) on the boundary, and b the right-hand side vector. Once all the values of α are determined by Eq. (16), the temperature and its derivatives at any point inside the domain or on the boundary can be easily evaluated. 3. Accurate calculation of 3D nearly singular integrals in anisotropic media 3.1. Nearly singular integrals The integrals arising in general 3D boundary element analysis can be expressed in the following form
∫
f
I= Γ
r
dΓ = 2α
1
∫
1
∫
−1
−1
f (s, t ) r 2α (s, t )
|J | dsdt ,
(17)
where Γ is, generally, a curved surface, s and t denote the local intrinsic coordinates, α = 1/2, 3/2 or 5/2 is a constant which indicates the orders of near singularities, |J | is the Jacobian of the transformation from the boundary surface to square [−1, 1] × [−1, 1], and the function f (s, t ) is a polynomial which usually consists of the shape functions used to interpolate the solution and the terms from taking the derivative of the fundamental solutions. To represent the geometry of the boundary, we introduce here, as usual in the BEM applications, an eight-node curved quadrilateral element, in terms of local shape functions xk (s, t ) =
8 ∑
j
Nj (s, t ) xk ,
k = 1, 2, 3,
(18)
j=1
(
j
j
j
)
where Ni (s, t ), as defined in Ref. [39], are the shape functions used to interpolate the geometry, while xj = x1 , x2 , x3 are eight nodal points of the quadrilateral element. Suppose xc is the projection of the near-boundary point y in the local intrinsic coordinates (see Fig. 1(a)), it is convenient to introduce, in each parameter space [−1, 1] × [−1, 1], the following polar coordinate system (ρ, θ) s = s0 + ρ cos (θ) ,
t = t0 + ρ sin (θ) ,
(19)
where (s0 , t0 ) is the local intrinsic coordinates of the projection point x . It is easy to obtain the relation dsdt = ρ dρ dθ . Generally, the integration domain should be split into a set of sub-triangles at the nearly singular point xc . When the nearly singular point is located inside the element, we split the square into four sub-triangles, as shown in Fig. 1(a). Three or two triangles can be obtained when the nearly singular point is located on one side or in a vertex of the square, as shown in Fig. 1(b) and (c). The nearly singular integral (17) can now be expressed as the following element-by-element summation c
I=
n ∫ ∑ m=1
θ2m θ1m
Rm (θ )
∫ 0
f (ρ, θ) r 2α (ρ, θ)
|J | ρ dρ dθ,
(20)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
5
where n is the total number of sub-triangles, θ1m , θ2 and [0, Rm (θ)] indicate the angular and radial span for each subtriangle, respectively (see Fig. 1). It should be noted that the explicit knowledge of s0 and t0 , the local intrinsic coordinates of the projection point xc , plays a key role in the present analysis. Their explicit expressions can be found in Ref. [39]. By using Taylor expansion in the neighborhood of the projection point xc , it is easy to establish the following formulae [38,45]
[
] m
xk − yk = xk − xck + xck − yk ⏐ ⏐ ) ∂ xk ⏐ ( c ( ) ∂ xk ⏐⏐ ⏐ = xk − yk + (s − s0 ) + (t − t0 ) + O ρ 2 ∂ s ⏐ s=s0 ∂ t ⏐ s=s0
t =t 0 ⎡t =t0 ⎤ ⏐ ⏐ ( ) ∂ xk ⏐ ∂ xk ⏐ ⏐ ⏐ = r0 nk (s0 , t0 ) + ρ ⎣ cos (θ) + sin (θ)⎦ + O ρ 2 ∂ s ⏐ s=s0 ∂ t ⏐ s=s0 t =t0
(21)
t =t 0
where k = 1, 2, 3, −1 ≤ s0 , t0 ≤ 1 indicate the local intrinsic coordinates of the projection point, nk denotes the component of the outward normal of the surface point, y = (y1 , y2 , y3 ) and x = (x1 , x2 , x3 ) denote the calculation and integration points, respectively. Let, for brevity,
⎡
⎤ ⏐ ⏐ ∂ xk ⏐ ∂ xk ⏐ ⏐ ⏐ Ak (θ) = ⎣ cos (θ) + sin (θ)⎦ , ∂ s ⏐ s=s0 ∂ t ⏐ s=s0 t =t0
k = 1, 2, or 3,
(22)
t =t 0
and A (θ) =
tij Ai (θ) Aj (θ),
√
(23)
the following equation can be easily obtained r 2 (ρ, θ) = tij (xi − yi ) xj − yj = r02 + ρ 2 A2 (θ) + O ρ 3 .
)
(
( )
(24)
From Eqs. (20) and (24), one has I=
n ∫ ∑ m=1
θ2m
Rm (θ )
∫
θ1m
f (ρ, θ)
( )]α |J | ρ dρ dθ.
(25)
r02 + ρ 2 A2 (θ) + O ρ 3
[
0
It can be seen that, if the value of r0 is very small, the integral (25) will present various orders of near singularities in the radial direction. The key to achieving high accuracy is then to find an algorithm to accurately calculate integral (25) for a small value of r0 , which will be analyzed in Section 3.2. 3.2. The sinh transformation for 3D nearly singular integrals As shown in Eq. (25), the nearly singular integrals under consideration all contain the term r02 + ρ 2 A2 (θ) + O ρ 3 . When faced with evaluating such integral, Johnston and Elliott [46] and Gu and Chen [37,39] suggested the following transformation with the form
( )
ρ = r0 sinh (k1 ξ1 − k2 ) ,
(26)
where k1 and k2 are chosen such that the transformation maps ρ ∈ [0, R (θ)] into ξ1 ∈ [−1, 1] so that the Gaussian quadrature can be applied. For this, the values of k1 and k2 can be calculated by m
⎛ k1 =
1 2
log ⎝
Rm (θ) r0
√
(
+
Rm (θ)
)2
r0
⎞ + 1⎠ and k2 = −k1 .
(27)
The Jacobian of transformation (26) is dρ = r0 k1 cosh (k1 ξ1 − k2 ) dξ1 .
(28)
Substituting (26) into the integral (25) yields I=
∫ n ∑ θ2m − θ1m k1 m=1
2
r02α−2
1
∫
−1
1
−1
f (ξ1 , ξ2 ) sinh (k1 ξ1 − k2 ) cosh (k1 ξ1 − k2 )
]α |J | dξ1 dξ2 ,
1 + sinh2 (k1 ξ1 − k2 ) A2 (ξ2 ) + O (ρ)
[
(29)
where the transformation in the angular direction
θ=
θ2m − θ1m
ξ2 +
θ2m + θ1m
2 2 has been already performed.
,
(30)
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
6
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 2. Steady state heat conduction in thin hollow ring torus.
Using the procedure described above, the kernels in the transformed integral (29) are now fully regular regardless of the value of r0 . Thus, the transformed integral can now be evaluated accurately by using the standard Gaussian-quadrature. It is interesting to notice that the Taylor series used in the above analysis is just to perform the asymptotic analysis of the kernel function and to show, in the sense of mathematics, the feasibility of the proposed regularization method. In the actual computation, the proposed transformation can be directly substituted into the nearly singular integral (20), without using the Taylor expansion and the above asymptotic analysis. 4. Numerical results and discussions Three benchmark heat conduction problems in 3D anisotropic media with thin shapes are presented to verify the methodology developed above. The effect of the proposed method as well as the stability, and the convergence with respect to the number of boundary elements are carefully investigated. Unless otherwise specified, the thermal conductivity tensor kij are chosen as kij = 12 + 12 δij , where δij is the Kronecker delta. For the ease of comparison and validation of the numerical results, the analytical solution is taken to be 1
1
3
1
1
x2 x3 + 2. (31) 2 2 2 2 2 The numerical values obtained will be compared with analytical solutions in terms of the relative error defined by u (x) =
x21 +
x22 −
[ Relative Error =
x23 −
M 1 ∑
M
k=1
(
x1 x2 + x1 x3 +
k k Inumerical − Iexact k Iexact
)2 ]1/2
,
(32)
k k where Inumerical and Iexact are the numerical and analytical solutions at the kth calculation point, respectively.
4.1. Test problem 1: heat conduction in a thin hollow ring torus Firstly, we consider the steady state heat conduction in a thin hollow ring torus, as shown in Fig. 2. In mathematics, a torus can be defined parametrically by x (φ, θ) = (f (φ) cos θ, f (φ) sin θ, r sin φ) ,
(33)
where f (φ) = 4 + cos (φ), and φ and θ vary in the interval [0, 2π ). The inner radius of the ring torus is a = 1 which is constant in this study, while the outer radius b changes from 1.1 to 1.000000001. It can be seen that the thickness of the structure (δ ) changes in the ranges of 10−1 ≤ δ ≤ 10−9 . In this 3D model, Dirichlet boundary conditions are specified on the outer surface of the ring torus, while the Neumann boundary conditions are specified at the inner surface of the ring torus. The actual values of the Dirichlet and Neumann boundary conditions can be easily obtained using the corresponding analytical solutions given by Eq. (31). Here, a total number of N = 300 evenly distributed boundary elements, in terms of the angle parameterization, are chosen on the whole surfaces, regardless of the thickness of the structure. Figs. 3 and 4 illustrate the relative error curves of u (x) and ∂ u (x) /∂ x1 calculated at the interior point (4 + (a + b) /2, 0, 0), as functions of different thickness-to-length ratio. Both the untransformed BEM and the sinh-transformed BEM proposed in this paper are employed for the purpose of comparisons. As can be seen in Figs. 3 and 4, when the thickness of the structure is not very small, both the BEM schemes are effective and can yield accurate numerical results. However, the untransformed BEM performs less satisfactorily as the thickness decreases, i.e., after the thickness is less than 1E−1. In contrast, the results obtained using the proposed method are excellently consistent with the exact solutions, even when the thickness is in the order of 10−9 . In addition, the number of BEM elements does not change across the entire range of thickness, and therefore, the solution time and memory requirements are quite modest. Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
7
Fig. 3. Relative error curves of temperatures as functions of different thicknesses of the structure.
Fig. 4. Relative error curves of fluxes as functions of different thicknesses of the structure.
4.2. Test problem 2: heat conduction in a thin hollow peanut-shaped domain Next, we consider steady state heat conduction in a thin hollow peanut-shaped domain, as shown in Fig. 5. The geometry of the outer and inner surfaces of the computational domain is defined parametrically as x (φ, θ) = ((f (φ) + r ) sin φ cos θ, (f (φ) + r ) sin φ sin θ, (f (φ) + r ) cos φ) ,
(34)
x (φ, θ) = (f (φ) sin φ cos θ, f (φ) sin φ sin θ, f (φ) cos φ) ,
(35)
√
√
respectively, where φ ∈ [0, π) , θ ∈ [0, 2π), f (φ) = cos (2φ) + 1.1 − sin2 (2φ) and r = 10−6 denoting the thickness of the structure. Here, Neumann boundary conditions are specified on the outer surface of the structure and Dirichlet boundary conditions on the inner surface of the structure. To investigate the convergence of method, Figs. 6 and 7 display the relative error surfaces of u (x) and ∂ u (x) /∂ x1 using N = 400 and N = 800 evenly distributed boundary elements (in terms of the angle parameterization), respectively. In the case of N = 400 boundary elements, there are 25 elements discretized along the θ -direction (along the peanut-shaped Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
8
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 5. Steady state heat conduction in a peanut-shaped domain.
a
b
10-4
0.012
8
0.01 Relative errors
Relative errors
6 0.008 0.006 0.004 0.002 30
0 30
4
2
0 30
20 20
10
10 0
30
20 20 10
0
10 0
0
Fig. 6. Relative error surfaces of temperatures, where the number of BEM elements is 400 (a) and 800 (b), respectively.
domain) and 16 elements along the φ -direction (around the peanut-shaped domain). In the case of N = 800 boundary elements, there are 40 elements discretized along the θ -direction and 20 elements along the φ -direction. The error surfaces were yielded at M = 30 × 30 calculation points uniformly-spaced over the surface x (φ, θ) = ((f (φ) + r /2) sin φ cos θ, (f (φ) + r /2) sin φ sin θ, (f (φ) + r /2) cos φ) .
(36) −7
It is noted that the distance between the calculation points to the boundary is in the orders of 10 . As can be seen from Figs. 6 and 7, the present method is stable, accurate, and rapidly convergent as the number of boundary elements increases. It can be also observed that the BEM results with only 400 boundary elements are quite accurate for this 3D case, so that the size of resulting system of linear algebraic equations is quite small. 4.3. Test problem 3: heat conduction in a thin hollow hyperboloid domain Finally, we consider the steady state heat conduction in a thin hollow hyperboloid domain, as shown in Fig. 8. The geometry of the outer surface is defined as x (θ , t ) =
(√
1 + t 2 cos θ,
√
)
1 + t 2 sin θ, t ,
(37)
where the two parameters θ and t vary between 0 ∼ 2π and −2 ∼ 2, respectively. The inner boundary is taken to be a curve similar to the outer surface. Here the distance between the inner and outer surfaces changes from 10−1 to 10−8 . In this problem, Dirichlet boundary conditions are prescribed on the outer surface, while Neumann boundary conditions are specified at the remaining surface. To show the numerical accuracy, M = 100 interior points uniformly-spaced over the surface between the inner and outer boundaries are chosen. Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
a
)
–
9
b 10-3
0.015
1.5
0.01 Relative errors
Relative errors
2
0.005
0 30
1 0.5 0 30
30
20 20 10
10 0
0
30
20 20 10
10 0
0
Fig. 7. Relative error surfaces of fluxes, where the number of BEM elements is 400 (a) and 800 (b), respectively.
Fig. 8. Steady state heat conduction in a thin hollow hyperboloid domain.
Here, both the BEM and the FEM are employed for the purpose of comparison. For the BEM model, the total number of N = 600 boundary elements is chosen, regardless of the thickness of the structure. The FEM discretization varies according to the structure thickness. In general, more finite elements are required as the thickness of the structure decreases, in order to maintain reasonable element aspect ratio. Fig. 9 illustrates the relative error curves of the u (x) using both the BEM and the FEM, as functions of different values of thicknesses. The number of BEM and FEM elements employed in the calculation is provided in Fig. 10. As can be seen from Fig. 9, the BEM results are excellently consistent with the analytical solutions, even when the thickness of the structure is as small as 10−8 . And the number of BEM elements does not change across the entire range of thickness. The FEM solution, however, shows a very different behavior. When the thickness of the structure is less than 0.01, the problem is easily solved using the FEM. However, the FEM solution requires significantly more effort when the thickness of the structure further decreases due to aspect ratio limitations. Hence we can conclude that, in comparison with the FEM, the proposed BEM could be considered a competitive alternative for solving thin structural problems. 5. Conclusions This paper investigates the applicability of the BEM for three-dimensional anisotropic heat conduction problems on thin domains. It is shown that the BEM can deal with ultra-thin structures very efficiently and accurately, as long as the Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
10
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
Fig. 9. Temperature relative errors by using the BEM and FEM for different thicknesses of the structure.
Fig. 10. Numbers of elements of the BEM and FEM required for different thicknesses of the structure.
nearly singular integrals existing in the BEM formulations are handled correctly. It can be seen that the BEM does not suffer from thickness or aspect ratio issues associated with FEM. Therefore, the number of BEM elements can be held constant as described in examples. The method proposed in this paper can be applied to other 3D anisotropic problems on thin domains in the analogous way. It must be pointed out that the proposed BEM approach to thin-bodies should be considered as a complement to the FEM. For structures that can be identified clearly as shells or bulky solids, the FEM should be used because of the efficiency FEM delivers. For thin structures where the solid models are difficult to obtain, the BEM can be used as an alternative, especially when high accuracy is desired, as in benchmarks. Acknowledgments The work described in this paper was supported by the National Basic Research Program of China (973 Project No. 2010CB832702), the National Natural Science Foundation of China (Nos. 11402075, 71571108), Projects of International (Regional) Cooperation and Exchanges of NSFC (Nos. 71611530712, 61661136002), the China Postdoctoral Science Foundation (Nos. 2015M570572, 2016T90608), and the Qingdao Postdoctoral Application Research Project (No. 15-9-1-49-jc). Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
11
References [1] J.F. Luo, Y.J. Liu, E.J. Berger, Analysis of two-dimensional thin structures (from micro- to nano-scales) using the boundary element method, Comput. Mech. 22 (1998) 404–412. [2] B.P. Lamichhane, S.G. Roberts, L. Stals, A mixed finite element discretisation of thin plate splines based on biorthogonal systems, J. Sci. Comput. 67 (2016) 20–42. [3] L. Marin, A meshless method for solving the Cauchy problem in three-dimensional elastostatics, Comput. Math. Appl. 50 (2005) 73–92. [4] A.H.D. Cheng, D.T. Cheng, Heritage and early history of the boundary element method, Eng. Anal. Bound. Elem. 29 (2005) 268–302. [5] M.G. Knight, L.C. Wrobel, J.L. Henshall, Micromechanical response of fibre-reinforced materials using the boundary element technique, Compos. Struct. 62 (2003) 341–352. [6] Y. Gu, W. Chen, C. Zhang, Stress analysis for thin multilayered coating systems using a sinh transformed boundary element method, Internat. J. Solids Struct. 50 (2013) 3460–3471. [7] R.E. Caflisch, Growth, Structure and pattern formation for thin films, J. Sci. Comput. 37 (2008) 3–17. [8] B.A. Szabó, G.J. Sahrmann, Hierarchic plate and shell models based on p-extension, Int. J. Numer. Meth. Eng. 26 (1988) 1855–1881. [9] G. Zboinski, Adaptive hpq finite element methods for the analysis of 3D-based models of complex structures. Part 1. Hierarchical modeling and approximations, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2913–2940. [10] M. Guiggiani, Computing principal-value integrals in 3-D BEM for time-harmonic elastodynamics—A direct approach, Commun. Appl. Numer. Methods 8 (1992) 141–149. [11] Y. Gu, W. Chen, C.Z. Zhang, Singular boundary method for solving plane strain elastostatic problems, Internat. J. Solids Struct. 48 (2011) 2549–2556. [12] J.T. Chen, Y.L. Chang, S.K. Kao, J. Jian, Revisit of the indirect boundary element method: Necessary and sufficient formulation, J. Sci. Comput. 65 (2015) 467–485. [13] L. Demkowicz, Asymptotic convergence in finite and boundary element methods: part 1: theoretical results, Comput. Math. Appl. 27 (1994) 69–84. [14] L. Demkowicz, Asymptotic convergence in finite and boundary element methods: part 2: The LBB constant for rigid and elastic scattering problems, Comput. Math. Appl. 28 (1994) 93–109. [15] V. Sladek, J. Sladek, M. Tanaka, Nonsingular BEM formulations for thin-walled structures and elastostatic crack problems, Acta Mech. 99 (1993) 173–190. [16] Y.J. Liu, Analysis of shell-like structures by the Boundary Element Method based on 3-D elasticity: formulation and verification, Int. J. Numer. Meth. Eng. 41 (1998) 541–558. [17] A. Nagarajan, S. Mukherjee, A mapping method for numerical evaluation of two-dimensional integrals with 1/r singularity, Comput. Mech. 12 (1993) 19–26. [18] M. Guiggiani, The evaluation of Cauchy principal value integrals in the boundary element method—a review, Math. Comput. Modelling 15 (1991) 175–184. [19] J.C. Lachat, J.O. Watson, Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics, Int. J. Numer. Meth. Eng. 10 (1976) 991–1005. [20] X.W. Gao, K. Yang, J. Wang, An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals, Eng. Anal. Bound. Elem. 32 (2008) 692–696. [21] X.W. Gao, Numerical evaluation of two-dimensional singular boundary integrals—Theory and Fortran code, J. Comput. Appl. Math. 188 (2006) 44–64. [22] G.S. Padhi, R.A. Shenoi, S.S.J. Moy, M.A. McCarthy, Analytic integration of kernel shape function product integrals in the boundary element method, Comput. Struct. 79 (2001) 1325–1333. [23] H.L. Zhou, Z.R. Niu, C.Z. Cheng, Z.W. Guan, Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for anisotropic potential problems, Comput. Struct. 86 (2008) 1656–1671. [24] Z.R. Niu, C.Z. Cheng, H.L. Zhou, Z.J. Hu, Analytic formulations for calculating nearly singular integrals in two-dimensional BEM, Eng. Anal. Bound. Elem. 31 (2007) 949–964. [25] J.C.F. Telles, A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals, Int. J. Numer. Meth. Eng. 24 (1987) 959–973. [26] Q. Huang, T.A. Cruse, Some notes on singular integral techniques in boundary element analysis, Int. J. Numer. Meth. Eng. 36 (1993) 2643–2659. [27] P.R. Johnston, D. Elliott, A sinh transformation for evaluating nearly singular boundary element integrals, Int. J. Numer. Meth. Eng. 62 (2005) 564–578. [28] P.R. Johnston, B.M. Johnston, D. Elliott, Using the iterated sinh transformation to evaluate two dimensional nearly singular boundary element integrals, Eng. Anal. Bound. Elem. 37 (2013) 708–718. [29] V. Sladek, J. Sladek, M. Tanaka, Optimal transformations of the integration variables in computation of singular integrals in BEM, Int. J. Numer. Meth. Eng. 47 (2000) 1263–1283. [30] Y.M. Zhang, Y. Gu, J.T. Chen, Boundary element analysis of the thermal behaviour in thin-coated cutting tools, Eng. Anal. Bound. Elem. 34 (2010) 775–784. [31] G. Xie, J. Zhang, X. Qin, G. Li, New variable transformations for evaluating nearly singular integrals in 2D boundary element method, Eng. Anal. Bound. Elem. 35 (2011) 811–817. [32] Y. Gu, Q. Hua, W. Chen, C. Zhang, Numerical evaluation of nearly hyper-singular integrals in the boundary element analysis, Comput. Struct. 167 (2016) 15–23. [33] H. Ma, N. Kamiya, A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two- and threedimensional elasticity, Comput. Mech. 29 (2002) 277–288. [34] Y. Gu, W. Chen, B. Zhang, W. Qu, Two general algorithms for nearly singular integrals in two dimensional anisotropic boundary element method, Comput. Mech. 53 (2014) 1223–1234. [35] B.M. Johnston, P.R. Johnston, D. Elliott, A new method for the numerical evaluation of nearly singular integrals on triangular elements in the 3D boundary element method, J. Comput. Appl. Math. 245 (2013) 148–161. [36] P.R. Johnston, Application of sigmoidal transformations to weakly singular and near-singular boundary element integrals, Int. J. Numer. Meth. Eng. 45 (1999) 1333–1348. [37] Y. Gu, W. Chen, C. Zhang, The sinh transformation for evaluating nearly singular boundary element integrals over high-order geometry elements, Eng. Anal. Bound. Elem. 37 (2013) 301–308. [38] J. Lv, Y. Miao, H. Zhu, The distance sinh transformation for the numerical evaluation of nearly singular integrals over curved surface elements, Comput. Mech. 53 (2014) 359–367. [39] Y. Gu, H. Gao, W. Chen, C. Zhang, A general algorithm for evaluating nearly singular integrals in anisotropic three-dimensional boundary element analysis, Comput. Methods Appl. Mech. Engrg. 308 (2016) 483–498. [40] Y. Gu, W. Chen, X.Q. He, Singular boundary method for steady-state heat conduction in three dimensional general anisotropic media, Int. J. Heat. Mass. Transfer 55 (2012) 4837–4848.
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.
12
Y. Gu et al. / Computers and Mathematics with Applications (
)
–
[41] Y.M. Zhang, Z.Y. Liu, J.T. Chen, Y. Gu, A novel boundary element approach for solving the anisotropic potential problems, Eng. Anal. Bound. Elem. 35 (2011) 1181–1189. [42] P.K. Banerjee, The Boundary Element Methods in Engineering, McGRAW-HILL Book Company Europe, 1994. [43] M. Guiggiani, P. Casalini, Direct computation of Cauchy principal value integrals in advanced boundary elements, Int. J. Numer. Meth. Eng. 24 (1987) 1711–1720. [44] V. Sladek, J. Sladek, Singular integrals and boundary elements, Comput. Methods Appl. Mech. Engrg. 157 (1998) 251–266. [45] M. Guiggiani, A. Gigante, A general algorithm for multidimensional Cauchy principal value integrals in the boundary element method, J. Appl. Mech. 57 (1990) 906–915. [46] B.M. Johnston, P.R. Johnston, D. Elliott, A sinh transformation for evaluating two-dimensional nearly singular boundary element integrals, Int. J. Numer. Meth. Eng 69 (2007) 1460–1479.
Please cite this article in press as: Y. Gu, et al., Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method, Computers and Mathematics with Applications (2017), http://dx.doi.org/10.1016/j.camwa.2017.08.030.