- Email: [email protected]
The generalized finite difference method for an inverse boundary value problem in three-dimensional thermo-elasticity
Advances in Engineering Software 131 (2019) 1–11
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
The generalized finite difference method for an inverse boundary value problem in three-dimensional thermo-elasticity
T
Hu Wena, Gu Yana, , Zhang Chuanzengb, He Xiaoqiaoc ⁎
a
School of Mathematics and Statistics, Qingdao University, Qingdao 266071, PR China Department of Civil Engineering, University of Siegen, Paul-Bonatz-Str. 9-11, D-57076, Siegen, Germany c Department of Civil and Architectural Engineering, City University of Hong Kong, Hong Kong b
ARTICLE INFO
ABSTRACT
Keywords: Generalized finite difference method Meshless method Inverse problems Three-dimensional coupled thermo-elasticity
In this study, a new framework for the efficient and accurate solutions for an inverse problem associated with three-dimensional (3D) coupled thermo-elasticity equation is presented. The ill-conditioned problem is solved here with the generalized finite difference method (GFDM) together with the first-order Tikhonov regularization technique and the L-curve criterion. The GFDM uses the Taylor series expansions and the moving least squares approximation to derive explicit formulae for the required partial derivatives of unknown variables. The method is truly meshless that can be applied for solving problems merely defined over irregular clouds of points. In addition, for problems involving complex geometries, a new distance criterion for adaptive selection of nodes in the GFDM simulations is proposed. Preliminary numerical experiments show that the regularized GFDM proposed in this study are very promising for accurate and efficient inverse thermo-elasticity simulations, even with a comparatively large level of noise added into the input data.
1. Introduction During many machining processes, the thermo-elastic behavior of solids would strongly influence the tool performance, such as tool wear, tool life as well as process quality [1–6]. Thus, having a clear understanding about the thermal and elastic properties of solids is very useful and important, especially in the fields of automobile, aircraft as well as cutting-tool manufacturing technologies. The well established and widely applied finite element method (FEM) offers without doubt many advantages in solving thermo-elasticity problems for its flexibilities in dealing with the geometry and loading type. The FEM itself, however, has also many inherent shortcomings especially in the re-meshing processes and when the elements become highly distortable [7–12]. The boundary element method (BEM) is an important alternative to alleviate the difficulty associated with domain discretization, due to the boundary-only discretizations and its semi-analytical nature [13–15]. However, for problems with general body forces, an extra domain integral equation is required which have to be dealt with using special attention and may weaken the advantages of the BEM [13, 16–20]. In this study, the numerical technique we choose to solve the thermo-elasticity problem is the generalized finite difference method (GFDM) which offers several advantages over the traditional discretization methods. In particular, the GFDM is meshless in the sense ⁎
that only an irregular clouds of points is required for the discretization of the problem under investigation [4, 21–27]. Unlike the BEM, no potentially troublesome singular or nearly singular integration is required in the GFDM. Unlike the FEM, the irregular clouds of points used in the GFDM need not have any connectivity and the method may be applied in almost exactly the same way to any other arbitrary and irregular domains since the data preparation required and the computational details involved are similar. These features make the method very easy to implement, in particular for problems in complex geometries and high dimensions. The main idea of the GFDM is to combine the Taylor series expansions and the moving-least squares (MLS) approximation to derive explicit formulae for the required partial derivatives of unknown variables. The derivatives of unknown variables at a node, and then, can be approximated by a linear combination of function values with respect to its neighboring nodes. The main idea of the method was proposed in the early eighties by Lizska and Orkisz [28, 29] and later essentially extended and improved by many other authors [30-33]. Prior to this study, this method has been successfully tried for 2D and 3D parabolic and hyperbolic equations [32, 34], third- and fourth-order partial differential equations [35], dynamic analysis of beams and plates [36], non-linear elliptic partial differential equations [37], and applied inverse problems [33, 38].
Corresponding author. E-mail address: [email protected] (Y. Gu).
https://doi.org/10.1016/j.advengsoft.2019.02.006 Received 2 October 2018; Received in revised form 5 February 2019; Accepted 25 February 2019 0965-9978/ © 2019 Elsevier Ltd. All rights reserved.
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Motivated by the rapidly growing interest in the area, this paper documents the first attempt to apply the GFDM for the solution of inverse problems associated with the 3D thermo-elasticity equations. In inverse problems, one or more of the data describing the direct problem is missing. To fully determine the process, additional data must be supplied, either other boundary conditions on the same accessible part of boundary or measurements at some internal points in the domain [38–40]. The inverse problems are, generally, difficult to solve numerically due to the fact that they are ill-posed in the sense that small errors in measured data may lead arbitrarily large changes in the numerical solution [41]. The resulting ill-conditioned algebraic equations is regularized here by employing the first-order Tikhonov regularization technique [42], while the optimal regularization parameter is determined by the L-curve criterion. A brief outline of the rest of the paper is as follows. In Section 2, the mathematical formulation for 3D coupled thermo-elasticity problems is briefly introduced. The methodology of the GFDM and its numerical implementation for general 3D problems are reviewed in Section 3. A new distance criterion, named as twenty regions criterion, is also introduced for the adaptive selection of nodes in the GFDM. The Tikhonov regularization method with the choice of the regularization parameter given by the L-curve criterion is presented in Section 4. In Section 5, three benchmark numerical examples involving complicated geometries are presented. Finally, some conclusions and remarks are provided in Section 6.
ij
ti (x ) =
c2
c2
2u 1 x 22
+
2u 1
x1 x2
2u 2 x12
+
2u 1
+ c2
x1 x3
+
2u 1 x 32
+ c1 2u 2
2u 2 x 22
x1 x2
+ c2
2u 2 x 32
+
2u 3 x12
+
x2 x3
2u 2
+ c2
+
+ c2
2u 3 x 22
2u 3
k
x12
2T
+
x 22
+
2T
x32
2 1
2 , 2
c2 =
2
,
=
1 2
uj (x ) ui (x ) + , xj xi
x
1,
(9)
x
1,
(10)
where the barred quantities T˜ , u˜ i , q˜ , and t˜i indicate the given values specified on the boundary, q(x) = −(k∇T(x)) • n(x) denotes the normal heat flux at the point x. In the above formulation, it can be seen that the boundary Γ1 is over-specified, whilst the remaining boundary Γ2 is under-specified since both the temperatures, displacements and tractions are unknown and have to be determined. In the case of the inverse BVP given by Eqs. (1)–(4), (9) and (10), the uniqueness of the solution was proved by Kozlov et al. [45], but the problem is still ill-posed. In order to solve the resulting ill-conditioned algebraic equations, suitable regularization procedures should be employed. 3. The generalized finite difference method (GFDM)
(1)
x2 x3
T = f2 (x ), x2
(2)
2u 3 x 32
T = f3 (x ), x3
Without loss of generality, let's consider the following general differential equation in the 3D domain [46]
(3)
a1
2u 3
+ c1
3.1. Explicit formulae in GFDM
2G (1 + ) , 1 2
U x1
+ a2 + a8
U x2
+ a3
2U x 1 x3
U x3
+ a9
+ a4
2U x2 x3
2U
x12
+ a5
= f (x ),
2U
x22
+ a6
2U
x32
+ a7
2U x1 x2
(11)
or for brevity
L2 [U ] = f (x ),
(12)
where L2[U] is a second-order partial differential operator, f(x) is known function, ai, i = 1, 2, ..., 9, are constants. In order to obtain the explicit GFDM formulae for partial differential equations, an irregular cloud of points is scattered in the computational domain, as shown in Fig. 1. For each given node x0, named as the central node, the m nearest nodes xi (i = 1, 2, ..., m), called the neighbors or supporting nodes, will be found within a prescribed distance dm from the central node x0, i.e., |xi − x0| ≤ dm. According to Gavete et al. [47], the concept of the ‘star’ then refers to the area of all supporting nodes in relation to the central node. Note that each node scattered inside the computational domain has an associated star assigned. Suppose U0 is the value of the function at the central node x0 and Ui, i = 1, 2, ..., m, are function values at the rest of the nodes inside the star. Then, expanding the values of Ui around the central point x0 using the Taylor series expansion, one has [34, 46]
(5)
are thermo-elastic constants in which ν, G, and α being Poisson's ratio, shear modulus and coefficient of linear thermal expansion, respectively. The kinematics of deformation is described by the linear strain tensor ij (x )
(8)
ij (x ) nj (x ),
ui (x ) = u˜ i (x ) ti (x ) = t˜i (x )
(4)
=
(7)
T = f1 (x ), x1
x1 x3
= g (x ),
1 1
T ij,
kk ij
2
and
where k > 0 is the thermal conductivity, fi(x), i = 1, 2, 3, and g(x) denote the known body forces and internal heat sources, respectively, and
c1 =
1
T (x ) = T˜ (x ) q (x ) = (k T (x ))·n (x )
and the steady-state heat conduction equation 2T
+
where nj(x) are the outward unit normal vector. In many practical situations, only a part of the boundary, say Γ1, is accessible for measurements, while the remaining boundary Γ2 is under-specified since no boundary data is available on it. In this case, additional measurements should be given on Γ1, and this corresponds to an inverse boundary value problem (BVP). In this study, we consider the inverse BVPs consisting of Eqs. (1)–(4) and the following thermal and mechanical boundary conditions
Consider a 3D domain Ω bounded by a surface Γ = ∂Ω which may consist of several surfaces, each being sufficiently smooth in the sense of Liapunov. We here assume that the surfaces consist of two parts, Γ = Γ1∪Γ2, where Γ1,Γ2 ≠ ∅ and Γ1∩Γ2 = ∅. The coupled equilibrium equations for 3D static thermo-elasticity problems, in terms of the temperature T(x) and displacements ui(x), i = 1, 2, 3, satisfying the following Navier–Lamé equations [4, 43, 44] 2u 1 x12
ij
where δij is the Kronecker delta. The boundary tractions ti are then defined in terms of the stresses as follows
2. Statement of the 3D thermo-elasticity problems
c1
= 2G
(6)
where sufficiently small displacements and displacement gradients are assumed. Here and in the following, the customary Einstein's notation for summation over repeated subscripts is applied. The stresses σij(x) are related to the strains εij(x) through Hooke's law by 2
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 1. An irregular cloud of points and the selection of stars for (a) 2D case; and (b) 3D case.
Ui = U0 + hi + hi k i
U0 x1
+ ki
2U 0 x 1 x2
U0 x2
+ li
+ h i li
U0 x3
2U 0 x 1 x3
+
1 2
+ ki l i
hi2
2U 0 x12
2U 0 x2 x3
2U 0 x22
+ ki2
+ ......,
2U 0 x32
+ li2
P=
i = 1, 2, ...,m .
ki = x i2
x 02 ,
li = x i3
m i=1
B (U ) = +
U0
Ui + hi
U0 x1
+ ki
U0 x2
U0 x3
+ li
+
+
6
x1 x2
+ hi l i
2U 0
x1 x3
+ k i li
2U 0
2
+8
di dm
U0 U0 U0 , , , x1 x2 x3
x2 x3
, A=
(21)
(22)
3
3
di dm
=
di
2U 0 , x12
2U 0 , x 22
2U 0 , x32
dm .
2U 0
x1 x2
,
hi2 ki2 li2 , , , h i ki, h i li , k i li , 2 2 2
k1 k2
k1 l1 k2 l2
T
2 1
h1 h2
2 2
km lm
2 m
k1 k2
hm k m
k1 l1 k2 l2 km lm (23)
PT WP ,
and
4
,
h1 h2
hm k m
(16)
b = P T W (U
2U 0
x1 x3
,
2U 0
x2 x3
i = 1, 2, ...m ,
(24)
U0 ).
According to Eqs. (22)–(24), the partial derivative vector DU as defined in Eq. (17) can now be expressed as
DU = A 1 b=A 1 PT W (U
U0)
m
=A
1 i=1
2 T i pi ,
2 T 1 p1 ,
,
2 T m pm
U0 U1 U2 . Um
(25)
According to the above analysis, for each node x0 inside the computational domain, the derivatives of the unknown function at x0 can now be approximated by a linear combination of function values at its neighboring nodes (see Eq. (25)). Substituting Eq. (25) into the partial differential equation L2[U] leads to
T
, (17)
pi = hi , ki, li,
T
where
2 i
where di = ‖xi − x0‖ is the distance between points x0 and xi, dm stands for the distance between nodes x0 and the farthest node inside the star. It is interesting to note that the weighting coefficients are inversely related to the distance from the corresponding node to the central node of the star, indicating that the Taylor series approximation is more important if the point xi is closer to the central node x0. Of course other choices of weighting function are possible, see Refs. [37, 46]. Let, for brevity,
DU =
(20)
U ),
ADU = b,
2U 0
di dm
U )T W (PDU + U0
then the following linear equation system can be obtained (for each node of the discretization)
ki2 2U0 2 x22
where ωi = ω(hi,ki,li) denotes the weighting coefficients at point xi. In our computations, the weighting coefficients are taken to be
=1
(19)
km lm
B (U ) = 0, {DU }
(15)
i
hm k m
.
where U ={U1, U2, ⋅⋅⋅, Um} , U0={U0, U0, ⋅⋅⋅, U0} and W = diag ( 12 , 22, , m2 ) . By minimizing the residual function B(U) with respect to the partial derivatives DU, i.e.,
li2 2U0 2 x32
+ h i ki
k1 l1 k2 l2
T
When the Taylor series (13) is truncated after the second-order derivatives, it is then possible to define the following residual function B(U) as hi2 2U0 2 x12
k1 k2
B (U ) = (PDU + U0
(14)
x 03.
h1 h2
Then the residual function B(U) defined in Eq. (15) can be expressed in the following matrix notation
Employing indicial notation for the coordinates of points x0 and xi, i.e., x 0 = (x 01, x 02 , x 03) and xi = (x i1, x i2 , x i3) , respectively, the coefficients in the Taylor series expansion (13) can be written as
x 01,
=
pm
(13)
hi = xi1
p1 p2
m
L2 [U0] = m 0 U0 +
(18)
mi Ui, i=1
and
(26)
where the coefficients m0 and mi (i = 1, 2, ..., m) can be immediately 3
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
obtained from Eq. (25). Repeating the same procedure to each one of the nodes inside the computational domain, a system of linear algebraic equations can then be established. On solving this equation system, the numerical results of the function U at each node inside the domain can be obtained. It is interesting to note that the concept of the star utilized in the GFDM yields a sparse matrix system, and therefore, the final linear algebraic equations can be solved quickly by using various sparse matrix solvers. 3.2. Criteria for adaptive selection of nodes inside each of the stars The location of nodes selected inside the star has a great influence on the GFDM results. In the early applications of the GFDM, researchers usually only considered the distance criterion in which m nearest nodes around the central node are selected. This criterion, however, may produce distorted (ill-conditioned) stars with unevenly distributed nodes which may lead to degenerated solutions, especially for problems with complicated geometries or problems involving very irregular clouds of nodes. To overcome the drawbacks associated with the distance criterion, Benito et al. [48] proposed four quadrants criterion and eight octants criterion for 2D and 3D cases, respectively. In 2D case, the area around the central node is divided into four sectors corresponding to quadrants of the Cartesian coordinates system originating at the central node. In each sector two or more nodes are selected according to the distance criterion. Similarly, in 3D case the volume around the central node is divided in eight sectors corresponding to octants of the Cartesian coordinates system. This approach is accurate, stable and corrects the problems associated with the traditional distance criterion. With this idea in mind, for 3D cases, we here introduce a similar distance criterion, named as twenty regions criterion, in order to further correct the effect of node density irregularity. It is interesting to note that the x-, y- and z-axes of a 3D Cartesian system divide the space into eight regions (or octants) and twelve semi-infinite planes each of which bounded by two half-axes (see Fig. 2). Therefore, in this twenty regions criterion, the space around the central node is split into twenty subregions, according to the aforementioned eight regions (octants) and twelve semi-infinite planes. In each sub-region two or more nodes are selected, the closest to the central node, as shown in Fig. 3. In order to avoid ill-conditioning clouds of nodes, a maximum distance is used as an indicator. If the distance between the new node and the central node is larger than this maximum distance, then the new node should not be
Fig. 3. The selection of two nodes in one of octants and semi-infinite planes, respectively.
added. This maximum distance is an important parameter involved in the control of the GFDM results and it is given as the scale of the domain, multiplied by a positive parameter α < 1. The scale of the domain can be indicated by the maximum distance between all the nodes of the mesh. 4. Regularization techniques for inverse problems As stated in Section 2, the inverse problem is ill-posed and we cannot use a direct approach, e.g. Gaussian elimination method, to solve the system of linear equations in a straightforward manner. Hence regularization techniques are necessary to filter out the influence of the noise. One widely used regularization technique for solving ill-conditioning problem is the Tikhonov regularization method, which considers a minimum of the following function
f (x ) = Ax
b
2 2
+
2
2
L(i ) x 2 , L(i)
R (N
i) × N ,
i = 0, 1, 2, ..., (27)
where λ is the regularization parameter, ‖ • ‖2denotes the Euclidean norm, L(i) is a matrix that defines a (semi) norm of solution vector in which the superscript (i) represents the ith derivative operator on L, for example, in the case of the first-order Tikhonov regularization method the matrix L(i) is given by
L(1) =
0 0
1
1 0
0 1 1
0 0 1 1
. (N 1) × N
(28)
N
Solving ∇fλ(x) = 0 for x ∈ R , we can obtain the Tikhonov regularized solution xλ of the Eq. (27) which is given as the solution of the regularized equation
(AT A +
2L (i)T L(i) ) x =AT b .
(29)
The performance of the Tikhonov method depends crucially on a suitable choice of the regularization parameter λ. The L-curve criterion is frequently used. The L-curve is a log-log plot of the norm of regularized solution (‖L(i)x‖2) versus the norm of the corresponding residual norm (‖Ax − b‖2), i.e., as a curve
(log Ax Fig. 2. The twenty regions criterion.
b 2 , log L(i) x 2 ),
(30)
parameterized by the regularization parameter. The horizontal part of 4
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 4. Geometry of the problem (a) and the configuration of the GFDM nodes distribution (b).
Fig. 5. Numerical results for displacements u1 on the under-specified boundary obtained using BL = 0.5 and various levels of data noise.
the L-curve corresponds to the index of how smooth the solution is treated, and the vertical part of the L-curve corresponds to the distance index between the predicted output and real output. The corner point of L-curve is a compromise between the regularization errors due to data smoothing and perturbation errors in measurements or other noise. Hence, the L-curve is really a trade-off curve between two quantities that both should be controlled and, according to the L-curve criterion, the optimal value of the regularization parameter is chosen at the “corner” of the L-curve, see Ref. [39].
(31)
T (x ) = cos(x1) + cos (x2) + cos (x3), for temperatures,
u1 (x ) = x1,
u2 (x ) = x2 ,
(32)
u3 (x ) = x3,
for displacements, and 2G (1 + ) 1 2
11 (x )
=
22 (x )
=
33 (x )
=
12 (x )
=
13 (x )
=
23 (x )
= 0,
1
1
2
T (x ) , (33)
for stresses. Unless otherwise specified, the thermo-elastic constants as introduced in Section 2 are taken to be k = 1 (thermal conductivity), ν = 0.2 (Poisson's ratio), G = 13889N/m2 (shear modulus), and α = 10−5 (coefficient of linear thermal expansion). In order to evaluate the performance of the numerical method, an L2 error norm is defined as follows:
5. Numerical results and discussions In this section, three benchmark numerical examples associated with 3D thermo-elasticity equations are presented to verify the methodologies developed above. The effect of regularization as well as the stability of the scheme with respect to the noise added into the data are carefully investigated. All the computations were implemented on a Core(TM) II laptop PC with a 2.90 GHz CPU, 16 G RAM and 500 G hard drive. For the ease of comparison and validation of the numerical results, the analytical solutions are given by
1/2
M
Global
Error =
[Inumer (k ) k=1
Iexact (k )]2
1/2
M
[Iexact (k )]2
/
,
k=1
(34) 5
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 6. Numerical results for temperatures T on the under-specified boundary obtained using BL = 0.5 and various levels of data noise.
where Inumerical (k ) and Iexact (k ) denote the numerical and analytical solutions at the kth calculated point, respectively, M is the number of calculation points tested. In our test cases, the boundary temperatures and displacements, prescribed on the over-specified boundary Γ1, have been perturbed as
b˜ = b 1 + rand ×
100
,
noise (3%) added into the input data. It is also shown that, as expected, the GFDM solutions converge to their corresponding exact solutions as the amount of noise decreases. A similar conclusion can be drawn from Fig. 7(a)−(d) which present the numerical results retrieved for tractions t1 on the under-specified boundary Γ2, with BL = 0.5 and various levels of noise added into the input data. Although not presented here for the sake of brevity, it should be noted that similar results have been obtained for other components of the displacement vector and stress tensor, as well as for the displacements and stresses in the cases with other combinations of the boundary conditions. It can be conclude that the GFDM, in conjunction with the Tikhonov regularization technique, can provide accurate and stable numerical solutions for this 3D case.
(35)
where b is the exact boundary data, rand ∈ [ − 1, 1] is a pseudorandom number which is generated using the MATLAB command ‘2*rand(•)−1′, and δ represents the level of noise added into the input data. Unless otherwise specified, the measure of the over-specified boundary is defined as Γ1 = BL × Γ in which
BL =
measure ( 1) , measure ( )
BL
(0, 1],
5.2. Test problem 2: Thermo-elastic analysis in a nuclear submarine
(36)
stands for the ratio parameter of the accessible boundary.
Next, we perform the 3D thermo-elastic analysis in a nuclear submarine. The geometry of the problem is shown in Fig. 8, with the global dimension 60 m in length, 8 m in width, and 9 m in height. To solve the problem numerically, a total number of 6259 irregularly distributed nodes, as shown in Fig. 9, are discretized inside the whole computational domain. Here, we investigate the influence of the length of the over-specified boundary segment Γ1 on the accuracy of the numerical solutions. To do so, Fig. 10 shows the relative error curves of the calculated displacements u1, temperatures T and tractions t1 retrieved on the whole surface of the computational domain, with 1% noise and various values of parameter BL. It can be seen from this figure that the numerical results retrieved on the under-specified boundary agree pretty well with their corresponding analytical solutions, and a small accessible part of the boundary, for example BL = 0.55, is sufficient for obtaining accurate numerical results. Moreover, it can be observed that, as expected, the accuracy of the numerical results improves significantly as the length of the accessible boundary increases, i.e. BL increases. Similar results have been obtained also for the other components of the displacements and stresses analyzed, but they are not presented here for the sake of brevity. Fig. 11 shows the L-curves for Tikhonov regularization techniques, with BL = 0.8 and 1% noise added into the data. The L-curve is a loglog plot of the norm of regularized solution versus the norm of the corresponding residual norm. The corner point of L-curve is a compromise between the regularization errors due to data smoothing and perturbation errors in measurements or other noise. According to the L-
5.1. Test problem 1: a five pieces fan blade for turbofan engines We first consider the thermo-elastic analysis in a five pieces fan blade as shown in Fig. 4. The principal dimension of the fan blade is 5 m in length, 1.5 m in width, and 5 m in height. For the numerical calculation, a total number of 4650 irregularly distributed nodes are discretized inside the whole computational domain. The geometry model and the automatic generation of nodes are realized here using the CAE software Hypermesh. The automatic generation of a mesh using the Hypermesh packages is often done in two stages: the generation of nodes and then the generation of elements. In the GFDM the second stage, usually more troublesome, may be omitted. This feature makes the method very easy to implement, in particular for problems in complex geometries and three dimensions. Furthermore, the GFDM may be applied in almost exactly the same way to any other arbitrary and irregular domains since the data preparation required and the computational details involved are similar. The over-specified boundary Γ1 is taken to be the left-half surface of the boundary ∂Ω with BL = 0.5, i.e., Γ1 = { − 2.5 ≤ y ≤ 0}. In Figs. 5 and 6 we illustrate the numerical results for displacements u1(x) and temperatures T(x) retrieved on the under-specified boundary Γ2, with various levels of noise added into the input data. The results presented in Figs. 5 and 6 illustrate that the GFDM solutions retrieved for both displacements and temperatures are in good agreement with their corresponding exact solutions, even with a relatively large amount of 6
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 7. Relative errors of the computed tractions t1 retrieved on the whole surface of the computational domain with BL = 0.5 and 0% noise (a), 1% noise (b), 2% noise (c), and 3% noise (d).
Fig. 8. The sketch of a nuclear submarine.
curve criterion, the optimal value of the regularization parameter is chosen at the “corner” of the L-curve.
implementation, a total number of 10,568 irregularly distributed nodes, as shown in Fig. 12(b), are discretized inside the whole computational domain. Fig. 13 illustrates the numerical results for temperatures T(x) retrieved on the under-specified boundary Γ2, with BL = 0.8 and various levels of noise added into the input data. The results presented in Fig. 13 illustrate that the GFDM solutions retrieved for temperatures are in good agreement with their corresponding exact solutions, even with a relatively large amount of noise (3%) added into the input data.
5.3. Test problem 3: thermo-elastic analysis in an airplane Finally, the inverse problem in an airplane is considered. The geometry of the problem is shown in Fig. 12(a) with dimension 33.2 m in length, 22.5 m in width, and 7.6 m in height. For the numerical 7
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 9. Configuration of the GFDM nodes distribution for the nuclear submarine.
Fig. 10. Relative error curves of the computed displacements, temperatures, and tractions with respect to the length of the accessible boundary BL.
Fig. 11. L-curves for Tikhonov regularization technique, with BL=0.8 and 1% noise added into the data.
8
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 12. Geometry of the problem (a) and the configuration of the GFDM nodes distribution (b).
Fig. 13. Numerical results for temperatures T on the under-specified boundary obtained using BL = 0.8 and various levels of data noise.
Fig. 14. The analytical solution profile of the displacement amplitudes |u(x)| on the surface of the airplane.
Figs. 14 and 15 respectively plot the exact and numerical solutions for ‘displacement amplitudes |u(x)|’ retrieved on the whole surface of the airplane. The ‘displacement amplitudes |u(x)|’ studied here is defined as
u (x ) =
u12 (x ) + u22 (x ) + u32 (x ) .
The numerical displacement amplitudes here are calculated with 2% noise added into the input data and BL = 0.8. As can be seen from these figures, the displacement amplitudes predicted by using the proposed scheme are in quite good agreement with the analytical solution. Although not presented, it is reported that numerous other numerical experiments have been performed and analogous conclusions
(37) 9
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 15. Numerical results for displacement amplitudes |u(x)| retrieved on the whole surface of the airplane, with BL = 0.8 and 2% noise added into the input data.
have been drawn. Overall, it can be concluded that the GFDM, in conjunction with the first-order Tikhonov regularization technique and the L-curve criterion, is accurate and stable with respect to decreasing the amount of noise added into the input data and convergent with respect to increasing the number of boundary nodes. In comparison with existing methods for solving numerically inverse problems in 3D thermo-elasticity, the proposed scheme could be considered as a competitive alternative.
2004;28:323–32. [2] Liu YJ, Shen L. A dual BIE approach for large-scale modelling of 3-D electrostatic problems with the fast multipole boundary element method. Int J Numer Methods Eng 2007;71:837–55. [3] Marin L, Karageorghis A, Lesnic D. Regularized MFS solution of inverse boundary value problems in three-dimensional steady-state linear thermoelasticity. Int J Solids Struct 2016;91:127–42. [4] Karageorghis A, Lesnic D, Marin L. The method of fundamental solutions for an inverse boundary value problem in static thermo-elasticity. Comput Struct 2014;135:32–9. [5] Al-shyyab AS, Darwish FH, Al-Nimr MA. Thermoelastic response of the solid boundaries of a semi-infinite rectangular micro-channel flow. Appl Therm Eng 2016;93:580–9. [6] Wredenberg F, Larsson P-L. On the effect of substrate deformation at scratching of soft thin film composites. Int J Mech Sci 2010;52:1008–14. [7] Rabizadeh E, Bagherzadeh AS, Rabczuk T. Goal-oriented error estimation and adaptive mesh refinement in dynamic coupled thermoelasticity. Comput Struct 2016;173:187–211. [8] Gu Y, He X, Chen W, Zhang C. Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method. Comput Math Appl 2018;75:33–44. [9] Boštjan M, Božidar Š. Application of the RBF collocation method to transient coupled thermoelasticity. Int J Numer Methods Heat Fluid Flow 2017;27:1064–77. [10] Fu Z-J, Xi Q, Chen W, Cheng AHD. A boundary-type meshless solver for transient heat conduction analysis of slender functionally graded materials with exponential variations. Comput Math Appl 2018;76:760–73. [11] Chen CS, Karageorghis A, Li Y. On choosing the location of the sources in the MFS. Numer Algorithms 2016;72:107–30. [12] Fu Z, Chen W, Wen P, Zhang C. Singular boundary method for wave propagation analysis in periodic structures. J Sound Vib 2018;425:170–88. [13] Cheng AHD, Cheng DT. Heritage and early history of the boundary element method. Eng Anal Bound Elem 2005;29:268–302. [14] Karageorghis A, Lesnic D, Marin L. A moving pseudo-boundary MFS for void detection in two-dimensional thermoelasticity. Int J Mech Sci 2014;88:276–88. [15] Šarler B, Kuhn G. Dual reciprocity boundary element method for convective-diffusive solid-liquid phase change problems, Part 1. Formulation. Eng Anal Bound Elem 1998;21:53–63. [16] Gao XW. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals. Comput Methods Appl Mech Eng 2010;199:2856–64. [17] Gu Y, Gao H, Chen W, Zhang C. A general algorithm for evaluating nearly singular integrals in anisotropic three-dimensional boundary element analysis. Comput Methods Appl Mech Eng 2016;308:483–98. [18] Zhang A, Gu Y, Hua Q, Chen W, Zhang C. A regularized singular boundary method for inverse Cauchy problem in three-dimensional elastostatics. Adv Appl Mathd Mech 2018;10:1459–77. [19] Alarcon E, Brebbia C, Dominguez J. The boundary element method in elasticity. Int J Mech Sci 1978;20:625–39. [20] Gu Y, Fan C-M, Xu R-P. Localized method of fundamental solutions for large-scale modelling of two-dimensional elasticity problems. Appl. Math. Lett. 2019:8–14. https://doi.org/10.1016/j.aml.2019.01.035. [21] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [22] Liu GR, Zhang GY, Wang YY, Zhong ZH, Li GY, Han X. A nodal integration technique for meshfree radial point interpolation method (NI-RPIM). Int J Solids Struct 2007;44:3840–60. [23] Gu Y, Chen W, Gao H, Zhang C. A meshless singular boundary method for threedimensional elasticity problems. Int J Numer Methods Eng 2016;107:109–26. [24] Sadat H, Dubus N, Gbahoué L, Sophy T. On the solution of heterogeneous heat conduction problems by a diffuse approximation meshless method. Numer Heat Tranf B-Fundam 2006;50(6):491–8. [25] Qu W, Chen W. Solution of two-dimensional stokes flow problems using improved
6. Concluding remarks This study documents the first attempt to apply the GFDM for the efficient and accurate solutions of inverse problems in 3D linear thermo-elasticity. The resulting ill-conditioned system of linear algebraic equations has been regularized by using the Tiknonov regularization technique, while the optimal regularization parameter was determined by the L-curve criterion. Our preliminary analysis and numerical experiments show that the GFDM methods are very promising for inverse problem simulations. However, further analyses and optimized implementation are required in order to fully explore the efficiency and accuracy of the new method. These include the detailed convergence order analyses of the method, optimized strategies for adaptive mesh refinements in both time and space, as well as the optimal choice of many different parameters. Results along these lines will be presented in the future. Finally, it must be pointed out that the proposed GFDM to 3D inverse thermo-elasticity problems should be considered as a complement to the popular FEM in the structural analysis. For structures that can be identified clearly as shells or bulky solids, the FEM should be used because of the widely available automatic meshing capabilities in the various FEM packages. For 3D engineering structures where the solid models are difficult to obtain, the proposed GFDM can be used as an alternative, especially when high accuracy is desired, as in many benchmark solutions. Acknowledgments The work described in this paper was supported by the National Natural Science Foundation of China (Nos. 11872220, 71571108, 11572112), Projects of International (Regional) Cooperation and Exchanges of NSFC (No. 71611530712), and the Natural Science Foundation of Shandong Province of China (Nos. ZR2017JL004, ZR2017BA003, ZR2017MF055). References [1] Copetti MIM, French DA. Numerical studies of the stability of steady-state solutions to a contact problem in coupled thermoelasticity. Appl Math Model
10
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al. singular boundary method. Adv Appl Math Mech 2015;7:13–30. [26] Qu W, Chen W, Gu Y. Fast multipole accelerated singular boundary method for the 3D Helmholtz equation in low frequency regime. Comput Math Appl 2015;70:679–90. [27] Sarler B, Vertnik R. Meshfree explicit local radial basis function collocation method for diffusion problems. Comput Math Appl 2006;51:1269–82. [28] Liszka T. An interpolation method for an irregular net of nodes. Int J Numer Methods Eng 1984;20:1599–612. [29] Liszka T, Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980;11:83–95. [30] Payre GMJ. Influence graphs and the generalized finite difference method. Comput Meth Appl Mech Eng 2007;196:1933–45. [31] Ureña F, Benito JJ, Salete E, Gavete L. A note on the application of the generalized finite difference method to seismic wave propagation in 2D. J Comput Appl Math 2012;236:3016–25. [32] Benito JJ, Urena F, Gavete L. Influence of several factors in the generalized finite difference method. Appl Math Model 2001;25:1039–53. [33] Gu Y, Wang L, Chen W, Zhang C, He X. Application of the meshless generalized finite difference method to inverse heat source problems. Int J Heat Mass Transf 2017;108(Part A):721–9. [34] Hua Q, Gu Y, Qu W, Chen W, Zhang C. A meshless generalized finite difference method for inverse Cauchy problems associated with three-dimensional inhomogeneous Helmholtz-type equations. Eng Anal Bound Elem 2017;82:162–71. [35] Ureña F, Salete E, Benito JJ, Gavete L. Solving third- and fourth-order partial differential equations using GFDM: application to solve problems of plates. Int J Comput Math 2012;89:366–76. [36] Gavete L, Urena F, Benito JJ, Salete E. A note on the dynamic analysis using the generalized finite difference method. J Comput Appl Math 2013;252:132–47. [37] Gavete L, Ureña F, Benito JJ, García A, Ureña M, Salete E. Solving second order
[38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]
11
non-linear elliptic partial differential equations using generalized finite difference method. J Comput Appl Math 2017;318:378–87. Fan CM, Huang YK, Li PW, Chiu CL. Application of the generlized finite-difference method to inverse biharmonic boundary-value problems. Numer Heat Tranf BFundam 2014;65:129–54. Marin L, Karageorghis A. The MFS for the Cauchy problem in two-dimensional steady-state linear thermoelasticity. Int J Solids Struct 2013;50:3387–98. Gu Y, Chen W, Zhang C, He X. A meshless singular boundary method for threedimensional inverse heat conduction problems in general anisotropic media. Int J Heat Mass Transf 2015;84:91–102. Karageorghis A, Lesnic D, Marin L. A survey of applications of the MFS to inverse problems. Inverse Probl Sci Eng 2011;19:309–36. Marin L, Elliott L, Ingham DB, Lesnic D. An iterative boundary element algorithm for a singular Cauchy problem in linear elasticity. Comput Mech 2002;28:479–88. Tanaka M, Matsumoto T, Moradi M. Application of boundary element method to 3D problems of coupled thermoelasticity. Eng Anal Bound Elem 1995;16:297–303. Zheng BJ, Gao XW, Yang K, Zhang CZ. A novel meshless local Petrov–Galerkin method for dynamic coupled thermoelasticity analysis under thermal and mechanical shock loading. Eng Anal Bound Elem 2015;60:154–61. Kozlov VA, Maz'ya VG, Fomin AV. Uniqueness of the solution to an inverse thermoelasticity problem. Comput Math Math Phys 2009;49:525–31. Gavete L, Benito JJ, Ureña F. Generalized finite differences for solving 3D elliptic and parabolic equations. Appl Math Model 2016;40:955–65. Gavete L, Gavete ML, Benito JJ. Improvements of generalized finite difference method and comparison with other meshless method. Appl Math Model 2003;27:831–47. Benito JJ, Urena F, Gavete L. Solving parabolic and hyperbolic equations by the generalized finite difference method. J Comput Appl Math 2007;209:208–33.
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
The generalized finite difference method for an inverse boundary value problem in three-dimensional thermo-elasticity
T
Hu Wena, Gu Yana, , Zhang Chuanzengb, He Xiaoqiaoc ⁎
a
School of Mathematics and Statistics, Qingdao University, Qingdao 266071, PR China Department of Civil Engineering, University of Siegen, Paul-Bonatz-Str. 9-11, D-57076, Siegen, Germany c Department of Civil and Architectural Engineering, City University of Hong Kong, Hong Kong b
ARTICLE INFO
ABSTRACT
Keywords: Generalized finite difference method Meshless method Inverse problems Three-dimensional coupled thermo-elasticity
In this study, a new framework for the efficient and accurate solutions for an inverse problem associated with three-dimensional (3D) coupled thermo-elasticity equation is presented. The ill-conditioned problem is solved here with the generalized finite difference method (GFDM) together with the first-order Tikhonov regularization technique and the L-curve criterion. The GFDM uses the Taylor series expansions and the moving least squares approximation to derive explicit formulae for the required partial derivatives of unknown variables. The method is truly meshless that can be applied for solving problems merely defined over irregular clouds of points. In addition, for problems involving complex geometries, a new distance criterion for adaptive selection of nodes in the GFDM simulations is proposed. Preliminary numerical experiments show that the regularized GFDM proposed in this study are very promising for accurate and efficient inverse thermo-elasticity simulations, even with a comparatively large level of noise added into the input data.
1. Introduction During many machining processes, the thermo-elastic behavior of solids would strongly influence the tool performance, such as tool wear, tool life as well as process quality [1–6]. Thus, having a clear understanding about the thermal and elastic properties of solids is very useful and important, especially in the fields of automobile, aircraft as well as cutting-tool manufacturing technologies. The well established and widely applied finite element method (FEM) offers without doubt many advantages in solving thermo-elasticity problems for its flexibilities in dealing with the geometry and loading type. The FEM itself, however, has also many inherent shortcomings especially in the re-meshing processes and when the elements become highly distortable [7–12]. The boundary element method (BEM) is an important alternative to alleviate the difficulty associated with domain discretization, due to the boundary-only discretizations and its semi-analytical nature [13–15]. However, for problems with general body forces, an extra domain integral equation is required which have to be dealt with using special attention and may weaken the advantages of the BEM [13, 16–20]. In this study, the numerical technique we choose to solve the thermo-elasticity problem is the generalized finite difference method (GFDM) which offers several advantages over the traditional discretization methods. In particular, the GFDM is meshless in the sense ⁎
that only an irregular clouds of points is required for the discretization of the problem under investigation [4, 21–27]. Unlike the BEM, no potentially troublesome singular or nearly singular integration is required in the GFDM. Unlike the FEM, the irregular clouds of points used in the GFDM need not have any connectivity and the method may be applied in almost exactly the same way to any other arbitrary and irregular domains since the data preparation required and the computational details involved are similar. These features make the method very easy to implement, in particular for problems in complex geometries and high dimensions. The main idea of the GFDM is to combine the Taylor series expansions and the moving-least squares (MLS) approximation to derive explicit formulae for the required partial derivatives of unknown variables. The derivatives of unknown variables at a node, and then, can be approximated by a linear combination of function values with respect to its neighboring nodes. The main idea of the method was proposed in the early eighties by Lizska and Orkisz [28, 29] and later essentially extended and improved by many other authors [30-33]. Prior to this study, this method has been successfully tried for 2D and 3D parabolic and hyperbolic equations [32, 34], third- and fourth-order partial differential equations [35], dynamic analysis of beams and plates [36], non-linear elliptic partial differential equations [37], and applied inverse problems [33, 38].
Corresponding author. E-mail address: [email protected] (Y. Gu).
https://doi.org/10.1016/j.advengsoft.2019.02.006 Received 2 October 2018; Received in revised form 5 February 2019; Accepted 25 February 2019 0965-9978/ © 2019 Elsevier Ltd. All rights reserved.
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Motivated by the rapidly growing interest in the area, this paper documents the first attempt to apply the GFDM for the solution of inverse problems associated with the 3D thermo-elasticity equations. In inverse problems, one or more of the data describing the direct problem is missing. To fully determine the process, additional data must be supplied, either other boundary conditions on the same accessible part of boundary or measurements at some internal points in the domain [38–40]. The inverse problems are, generally, difficult to solve numerically due to the fact that they are ill-posed in the sense that small errors in measured data may lead arbitrarily large changes in the numerical solution [41]. The resulting ill-conditioned algebraic equations is regularized here by employing the first-order Tikhonov regularization technique [42], while the optimal regularization parameter is determined by the L-curve criterion. A brief outline of the rest of the paper is as follows. In Section 2, the mathematical formulation for 3D coupled thermo-elasticity problems is briefly introduced. The methodology of the GFDM and its numerical implementation for general 3D problems are reviewed in Section 3. A new distance criterion, named as twenty regions criterion, is also introduced for the adaptive selection of nodes in the GFDM. The Tikhonov regularization method with the choice of the regularization parameter given by the L-curve criterion is presented in Section 4. In Section 5, three benchmark numerical examples involving complicated geometries are presented. Finally, some conclusions and remarks are provided in Section 6.
ij
ti (x ) =
c2
c2
2u 1 x 22
+
2u 1
x1 x2
2u 2 x12
+
2u 1
+ c2
x1 x3
+
2u 1 x 32
+ c1 2u 2
2u 2 x 22
x1 x2
+ c2
2u 2 x 32
+
2u 3 x12
+
x2 x3
2u 2
+ c2
+
+ c2
2u 3 x 22
2u 3
k
x12
2T
+
x 22
+
2T
x32
2 1
2 , 2
c2 =
2
,
=
1 2
uj (x ) ui (x ) + , xj xi
x
1,
(9)
x
1,
(10)
where the barred quantities T˜ , u˜ i , q˜ , and t˜i indicate the given values specified on the boundary, q(x) = −(k∇T(x)) • n(x) denotes the normal heat flux at the point x. In the above formulation, it can be seen that the boundary Γ1 is over-specified, whilst the remaining boundary Γ2 is under-specified since both the temperatures, displacements and tractions are unknown and have to be determined. In the case of the inverse BVP given by Eqs. (1)–(4), (9) and (10), the uniqueness of the solution was proved by Kozlov et al. [45], but the problem is still ill-posed. In order to solve the resulting ill-conditioned algebraic equations, suitable regularization procedures should be employed. 3. The generalized finite difference method (GFDM)
(1)
x2 x3
T = f2 (x ), x2
(2)
2u 3 x 32
T = f3 (x ), x3
Without loss of generality, let's consider the following general differential equation in the 3D domain [46]
(3)
a1
2u 3
+ c1
3.1. Explicit formulae in GFDM
2G (1 + ) , 1 2
U x1
+ a2 + a8
U x2
+ a3
2U x 1 x3
U x3
+ a9
+ a4
2U x2 x3
2U
x12
+ a5
= f (x ),
2U
x22
+ a6
2U
x32
+ a7
2U x1 x2
(11)
or for brevity
L2 [U ] = f (x ),
(12)
where L2[U] is a second-order partial differential operator, f(x) is known function, ai, i = 1, 2, ..., 9, are constants. In order to obtain the explicit GFDM formulae for partial differential equations, an irregular cloud of points is scattered in the computational domain, as shown in Fig. 1. For each given node x0, named as the central node, the m nearest nodes xi (i = 1, 2, ..., m), called the neighbors or supporting nodes, will be found within a prescribed distance dm from the central node x0, i.e., |xi − x0| ≤ dm. According to Gavete et al. [47], the concept of the ‘star’ then refers to the area of all supporting nodes in relation to the central node. Note that each node scattered inside the computational domain has an associated star assigned. Suppose U0 is the value of the function at the central node x0 and Ui, i = 1, 2, ..., m, are function values at the rest of the nodes inside the star. Then, expanding the values of Ui around the central point x0 using the Taylor series expansion, one has [34, 46]
(5)
are thermo-elastic constants in which ν, G, and α being Poisson's ratio, shear modulus and coefficient of linear thermal expansion, respectively. The kinematics of deformation is described by the linear strain tensor ij (x )
(8)
ij (x ) nj (x ),
ui (x ) = u˜ i (x ) ti (x ) = t˜i (x )
(4)
=
(7)
T = f1 (x ), x1
x1 x3
= g (x ),
1 1
T ij,
kk ij
2
and
where k > 0 is the thermal conductivity, fi(x), i = 1, 2, 3, and g(x) denote the known body forces and internal heat sources, respectively, and
c1 =
1
T (x ) = T˜ (x ) q (x ) = (k T (x ))·n (x )
and the steady-state heat conduction equation 2T
+
where nj(x) are the outward unit normal vector. In many practical situations, only a part of the boundary, say Γ1, is accessible for measurements, while the remaining boundary Γ2 is under-specified since no boundary data is available on it. In this case, additional measurements should be given on Γ1, and this corresponds to an inverse boundary value problem (BVP). In this study, we consider the inverse BVPs consisting of Eqs. (1)–(4) and the following thermal and mechanical boundary conditions
Consider a 3D domain Ω bounded by a surface Γ = ∂Ω which may consist of several surfaces, each being sufficiently smooth in the sense of Liapunov. We here assume that the surfaces consist of two parts, Γ = Γ1∪Γ2, where Γ1,Γ2 ≠ ∅ and Γ1∩Γ2 = ∅. The coupled equilibrium equations for 3D static thermo-elasticity problems, in terms of the temperature T(x) and displacements ui(x), i = 1, 2, 3, satisfying the following Navier–Lamé equations [4, 43, 44] 2u 1 x12
ij
where δij is the Kronecker delta. The boundary tractions ti are then defined in terms of the stresses as follows
2. Statement of the 3D thermo-elasticity problems
c1
= 2G
(6)
where sufficiently small displacements and displacement gradients are assumed. Here and in the following, the customary Einstein's notation for summation over repeated subscripts is applied. The stresses σij(x) are related to the strains εij(x) through Hooke's law by 2
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 1. An irregular cloud of points and the selection of stars for (a) 2D case; and (b) 3D case.
Ui = U0 + hi + hi k i
U0 x1
+ ki
2U 0 x 1 x2
U0 x2
+ li
+ h i li
U0 x3
2U 0 x 1 x3
+
1 2
+ ki l i
hi2
2U 0 x12
2U 0 x2 x3
2U 0 x22
+ ki2
+ ......,
2U 0 x32
+ li2
P=
i = 1, 2, ...,m .
ki = x i2
x 02 ,
li = x i3
m i=1
B (U ) = +
U0
Ui + hi
U0 x1
+ ki
U0 x2
U0 x3
+ li
+
+
6
x1 x2
+ hi l i
2U 0
x1 x3
+ k i li
2U 0
2
+8
di dm
U0 U0 U0 , , , x1 x2 x3
x2 x3
, A=
(21)
(22)
3
3
di dm
=
di
2U 0 , x12
2U 0 , x 22
2U 0 , x32
dm .
2U 0
x1 x2
,
hi2 ki2 li2 , , , h i ki, h i li , k i li , 2 2 2
k1 k2
k1 l1 k2 l2
T
2 1
h1 h2
2 2
km lm
2 m
k1 k2
hm k m
k1 l1 k2 l2 km lm (23)
PT WP ,
and
4
,
h1 h2
hm k m
(16)
b = P T W (U
2U 0
x1 x3
,
2U 0
x2 x3
i = 1, 2, ...m ,
(24)
U0 ).
According to Eqs. (22)–(24), the partial derivative vector DU as defined in Eq. (17) can now be expressed as
DU = A 1 b=A 1 PT W (U
U0)
m
=A
1 i=1
2 T i pi ,
2 T 1 p1 ,
,
2 T m pm
U0 U1 U2 . Um
(25)
According to the above analysis, for each node x0 inside the computational domain, the derivatives of the unknown function at x0 can now be approximated by a linear combination of function values at its neighboring nodes (see Eq. (25)). Substituting Eq. (25) into the partial differential equation L2[U] leads to
T
, (17)
pi = hi , ki, li,
T
where
2 i
where di = ‖xi − x0‖ is the distance between points x0 and xi, dm stands for the distance between nodes x0 and the farthest node inside the star. It is interesting to note that the weighting coefficients are inversely related to the distance from the corresponding node to the central node of the star, indicating that the Taylor series approximation is more important if the point xi is closer to the central node x0. Of course other choices of weighting function are possible, see Refs. [37, 46]. Let, for brevity,
DU =
(20)
U ),
ADU = b,
2U 0
di dm
U )T W (PDU + U0
then the following linear equation system can be obtained (for each node of the discretization)
ki2 2U0 2 x22
where ωi = ω(hi,ki,li) denotes the weighting coefficients at point xi. In our computations, the weighting coefficients are taken to be
=1
(19)
km lm
B (U ) = 0, {DU }
(15)
i
hm k m
.
where U ={U1, U2, ⋅⋅⋅, Um} , U0={U0, U0, ⋅⋅⋅, U0} and W = diag ( 12 , 22, , m2 ) . By minimizing the residual function B(U) with respect to the partial derivatives DU, i.e.,
li2 2U0 2 x32
+ h i ki
k1 l1 k2 l2
T
When the Taylor series (13) is truncated after the second-order derivatives, it is then possible to define the following residual function B(U) as hi2 2U0 2 x12
k1 k2
B (U ) = (PDU + U0
(14)
x 03.
h1 h2
Then the residual function B(U) defined in Eq. (15) can be expressed in the following matrix notation
Employing indicial notation for the coordinates of points x0 and xi, i.e., x 0 = (x 01, x 02 , x 03) and xi = (x i1, x i2 , x i3) , respectively, the coefficients in the Taylor series expansion (13) can be written as
x 01,
=
pm
(13)
hi = xi1
p1 p2
m
L2 [U0] = m 0 U0 +
(18)
mi Ui, i=1
and
(26)
where the coefficients m0 and mi (i = 1, 2, ..., m) can be immediately 3
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
obtained from Eq. (25). Repeating the same procedure to each one of the nodes inside the computational domain, a system of linear algebraic equations can then be established. On solving this equation system, the numerical results of the function U at each node inside the domain can be obtained. It is interesting to note that the concept of the star utilized in the GFDM yields a sparse matrix system, and therefore, the final linear algebraic equations can be solved quickly by using various sparse matrix solvers. 3.2. Criteria for adaptive selection of nodes inside each of the stars The location of nodes selected inside the star has a great influence on the GFDM results. In the early applications of the GFDM, researchers usually only considered the distance criterion in which m nearest nodes around the central node are selected. This criterion, however, may produce distorted (ill-conditioned) stars with unevenly distributed nodes which may lead to degenerated solutions, especially for problems with complicated geometries or problems involving very irregular clouds of nodes. To overcome the drawbacks associated with the distance criterion, Benito et al. [48] proposed four quadrants criterion and eight octants criterion for 2D and 3D cases, respectively. In 2D case, the area around the central node is divided into four sectors corresponding to quadrants of the Cartesian coordinates system originating at the central node. In each sector two or more nodes are selected according to the distance criterion. Similarly, in 3D case the volume around the central node is divided in eight sectors corresponding to octants of the Cartesian coordinates system. This approach is accurate, stable and corrects the problems associated with the traditional distance criterion. With this idea in mind, for 3D cases, we here introduce a similar distance criterion, named as twenty regions criterion, in order to further correct the effect of node density irregularity. It is interesting to note that the x-, y- and z-axes of a 3D Cartesian system divide the space into eight regions (or octants) and twelve semi-infinite planes each of which bounded by two half-axes (see Fig. 2). Therefore, in this twenty regions criterion, the space around the central node is split into twenty subregions, according to the aforementioned eight regions (octants) and twelve semi-infinite planes. In each sub-region two or more nodes are selected, the closest to the central node, as shown in Fig. 3. In order to avoid ill-conditioning clouds of nodes, a maximum distance is used as an indicator. If the distance between the new node and the central node is larger than this maximum distance, then the new node should not be
Fig. 3. The selection of two nodes in one of octants and semi-infinite planes, respectively.
added. This maximum distance is an important parameter involved in the control of the GFDM results and it is given as the scale of the domain, multiplied by a positive parameter α < 1. The scale of the domain can be indicated by the maximum distance between all the nodes of the mesh. 4. Regularization techniques for inverse problems As stated in Section 2, the inverse problem is ill-posed and we cannot use a direct approach, e.g. Gaussian elimination method, to solve the system of linear equations in a straightforward manner. Hence regularization techniques are necessary to filter out the influence of the noise. One widely used regularization technique for solving ill-conditioning problem is the Tikhonov regularization method, which considers a minimum of the following function
f (x ) = Ax
b
2 2
+
2
2
L(i ) x 2 , L(i)
R (N
i) × N ,
i = 0, 1, 2, ..., (27)
where λ is the regularization parameter, ‖ • ‖2denotes the Euclidean norm, L(i) is a matrix that defines a (semi) norm of solution vector in which the superscript (i) represents the ith derivative operator on L, for example, in the case of the first-order Tikhonov regularization method the matrix L(i) is given by
L(1) =
0 0
1
1 0
0 1 1
0 0 1 1
. (N 1) × N
(28)
N
Solving ∇fλ(x) = 0 for x ∈ R , we can obtain the Tikhonov regularized solution xλ of the Eq. (27) which is given as the solution of the regularized equation
(AT A +
2L (i)T L(i) ) x =AT b .
(29)
The performance of the Tikhonov method depends crucially on a suitable choice of the regularization parameter λ. The L-curve criterion is frequently used. The L-curve is a log-log plot of the norm of regularized solution (‖L(i)x‖2) versus the norm of the corresponding residual norm (‖Ax − b‖2), i.e., as a curve
(log Ax Fig. 2. The twenty regions criterion.
b 2 , log L(i) x 2 ),
(30)
parameterized by the regularization parameter. The horizontal part of 4
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 4. Geometry of the problem (a) and the configuration of the GFDM nodes distribution (b).
Fig. 5. Numerical results for displacements u1 on the under-specified boundary obtained using BL = 0.5 and various levels of data noise.
the L-curve corresponds to the index of how smooth the solution is treated, and the vertical part of the L-curve corresponds to the distance index between the predicted output and real output. The corner point of L-curve is a compromise between the regularization errors due to data smoothing and perturbation errors in measurements or other noise. Hence, the L-curve is really a trade-off curve between two quantities that both should be controlled and, according to the L-curve criterion, the optimal value of the regularization parameter is chosen at the “corner” of the L-curve, see Ref. [39].
(31)
T (x ) = cos(x1) + cos (x2) + cos (x3), for temperatures,
u1 (x ) = x1,
u2 (x ) = x2 ,
(32)
u3 (x ) = x3,
for displacements, and 2G (1 + ) 1 2
11 (x )
=
22 (x )
=
33 (x )
=
12 (x )
=
13 (x )
=
23 (x )
= 0,
1
1
2
T (x ) , (33)
for stresses. Unless otherwise specified, the thermo-elastic constants as introduced in Section 2 are taken to be k = 1 (thermal conductivity), ν = 0.2 (Poisson's ratio), G = 13889N/m2 (shear modulus), and α = 10−5 (coefficient of linear thermal expansion). In order to evaluate the performance of the numerical method, an L2 error norm is defined as follows:
5. Numerical results and discussions In this section, three benchmark numerical examples associated with 3D thermo-elasticity equations are presented to verify the methodologies developed above. The effect of regularization as well as the stability of the scheme with respect to the noise added into the data are carefully investigated. All the computations were implemented on a Core(TM) II laptop PC with a 2.90 GHz CPU, 16 G RAM and 500 G hard drive. For the ease of comparison and validation of the numerical results, the analytical solutions are given by
1/2
M
Global
Error =
[Inumer (k ) k=1
Iexact (k )]2
1/2
M
[Iexact (k )]2
/
,
k=1
(34) 5
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 6. Numerical results for temperatures T on the under-specified boundary obtained using BL = 0.5 and various levels of data noise.
where Inumerical (k ) and Iexact (k ) denote the numerical and analytical solutions at the kth calculated point, respectively, M is the number of calculation points tested. In our test cases, the boundary temperatures and displacements, prescribed on the over-specified boundary Γ1, have been perturbed as
b˜ = b 1 + rand ×
100
,
noise (3%) added into the input data. It is also shown that, as expected, the GFDM solutions converge to their corresponding exact solutions as the amount of noise decreases. A similar conclusion can be drawn from Fig. 7(a)−(d) which present the numerical results retrieved for tractions t1 on the under-specified boundary Γ2, with BL = 0.5 and various levels of noise added into the input data. Although not presented here for the sake of brevity, it should be noted that similar results have been obtained for other components of the displacement vector and stress tensor, as well as for the displacements and stresses in the cases with other combinations of the boundary conditions. It can be conclude that the GFDM, in conjunction with the Tikhonov regularization technique, can provide accurate and stable numerical solutions for this 3D case.
(35)
where b is the exact boundary data, rand ∈ [ − 1, 1] is a pseudorandom number which is generated using the MATLAB command ‘2*rand(•)−1′, and δ represents the level of noise added into the input data. Unless otherwise specified, the measure of the over-specified boundary is defined as Γ1 = BL × Γ in which
BL =
measure ( 1) , measure ( )
BL
(0, 1],
5.2. Test problem 2: Thermo-elastic analysis in a nuclear submarine
(36)
stands for the ratio parameter of the accessible boundary.
Next, we perform the 3D thermo-elastic analysis in a nuclear submarine. The geometry of the problem is shown in Fig. 8, with the global dimension 60 m in length, 8 m in width, and 9 m in height. To solve the problem numerically, a total number of 6259 irregularly distributed nodes, as shown in Fig. 9, are discretized inside the whole computational domain. Here, we investigate the influence of the length of the over-specified boundary segment Γ1 on the accuracy of the numerical solutions. To do so, Fig. 10 shows the relative error curves of the calculated displacements u1, temperatures T and tractions t1 retrieved on the whole surface of the computational domain, with 1% noise and various values of parameter BL. It can be seen from this figure that the numerical results retrieved on the under-specified boundary agree pretty well with their corresponding analytical solutions, and a small accessible part of the boundary, for example BL = 0.55, is sufficient for obtaining accurate numerical results. Moreover, it can be observed that, as expected, the accuracy of the numerical results improves significantly as the length of the accessible boundary increases, i.e. BL increases. Similar results have been obtained also for the other components of the displacements and stresses analyzed, but they are not presented here for the sake of brevity. Fig. 11 shows the L-curves for Tikhonov regularization techniques, with BL = 0.8 and 1% noise added into the data. The L-curve is a loglog plot of the norm of regularized solution versus the norm of the corresponding residual norm. The corner point of L-curve is a compromise between the regularization errors due to data smoothing and perturbation errors in measurements or other noise. According to the L-
5.1. Test problem 1: a five pieces fan blade for turbofan engines We first consider the thermo-elastic analysis in a five pieces fan blade as shown in Fig. 4. The principal dimension of the fan blade is 5 m in length, 1.5 m in width, and 5 m in height. For the numerical calculation, a total number of 4650 irregularly distributed nodes are discretized inside the whole computational domain. The geometry model and the automatic generation of nodes are realized here using the CAE software Hypermesh. The automatic generation of a mesh using the Hypermesh packages is often done in two stages: the generation of nodes and then the generation of elements. In the GFDM the second stage, usually more troublesome, may be omitted. This feature makes the method very easy to implement, in particular for problems in complex geometries and three dimensions. Furthermore, the GFDM may be applied in almost exactly the same way to any other arbitrary and irregular domains since the data preparation required and the computational details involved are similar. The over-specified boundary Γ1 is taken to be the left-half surface of the boundary ∂Ω with BL = 0.5, i.e., Γ1 = { − 2.5 ≤ y ≤ 0}. In Figs. 5 and 6 we illustrate the numerical results for displacements u1(x) and temperatures T(x) retrieved on the under-specified boundary Γ2, with various levels of noise added into the input data. The results presented in Figs. 5 and 6 illustrate that the GFDM solutions retrieved for both displacements and temperatures are in good agreement with their corresponding exact solutions, even with a relatively large amount of 6
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 7. Relative errors of the computed tractions t1 retrieved on the whole surface of the computational domain with BL = 0.5 and 0% noise (a), 1% noise (b), 2% noise (c), and 3% noise (d).
Fig. 8. The sketch of a nuclear submarine.
curve criterion, the optimal value of the regularization parameter is chosen at the “corner” of the L-curve.
implementation, a total number of 10,568 irregularly distributed nodes, as shown in Fig. 12(b), are discretized inside the whole computational domain. Fig. 13 illustrates the numerical results for temperatures T(x) retrieved on the under-specified boundary Γ2, with BL = 0.8 and various levels of noise added into the input data. The results presented in Fig. 13 illustrate that the GFDM solutions retrieved for temperatures are in good agreement with their corresponding exact solutions, even with a relatively large amount of noise (3%) added into the input data.
5.3. Test problem 3: thermo-elastic analysis in an airplane Finally, the inverse problem in an airplane is considered. The geometry of the problem is shown in Fig. 12(a) with dimension 33.2 m in length, 22.5 m in width, and 7.6 m in height. For the numerical 7
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 9. Configuration of the GFDM nodes distribution for the nuclear submarine.
Fig. 10. Relative error curves of the computed displacements, temperatures, and tractions with respect to the length of the accessible boundary BL.
Fig. 11. L-curves for Tikhonov regularization technique, with BL=0.8 and 1% noise added into the data.
8
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 12. Geometry of the problem (a) and the configuration of the GFDM nodes distribution (b).
Fig. 13. Numerical results for temperatures T on the under-specified boundary obtained using BL = 0.8 and various levels of data noise.
Fig. 14. The analytical solution profile of the displacement amplitudes |u(x)| on the surface of the airplane.
Figs. 14 and 15 respectively plot the exact and numerical solutions for ‘displacement amplitudes |u(x)|’ retrieved on the whole surface of the airplane. The ‘displacement amplitudes |u(x)|’ studied here is defined as
u (x ) =
u12 (x ) + u22 (x ) + u32 (x ) .
The numerical displacement amplitudes here are calculated with 2% noise added into the input data and BL = 0.8. As can be seen from these figures, the displacement amplitudes predicted by using the proposed scheme are in quite good agreement with the analytical solution. Although not presented, it is reported that numerous other numerical experiments have been performed and analogous conclusions
(37) 9
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al.
Fig. 15. Numerical results for displacement amplitudes |u(x)| retrieved on the whole surface of the airplane, with BL = 0.8 and 2% noise added into the input data.
have been drawn. Overall, it can be concluded that the GFDM, in conjunction with the first-order Tikhonov regularization technique and the L-curve criterion, is accurate and stable with respect to decreasing the amount of noise added into the input data and convergent with respect to increasing the number of boundary nodes. In comparison with existing methods for solving numerically inverse problems in 3D thermo-elasticity, the proposed scheme could be considered as a competitive alternative.
2004;28:323–32. [2] Liu YJ, Shen L. A dual BIE approach for large-scale modelling of 3-D electrostatic problems with the fast multipole boundary element method. Int J Numer Methods Eng 2007;71:837–55. [3] Marin L, Karageorghis A, Lesnic D. Regularized MFS solution of inverse boundary value problems in three-dimensional steady-state linear thermoelasticity. Int J Solids Struct 2016;91:127–42. [4] Karageorghis A, Lesnic D, Marin L. The method of fundamental solutions for an inverse boundary value problem in static thermo-elasticity. Comput Struct 2014;135:32–9. [5] Al-shyyab AS, Darwish FH, Al-Nimr MA. Thermoelastic response of the solid boundaries of a semi-infinite rectangular micro-channel flow. Appl Therm Eng 2016;93:580–9. [6] Wredenberg F, Larsson P-L. On the effect of substrate deformation at scratching of soft thin film composites. Int J Mech Sci 2010;52:1008–14. [7] Rabizadeh E, Bagherzadeh AS, Rabczuk T. Goal-oriented error estimation and adaptive mesh refinement in dynamic coupled thermoelasticity. Comput Struct 2016;173:187–211. [8] Gu Y, He X, Chen W, Zhang C. Analysis of three-dimensional anisotropic heat conduction problems on thin domains using an advanced boundary element method. Comput Math Appl 2018;75:33–44. [9] Boštjan M, Božidar Š. Application of the RBF collocation method to transient coupled thermoelasticity. Int J Numer Methods Heat Fluid Flow 2017;27:1064–77. [10] Fu Z-J, Xi Q, Chen W, Cheng AHD. A boundary-type meshless solver for transient heat conduction analysis of slender functionally graded materials with exponential variations. Comput Math Appl 2018;76:760–73. [11] Chen CS, Karageorghis A, Li Y. On choosing the location of the sources in the MFS. Numer Algorithms 2016;72:107–30. [12] Fu Z, Chen W, Wen P, Zhang C. Singular boundary method for wave propagation analysis in periodic structures. J Sound Vib 2018;425:170–88. [13] Cheng AHD, Cheng DT. Heritage and early history of the boundary element method. Eng Anal Bound Elem 2005;29:268–302. [14] Karageorghis A, Lesnic D, Marin L. A moving pseudo-boundary MFS for void detection in two-dimensional thermoelasticity. Int J Mech Sci 2014;88:276–88. [15] Šarler B, Kuhn G. Dual reciprocity boundary element method for convective-diffusive solid-liquid phase change problems, Part 1. Formulation. Eng Anal Bound Elem 1998;21:53–63. [16] Gao XW. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals. Comput Methods Appl Mech Eng 2010;199:2856–64. [17] Gu Y, Gao H, Chen W, Zhang C. A general algorithm for evaluating nearly singular integrals in anisotropic three-dimensional boundary element analysis. Comput Methods Appl Mech Eng 2016;308:483–98. [18] Zhang A, Gu Y, Hua Q, Chen W, Zhang C. A regularized singular boundary method for inverse Cauchy problem in three-dimensional elastostatics. Adv Appl Mathd Mech 2018;10:1459–77. [19] Alarcon E, Brebbia C, Dominguez J. The boundary element method in elasticity. Int J Mech Sci 1978;20:625–39. [20] Gu Y, Fan C-M, Xu R-P. Localized method of fundamental solutions for large-scale modelling of two-dimensional elasticity problems. Appl. Math. Lett. 2019:8–14. https://doi.org/10.1016/j.aml.2019.01.035. [21] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [22] Liu GR, Zhang GY, Wang YY, Zhong ZH, Li GY, Han X. A nodal integration technique for meshfree radial point interpolation method (NI-RPIM). Int J Solids Struct 2007;44:3840–60. [23] Gu Y, Chen W, Gao H, Zhang C. A meshless singular boundary method for threedimensional elasticity problems. Int J Numer Methods Eng 2016;107:109–26. [24] Sadat H, Dubus N, Gbahoué L, Sophy T. On the solution of heterogeneous heat conduction problems by a diffuse approximation meshless method. Numer Heat Tranf B-Fundam 2006;50(6):491–8. [25] Qu W, Chen W. Solution of two-dimensional stokes flow problems using improved
6. Concluding remarks This study documents the first attempt to apply the GFDM for the efficient and accurate solutions of inverse problems in 3D linear thermo-elasticity. The resulting ill-conditioned system of linear algebraic equations has been regularized by using the Tiknonov regularization technique, while the optimal regularization parameter was determined by the L-curve criterion. Our preliminary analysis and numerical experiments show that the GFDM methods are very promising for inverse problem simulations. However, further analyses and optimized implementation are required in order to fully explore the efficiency and accuracy of the new method. These include the detailed convergence order analyses of the method, optimized strategies for adaptive mesh refinements in both time and space, as well as the optimal choice of many different parameters. Results along these lines will be presented in the future. Finally, it must be pointed out that the proposed GFDM to 3D inverse thermo-elasticity problems should be considered as a complement to the popular FEM in the structural analysis. For structures that can be identified clearly as shells or bulky solids, the FEM should be used because of the widely available automatic meshing capabilities in the various FEM packages. For 3D engineering structures where the solid models are difficult to obtain, the proposed GFDM can be used as an alternative, especially when high accuracy is desired, as in many benchmark solutions. Acknowledgments The work described in this paper was supported by the National Natural Science Foundation of China (Nos. 11872220, 71571108, 11572112), Projects of International (Regional) Cooperation and Exchanges of NSFC (No. 71611530712), and the Natural Science Foundation of Shandong Province of China (Nos. ZR2017JL004, ZR2017BA003, ZR2017MF055). References [1] Copetti MIM, French DA. Numerical studies of the stability of steady-state solutions to a contact problem in coupled thermoelasticity. Appl Math Model
10
Advances in Engineering Software 131 (2019) 1–11
W. Hu, et al. singular boundary method. Adv Appl Math Mech 2015;7:13–30. [26] Qu W, Chen W, Gu Y. Fast multipole accelerated singular boundary method for the 3D Helmholtz equation in low frequency regime. Comput Math Appl 2015;70:679–90. [27] Sarler B, Vertnik R. Meshfree explicit local radial basis function collocation method for diffusion problems. Comput Math Appl 2006;51:1269–82. [28] Liszka T. An interpolation method for an irregular net of nodes. Int J Numer Methods Eng 1984;20:1599–612. [29] Liszka T, Orkisz J. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980;11:83–95. [30] Payre GMJ. Influence graphs and the generalized finite difference method. Comput Meth Appl Mech Eng 2007;196:1933–45. [31] Ureña F, Benito JJ, Salete E, Gavete L. A note on the application of the generalized finite difference method to seismic wave propagation in 2D. J Comput Appl Math 2012;236:3016–25. [32] Benito JJ, Urena F, Gavete L. Influence of several factors in the generalized finite difference method. Appl Math Model 2001;25:1039–53. [33] Gu Y, Wang L, Chen W, Zhang C, He X. Application of the meshless generalized finite difference method to inverse heat source problems. Int J Heat Mass Transf 2017;108(Part A):721–9. [34] Hua Q, Gu Y, Qu W, Chen W, Zhang C. A meshless generalized finite difference method for inverse Cauchy problems associated with three-dimensional inhomogeneous Helmholtz-type equations. Eng Anal Bound Elem 2017;82:162–71. [35] Ureña F, Salete E, Benito JJ, Gavete L. Solving third- and fourth-order partial differential equations using GFDM: application to solve problems of plates. Int J Comput Math 2012;89:366–76. [36] Gavete L, Urena F, Benito JJ, Salete E. A note on the dynamic analysis using the generalized finite difference method. J Comput Appl Math 2013;252:132–47. [37] Gavete L, Ureña F, Benito JJ, García A, Ureña M, Salete E. Solving second order
[38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48]
11
non-linear elliptic partial differential equations using generalized finite difference method. J Comput Appl Math 2017;318:378–87. Fan CM, Huang YK, Li PW, Chiu CL. Application of the generlized finite-difference method to inverse biharmonic boundary-value problems. Numer Heat Tranf BFundam 2014;65:129–54. Marin L, Karageorghis A. The MFS for the Cauchy problem in two-dimensional steady-state linear thermoelasticity. Int J Solids Struct 2013;50:3387–98. Gu Y, Chen W, Zhang C, He X. A meshless singular boundary method for threedimensional inverse heat conduction problems in general anisotropic media. Int J Heat Mass Transf 2015;84:91–102. Karageorghis A, Lesnic D, Marin L. A survey of applications of the MFS to inverse problems. Inverse Probl Sci Eng 2011;19:309–36. Marin L, Elliott L, Ingham DB, Lesnic D. An iterative boundary element algorithm for a singular Cauchy problem in linear elasticity. Comput Mech 2002;28:479–88. Tanaka M, Matsumoto T, Moradi M. Application of boundary element method to 3D problems of coupled thermoelasticity. Eng Anal Bound Elem 1995;16:297–303. Zheng BJ, Gao XW, Yang K, Zhang CZ. A novel meshless local Petrov–Galerkin method for dynamic coupled thermoelasticity analysis under thermal and mechanical shock loading. Eng Anal Bound Elem 2015;60:154–61. Kozlov VA, Maz'ya VG, Fomin AV. Uniqueness of the solution to an inverse thermoelasticity problem. Comput Math Math Phys 2009;49:525–31. Gavete L, Benito JJ, Ureña F. Generalized finite differences for solving 3D elliptic and parabolic equations. Appl Math Model 2016;40:955–65. Gavete L, Gavete ML, Benito JJ. Improvements of generalized finite difference method and comparison with other meshless method. Appl Math Model 2003;27:831–47. Benito JJ, Urena F, Gavete L. Solving parabolic and hyperbolic equations by the generalized finite difference method. J Comput Appl Math 2007;209:208–33.