On applying polyharmonic radial basis functions to solve 3D uncoupled anisotropic thermoelasticity problems using Boundary Element and Dual Reciprocity Method

Abstract

This paper presents the uncoupled anisotropic thermoelasticity formulation for Boundary Element using the Dual Reciprocity Method. The proposed formulation demands two interpolation functions, one for the particular solution of the elastic problem and another to interpolate the derivatives of the thermal field. For the last one, the shape parameter of multiquadric Radial Basis Functions (RBF) was demonstrated to be a problem to be solved. Some polyharmonic RBFs were tested and a modified function was proposed. Singularity at derivatives of some functions made them impossible to be used and some of them demonstrated great sensibility to saturation. The use of internal points was mandatory under temperature fields irregularly distributed. Few polyharmonic functions were suitable to be used for thermoelasticity and low order functions presented reliable solutions, with the proposed modified RBF presenting slightly better performance over the previously used function.

Keywords
Boundary Element; Dual Reciprocity; Radial Basis Functions; Thermoelasticity

Graphical Abstract

1 INTRODUCTION

Thermoelasticity theory examines the impact of a temperature field on an elastic body, resulting in strains and stresses. This theory finds application in various engineering scenarios, including environmental exposure of buildings, computer processor chips, airplane structures and turbines, components of automotive engines, residual stresses from thermal treatments or cooling processes in cast iron, and the behavior of 3D printed materials, among other instances. Stresses and strains produced by thermal fields may be significant and produce early failure if not considered.

The Boundary Element Method (BEM) has been extensively applied to thermoelasticity problems. If a steady state temperature field is considered, the potential field can be solved independently of the elastic problem, which is a class of problems called uncoupled thermoelasticity. After that, the thermal effects are included in elasticity and they appear at the Boundary Integral Equation (BIE) as a domain integral, like a body force. To solve it without discretizing the domain, Rizzo and Shippy (1977) presented a simplification able to transform the volume integral into a surface one, but it demands that the second order derivative of the thermal effects should be a constant value, limiting its application.

Gao (2003) applied the Radial Integration Method (RIM) for 2D and 3D isotropic thermoelasticity. For these cases, a thermal field was imposed as a known function and the thermoelasticity BIE was easily solved. However, if the thermal field is an unknown function, the RIM is more complicated to use and may be unattractive.

Another possibility is the use of the Dual Reciprocity Method (DRM) proposed by Brebbia and Nardini (1983), which consists of approximating the kernel of the domain integral by means of the Radial Basis Functions (RBF). This makes this technique capable of solving any domain integral, being a generic solution technique, but it demands the inversion of a fully populated interpolant matrix (Goldberg et al. 1999). Kögl and Gaul (2003) applied it to thermoelasticity and demonstrated that the gradient field that appears in volume integrals can be approximated by the derivative of interpolation functions for the temperature field.

Theoretically, DRM can approximate any kernel of a domain integral and solve it, but the quality of the RBF as interpolant is an important matter, and some kind of functions are preferred to some problems. For example, the thin plate spline of Table 1 is suitable for two-dimensional problems, while the linear function is adequate to three-dimensional ones (Goldberg et al., 1998). The polyharmonic spline functions (PHS) are the simplest RBF and depend only on the Euclidean distance, but the interpolation matrix may become ill conditioned as points used for interpolation become close to each other (Goldberg et al. 1999). Also, not all integer values of $m$ can be successfully applied, as indicated in Table 1 (Goldberg et al. 1998).

Table 1
Common RBF types.

Piecewise smooth RBF’s
Type of function	Radial function	comments
Polyharmonic spline (PHS)	$r^m$	$m=\mathrm{1,3},5,\dots$
	$r^m\ln \left(r\right)$	$m=\mathrm{2,4},6,\dots$
Linear	$r$	particular PHS case
Cubic	$r^3$	particular PHS case
Thin plate spline (TPS)	$r^2\ln \left(r\right)$	particular PHS case
Infinitely smooth RBF’s
Type of function	Radial function	comments
Gaussian (GA)	$e^{-{\left(ϵr\right)}^2}$	$ϵ>0$
Multiquadric (MQ)	$\sqrt{1+{\left(ϵr\right)}^2}$	$ϵ>0$
Inverse Quadratic (IQ)	${\left(1+{\left(ϵr\right)}^2\right)}^{-1}$	$ϵ>0$
Inverse Multiquadric (IMQ)	${\left(\sqrt{1+{\left(ϵr\right)}^2}\right)}^{-1}$	$ϵ>0$

Karur and Ramachandran (1994) examined the convergence on bidimensional domains of a linear, a thin plate spline (TPS) and a linear with a constant and estimated the condition number of the matrix generated by interpolation, although no relation they found among them. Goldberg et al. (1996) highlighted that TPS is better interpolation than cubic, but since they have linear convergence rates, they studied the multiquadric (MQ) first studied by Franke (1982), which has exponential convergence. When $ϵ$increases from zero, it increases the accuracy until an optimal value, and after that, the accuracy decreases fast. One of the reasons pointed out by the authors for this is that as $ϵ$ becomes larger, the interpolation matrix becomes ill-conditioned. Goldberg et al. (1996) applied cross validation to calculate a good shape parameter. Carlson and Foley (1991) found out that the magnitude of the body force to be interpolated strongly influences the shape parameter value, while the number of data points and their location had little influence. Another variation of this problem was proposed by Kansa (1990) with a scheme using variable shape parameters.

Not only the MQ, but all infinitely smooth RBF in Table 1 are dependent on the shape parameter. A list of some recent works on finding the optimal shape parameter is given in Pooladi and Larsson (2023), but independently of the technique used, the problem of calculating an optimal shape parameter means extra computational cost, and depending on the size of the problem, it can significantly increase the computation time.

The use of interpolation functions to solve differential equations is an important matter in the Method of Fundamental Solutions (MFS) and other meshless methods and many developments have been done. Tolstykh and Shirobokov (2003) proposed that interpolation should use only some points around the point of interest, instead of using the known values on all domain of the problem. In other words, a local approximation is proposed instead of a global approximation, generating a sparse interpolation matrix, making it easier to invert. This method was named RBF-FD and can use any of the functions in Table 1. Flyer et al. (2016a) experimented to augment the RBF functions with polynomials and noticed that for PHS, the high condition number of interpolation matrix was not a problem in accuracy and the use of high order polynomials improved accuracy until a certain level. Flyer et al. (2016b) found out that this limit on the polynomial degree depends on the domain dimension and the number of nodes used in interpolation.

The augmented polynomials RBF are mentioned many times in papers, as Goldberg et al. (1999), but only first order polynomials are cited. Application of the limit found by Flyer et al. (2016b) would greatly increase the size of the interpolation matrix, greatly improving the computational cost for matrix inversion.

Using BEM and DRM to solve decoupled thermoelasticity as presented by Kögl and Gaul (2003) demands two interpolation functions: one for the particular solution of elasticity and another for interpolating the temperature field to approximate its gradient. For the second function, Kögl and Gaul (2003) only present Function A of Table 2 without explaining the reasons for this choice. After two decades of this publication, the papers usually cite the precursor work and their successful application of DRM to thermoelasticity, and, to the authors knowledge, only Benedetti (2023) has applied their formulation to study polycrystalline materials, using the same formulation and functions.

Table 2
Selected Radial Basis Functions.

Piecewise smooth RBF’s
Identification	Function $\mathit{\psi }$	Derivative ${\mathit{\psi }}_{,\mathit{j}}$
Function A	$1+r^2+r^3$	$\left(2+3r\right)r_j$
Function B	$1+r^3$	$3r{r}_j$
Function C	$1+r^5$	$5r^3{r}_j$
Function D	$1+r^2\mathrm{l}\mathrm{n}\left(r+1\right)$	$\left[2\ln \left(r+1\right)+\frac{r}{r+1}\right]r_j$
Function E	${1+r}^3\mathrm{l}\mathrm{n}\left(r+1\right)$	$\left[3r\ln \left(r+1\right)+\frac{r^2}{r+1}\right]r_j$

As it has been presented in the last paragraphs, there are many studies concerning different RBF and different applications, and how they affect the results, but no study concerning different interpolation functions for the derivative of the temperature field at different conditions. The aim of this work is to investigate how some of the RBF listed in Table 1 affect the results of uncoupled thermoelasticity solution considering isotropic and anisotropic three-dimensional problems, making use of the generality formulation of DRM which is able to solve any source term (Partridge, Brebbia and Wrobel, 1992). The goal is to find at least one simple RBF function that can be used to solve any thermoelasticity problem.

This paper is organized as follows: Section 2 presents the thermoelasticity formulation and the use of the Dual Reciprocity Method to solve the volume integrals with thermoeffects. Section 3 presents the problems used for numerical simulations and the reference results for comparison, as well as the Radial Basis Functions used and the reasons for their choice. Section 4 shows the results of the numerical simulations and discusses them. Section 5 concludes the paper.

2 THERMOELASTICITY

2.1 Boundary Integral Equation

For uncoupled thermoelasticity, the thermal problem can be solved first and then, the temperature field effects can be applied to elasticity (Kögl and Gaul, 2003). In this analysis, the potential field is solved using the standard isotropic formulation given by Brebbia and Dominguez (1992). Anisotropic potential problems are solved using the Direct Domain Mapping technique presented by Shiah and Tan (2004), which transforms any anisotropic geometry into an equivalent isotropic geometry, by deforming it. The potential field is solved in this new deformed and isotropic domain, and then, the results are transported back to the original anisotropic geometry.

Once the temperature field is solved, thermal effects are added to the elasticity problem and the equilibrium equation becomes as is given in Equation 1.

${\sigma }_{ij}=C_{ijkl}{ϵ}_{kl}-{\gamma }_{ij}\theta$,(1)

where the subscripts mean summation, $C_{ijkl}$ is the stiffness tensor, ${\gamma }_{ij}$ is the thermoelastic tensor calculated by Equation 2and $\theta$ is the temperature difference from the reference temperature $T_0$ and temperature analysis $T$:

${\gamma }_{ij}\mathrm{}=\mathrm{}{C}_{ijkl}{\alpha }_{kl}$,(2)

where ${\alpha }_{kl}$ is the thermal expansion tensor. Since deformations are small, they can be expressed in terms of the displacement derivatives as

${ϵ}_{kl}=\frac{1}{2}\left(u_{k,l}+u_{l,k}\right)$.(3)

The equilibrium equation is given by:

${\sigma }_{ji,j}+\rho b_i=0$ .(4)

Applying Equations 3 and 4 into 1, a partial differential equation appears:

$C_{ijkl}{u}_{k,lj}={\gamma }_{ij}{\theta }_{,j}-\rho b_i$ (5)

Right hand side of Equation 5 represents the body forces, where the first is a thermal effect and the second an elastic effect. Defining the following differential operator:

${\mathcal{L}}_{ik}=C_{ijkl}\frac{{\partial }^2}{\partial x_l\partial x_j}$,(6)

and using it into Equation 5, one obtains

${\mathcal{L}}_{ik}{u}_k={\gamma }_{ij}{\theta }_{,j}-\rho b_i$.(7)

Equation 7 can be integrated using an interpolation function $u_{mi}^{*}$ which is the elasticity fundamental solution. Integrating by parts, using Betti’s reciprocal theorem and Cauchy’s Principal Value definition, the inverse integral statement is obtained in Equation 8.

${\delta }_{mi}c\left(\xi \right)u_i\left(\xi \right)+{\int }_{\mathrm{\Gamma }}^{\mathrm{}}{t}_{mi}^{\mathrm{*}}{u}_{i}d\mathrm{\Gamma }={\int }_{\mathrm{\Gamma }}^{\mathrm{}}{u}_{mi}^{\mathrm{*}}{t}_{i}d\mathrm{\Gamma }+{\int }_{\mathrm{\Gamma }}^{\mathrm{}}{u}_{mi}^{\mathrm{*}}{\gamma }_{ij}{n}_j\theta d\mathrm{\Gamma }-{\int }_{\mathrm{\Omega }}^{\mathrm{}}\left({\gamma }_{ij}{\theta }_{,j}-\rho b_i\right)u_{mi}^{\mathrm{*}}d\mathrm{\Omega }$ (8)

where $t_{mi}^{*}$ is the derivative of $u_{mi}^{*}$ given by:

$t_{mi}^{\mathrm{*}}=C_{ijkl}{u}_{mi,j}^{\mathrm{*}}{n}_l$,(9)

and $c\left(\xi \right)$ is 0.5 for smooth surfaces. The second boundary integral on right hand side of Equation 8 is the effect of temperature and it is easily solved since it is already a boundary integral. The gradient thermal field effect to elasticity is inside the domain integral and is considered a body force of thermal effect.

2.2 Dual Reciprocity Method

The expression in parenthesis in Equation 8 is approximated by a sum of Radial Basis Functions:

$f_i=-\rho b_i+{\gamma }_{ij}{\theta }_{,j}\mathrm{}\approx \mathrm{}{\sum }_{q=1}^N{\varphi }_{in}^q{\alpha }_n^q$,(10)

with ${\varphi }_{in}^q={\varphi }_{in}^q\left(r\right)$, where $r$ is the Euclidean distance. Equation (10) can be written in matrix form as:

$\left\{\mathit{f}\right\}=\left[\mathit{F}\right]\left\{\mathit{\alpha }\right\}$,(11)

which can be solved for $\left\{\mathit{\alpha }\right\}$. Substituting Equation (10) into the domain integral of Equation (8), using the differential operator defined in Equation (6), and integrating by parts, the domain integral is transformed into boundary integrals given by Equation (12) (for more details, see Brebbia and Dominguez (1992) or Kögl and Gaul (2003)):

$I={\int }_{\mathrm{\Omega }}^{\mathrm{}}\left(-\rho b_i+{\gamma }_{ij}{\theta }_{,j}\right)u_{mi}^{\mathrm{*}}d\mathrm{\Omega }\mathrm{}~{\sum }_{q=1}^N\mathrm{}\mathrm{}\left\{{\int }_{\mathrm{\Gamma }}^{\mathrm{}}{u}_{mi}^{\mathrm{*}}{t}_{in}^{q}d\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}^{\mathrm{}}{t}_{mi}^{\mathrm{*}}{u}_{in}^{q}d\mathrm{\Gamma }-{\delta }_{mi}c\left(\xi \right)u_{in}^q\left(\xi \right)\right\}{\alpha }_n^q$,(12)

where

$t_{in}^q=C_{ijkl}{u}_{kn,j}^q{n}_l$.(13)

To obtain the coefficients ${\alpha }_n^q$ used in Equation (12), Kögl and Gaul (2003) consider that the temperature gradient ${\theta }_{,j}$is not a known function at the entire domain, but it is only known at some points. To evaluate this term, the temperature field is interpolated by a sum of RBF’s given by:

$\theta =\mathrm{}{\sum }_{q=1}^N{\psi }^q{\zeta }^q$,(14)

where ${\psi }^q$ are the temperature interpolation functions. Equation (14) generates the matrix system of Equation (15), which can be inverted to calculate coefficients ${\zeta }^q$.

$\left\{\mathit{\theta }\right\}=\left[\mathit{E}\right]\left\{\mathit{\zeta }\right\}$ (15)

Now, Equation 14 is derived with respect to cartesian coordinates:

${\theta }_{,j}=\mathrm{}{\sum }_{q=1}^N{\psi }_{,j}^q{\zeta }^q$ (16)

With the coefficients ${\zeta }^q$ calculated previously, the temperature gradient ${\theta }_{,j}$ is interpolated by Equation 16 and can be applied in Equation 10 to solve the coefficients ${\alpha }_n^q$. So, this approach presented by Kögl and Gaul (2003) needs to use two RBFs: the first is the particular solution ${\varphi }_{in}^q$ of DRM, and the second is the interpolation function for temperature field ${\psi }^q$ and its derivative ${\psi }_{,j}^q$.

It is a natural expectation that the previous studies presented in Section 1 apply to functions ${\psi }^q$. Nevertheless, they are used to interpolate a derivative of potential field by using the potential field, in a way that the derivative of the RBF plays the main role in producing the derivative of the potential field, and some specific situations may happen. For example, Mai-Duy and Tran-Cong (2003) showed that the use of a good interpolation function for the primitive does not guarantee the derivative will also be a good interpolant. Moreover, Kögl and Gaul (2003) presented this formulation and chose the interpolant function although no reason for the choice of the function had been presented.

In this paper, thermoelasticity was implemented in BESLE, an opensource MPI parallelized software developed by Galvis et al. (2021) that uses linear triangular elements, elasticity fundamental solutions of Barnett-Lothe tensor expressed by Fourier series (Tan et al., 2013) and particular solution ${\varphi }_{in}^q$ and derivatives expressed in Equations 17 to 19.

${\varphi }_{in}^q=\mathrm{}{\delta }_{in}^q\left(r^2+r^3\right)$ (17)

${\varphi }_{in,l}^q=\mathrm{}{\delta }_{in}^q\left(2r+{3r}^2\right)r_{,l}$ (18)

${\varphi }_{in,lm}^q=\mathrm{}{\delta }_{in}^q\left[\left(2+3r\right){\delta }_{lm}+3r{r}_{,l}{r}_{,m}\right]$ (19)

3 NUMERICAL EXPERIMENTS

3.1 Problems

The numerical experiments were done using several RBFs in three different problems. The first one was presented and solved by Gao (2003) and consists of an isotropic plate with imposed temperature given by Equation 20, with boundary conditions and material properties given in Figure 1. The plate width is normal to the figure and is equal to 1. Faces normal to $z$ axis have restricted displacements at $z$ direction. The analytical solution is given in Equation 21. Meshes of 80, 180, 320, 500, 720, 980, 1280, 1620, and 2000 elements were used, combined with the use of 0, 1, 3, 20, 63, 144, 275, 468, 735, 1088, 1539, 2100, and 2783 internal points.

Figure 1
Problem 1 - isotropic plate subjected to quadratic thermal condition.

$\theta \left(y\right)=40y^2+60y$ (20)

$u_y=k\left(\frac{1+\nu }{1-\nu }\right)\left(\frac{40}{3}{y}^3-30y^2+\frac{55}{6}\right)$ (21)

The second problem is an orthotropic unit cube with one fixed face and all the others free to deform. The fixed face has a thermal boundary condition $\theta =0$ while its opposite face is at $\theta =100$, and the other four faces are adiabatic, as shown in Figure 2. The displacement $u_x$, $u_y$ and $u_z$ are calculated and compared with ANSYS solutions using 1728 quadratic cubic elements. Tensors of thermal conductivity, expansion, and elasticity are given in Equations 22 to 24. Note that the stiffness matrix is given in terms of Voigt’s notation. This problem was solved with meshes of 48, 108, 192, 432, 768, 972, and 1758 elements using 0, 1, 8, 27, 64, 125, 216, 343, 512, 729, 1000, 1331, 1728, and 2197 internal points.

Figure 2
Problem 2 - orthotropic cube subjected to variable temperature.

$\kappa =\left(\begin{array}{ccc}1& 0& 0\\ 0& 2& 0\\ 0& 0& 3\end{array}\right)\mathrm{}W/m\mathrm{º}C$ (22)

$\alpha =\left(\begin{array}{ccc}2& 0& 0\\ 0& 2.5& 0\\ 0& 0& 3\end{array}\right)\mathrm{}\times {10}^{-5}\mathrm{}\mathrm{º}{C}^{-1}$ (23)

$C=\left(\begin{array}{c}\begin{array}{cccccc}1334.5& 738.8& 619& 0& 0& 0\\ 738.8& 1701.2& 646.9& 0& 0& 0\\ 619& 646.9& 1147.4& 0& 0& 0\\ 0& 0& 0& 800& 0& 0\\ 0& 0& 0& 0& 1200& 0\\ 0& 0& 0& 0& 0& 1000\end{array}\end{array}\right)\mathrm{}Pa$ (24)

The third problem was solved by Kögl and Gaul (2003) and consists of the L shape geometry shown in Figure 3 with anisotropic material which elastic and thermoelastic tensors are given in Equations 25 and 26, respectively. Thermal and mechanical boundary conditions are also given in Figure 3. All other faces are elastically free and insulated. Kögl and Gaul (2003) presented only $u_z$ displacements at the line $\left(50,0,z\right)$ which are used as reference. To be compatible with the reference solution mesh, few meshes could be used and they have 28, 112, 448, 700, and 2800 elements using 0, 3, 32, 99, 224, 438, 720, 1127, 1664, 2474, 3200, 4235, 5472, and 6929 internal points.

Figure 3
Problem 3 – anisotropic material subjected to bidimensional temperature field.

$C=\left(\begin{array}{c}\begin{array}{cccccc}430.1& 130.4& 18.2& 0& 0& 201.3\\ 130.4& 116.7& 21.0& 0& 0& 70.1\\ 18.2& 21.0& 73.6& 0& 0& 2.4\\ 0& 0& 0& 19.8& -8.0& 0\\ 0& 0& 0& -8.0& 29.1& 0\\ 201.3& 70.1& 2.4& 0& 0& 147.3\end{array}\end{array}\right)\mathrm{}GPa$ (25)

$\gamma =\left(\begin{array}{c}\begin{array}{ccc}1.01& 2.00& 0\\ 2.00& 1.48& 0\\ 0& 0& 7.52\end{array}\end{array}\right)\mathrm{}{10}^6\mathrm{}N/K.m^2$ (26)

For each problem listed, the displacements are calculated by combining the different refined meshes with different amounts of internal points used for the DRM procedure. The surface mesh is regular, i.e., all the elements on the boundary surface are of the same size. Also, internal points, when used, are regularly spaced from each other and the boundary.

3.2 Temperature Interpolation Functions

Before applying the RBF for temperature interpolation on problems presented previously, some considerations are made about the derivatives of the temperature field calculated by Equation 16 using the derivatives of the interpolation functions ${\psi }_{,j}^q$. The derivative of the first RBF in Table 1 is given by:

${\psi }_{,j}=\mathrm{}{\left(r^m\right)}_{,j}=m{r}^{m-1}{r}_{,j}=m{r}^{m-1}\left(\frac{r_j}{r}\right)=m{r}^{m-2}{r}_j$,(27)

with $r_j$ being the difference of cartesian coordinates in $j$ direction. So, for $m=1$, the derivative gets a singularity when $r$ goes to $0$, which makes it unable to be used in Equation 16. For $m>1$, despite no singularity, the literature indicates that only odd values of $m$ should be used. In fact, Figure 4 shows the results obtained for the orthotropic cube using the Kögl and Gaul (2003) temperature interpolation function, $m=2$ and $m=4$. Different meshes were used, all of them with results like these. Also, the displacements have order of ${10}^{28}$, being obviously nonrepresentative of the solution. Even coefficients demonstrated themselves useless for interpolation purpose.

Figure 4
Deformed geometry obtained for Problem 2 using different RBF for temperature interpolation. Left: $\psi =\mathrm{}1+r^2+r^3$; center: $\psi =\mathrm{}1+r^2$; right: $\psi =\mathrm{}1+r^4$.

The singularity in derivatives also appears when using the second RBF function presented in Table 1:

${\psi }_{,j}=\mathrm{}{\left[r^m\ln \left(r\right)\right]}_{,j}=m{r}^{m-1}\ln \left(r\right)r_{,j}\mathrm{}+\mathrm{}{r}^m\left(\frac{1}{r}\right)r_{,j}\mathrm{}=\left[m\ln \left(r\right)+1\right]r^{m-2}{r}_j$.(28)

It is clear from Equation 28 that singularity will happen when $r$ goes to $0$ independently of the value of $m$. This singularity can be overcome if a small change is made in the function, as shown in Equation 29, with derivatives given in Equation 30.

$\psi =r^m\ln \left(r+1\right)$ (29)

${\psi }_{,j}=\mathrm{}{\left[r^m\ln \left(r+1\right)\right]}_{,j}=m{r}^{m-1}\ln \left(r+1\right)r_{,j}\mathrm{}+\mathrm{}{r}^m\frac{1}{r+1}{r}_{,j}\mathrm{}=\left[m{r}^{m-2}\ln \left(r+1\right)+\frac{r^{m-1}}{r+1}\right]r_j$ (30)

Now, singularity is avoided for $m>1$. For the infinitely smooth functions of Table 1, the presence of the shape parameter is a considerable problem to be solved. In an attempt to use the MQ function, a range of $ϵ$ values from near zero to greater than ${10}^4$ were tried on problem 1, but the results were unable to minimally represent the analytical solution, demonstrating how much the results are dependent on the shape parameter and the importance of calculating it. Because of the extra computational cost necessary to find the optimal $ϵ$, the infinitely smooth functions do not make part of the tests done in this paper.

With these considerations, the functions shown in Table 2 and their derivatives were selected and applied to the problems presented previously. All the functions were added by 1 for completeness (Karur and Ramachandran, 1994). Each of them was applied to the problems presented in Section 3.1. The first function of Table 2 was used by Kögl and Gaul (2003), and it is used as a benchmark.

4 RESULTS

4.1 Function A

This section presents the results obtained using Function A. Figure 5 is a log-log graph of the L₂ error calculated for displacement $u_y$ at the first problem. The first point of each series data are the results obtained without using internal points, and each series is defined by the number of elements of the mesh. The results show clearly that mesh refinement was the main mechanism for improving accuracy to a certain level. The L₂ value achieved without using internal points of the mesh with 720 is the lower, demonstrating that convergence reached a maximum. This phenomenon is called saturation, when the increasing of points used for interpolation does not help in improving accuracy.

Figure 5
L₂ error for displacements $u_y$ for problem 1 using Function A.

For course meshes, the use of internal points resulted in the loss of accuracy. For more refined meshes, a small improvement in results can be seen with the increasing number of internal points and then the accuracy drops, showing that saturation has been achieved. The more refined a mesh is, more internal points are needed to reach the optimal condition. Besides this isotropic problem has a quadratic temperature profile given in Equation 20, the use of internal points showed a small improvement in displacement results.

Figure 6 shows the average relative difference to ANSYS calculated for displacement $u_x$ for the orthotropic cube. Displacements $u_y$ and $u_z$ are not presented because they showed no significant difference to the $u_x$ displacements. As the mesh is refined, the results approach the reference and the use of internal points improved them to a certain level, with asymptotic behavior. Different from the isotropic problem, using internal points only improved the accuracy of the results, with no concern in finding an optimal condition.

Figure 6
Relative diference for displacements $u_x$ for problem 2 using Function A.

Figure 7 shows the relative difference for displacements $u_z$ for the third problem. In this anisotropic condition, mesh refinement alone was unable to give satisfactory accuracy, and using internal points was essential to achieve it. Also, mesh refinement alone produced results even worse. It is possible to see the saturation is reached when internal points are used and that this happens with a number of internal points that is dependent on the mesh size. The best result of the mesh with 28 elements was achieved using 99 internal points, while the best result of the meshes with 112 and 448 elements were obtained using 3200 internal points. For the meshes with 700 and 2800 elements, the number of internal points tested was not enough to achieve saturation.

Figure 7
Relative diference for displacements $u_z$ for problem 3 along line $\left(50,\mathrm{}0,\mathrm{}z\right)$ using Function A.

4.2 Function B

As commented in Section 3.2, the function $\psi =1+r^2$ does not produce good results, but the even coefficient is presented in Function A. So, Function B is obtained by deleting it from Function A, and it is expected to achieve better accuracy by doing this. The results obtained for the problems are given in Figures 8, 9, and 10.

Figure 8
L₂ error for displacements $u_y$ for problem 1 using Function B.

Figure 9
Relative diference for displacements $u_x$ for problem 2 using Function B.

Figure 10
Relative diference for displacements $u_z$ for problem 3 along line $\left(50,\mathrm{}0,\mathrm{}z\right)$ using Function B.

Surprisingly, the results have no significant difference from the ones obtained using Function A. So, there is no advantage in using a more complete radial basis, but also no disadvantage in using the combined terms with exponent 2 and 3.

4.3 Function C

Function C is an attempt to use a high order RBF function. Figure 11 shows the results obtained in the first problem. As the mesh is refined, the accuracy of results increases, but at a certain degree of refinement they diverge, and the use of internal points doesn’t make any improvement at the solution. At the mesh with 1280 elements, despite the good results without using internal points, bad results appear suddenly when using 735 internal points. For the meshes with 1620 and 2000 elements, saturation is achieved, and Function C behaves in an unpredictable way.

Figure 11
L₂ error for displacements $u_y$ for problem 1 using Function C.

This bad behavior repeats at orthotropic and anisotropic problems, as shown in Figures 12 and 13. After saturation is reached, the results are no longer representative of the solution. In Figure 12, the mesh with 768 elements had a suddenly bad result using 1331 points, showing that the distance between the interpolation points is not the only factor for this function to produce good or bad results.

Figure 12
Relative diference for displacements $u_x$ for problem 2 using Function C.

Figure 13
Relative diference for displacements $u_z$ for problem 3 along line $\left(50,\mathrm{}0,\mathrm{}z\right)$ using Function C.

Figure 14 illustrates the effect of saturation. The left cube is a good result obtained with a mesh of 972 elements and 1331 internal points. When using 1728 or 2197 internal points, the graphical difference is difficult to be seen in this kind of view, because the relative difference lies below 3%. For this reason, it is presented the deformation of the mesh of 1758 elements without internal points at the center of the Figure 14, which starts to lose symmetry despite of keeping the general shape. Using only 1 internal point is enough to change the shape and destroy the symmetry.

Figure 14
Results for displacements$u_x$. Left: mesh of 972 elements and 1331 internal points; center: mesh of 1758 elements without internal points; right: mesh of 1758 elements and 1 internal point.

4.4 Function D

Figure 15 shows the average error of $u_y$ displacements by using Function D to solve the first problem and few differences appear in the results of Function A. The performance of Function D without using internal points was a little better than Function A after some degree of refinement. Both functions experiment stabilization of accuracy with mesh refinement, but Function D shows slightly better results. Concerning the use of internal points, both functions improved results before saturation with refined mesh, and then the accuracy was reduced. It is interesting to note that errors for using internal points after saturation are higher than the ones with no internal points.

Figure 15
L₂ error for displacements $u_y$ for problem 1 using Function D.

For the second problem, as shown in Figure 16, Function D also performed better than Function A when no internal points are used and both functions have the tendency of very slow increase accuracy with the increasing number of internal points, at all meshes, although this better accuracy is not significative for refined meshes.

Figure 16
Relative diference for displacements $u_x$ for problem 2 using Function D.

Figure 17 shows the relative difference in displacements for the third problem. Like Function A, the use of internal points was also mandatory to achieve good results, but the accuracy of Function D is a little slower specially for the two most refined meshes, presenting similar accuracies if huge amounts of internal points are used.

Figure 17
Relative diference for displacements $u_z$ for problem 3 along line $\left(50,\mathrm{}0,\mathrm{}z\right)$ using Function D.

4.5 Function E

The saturation problem with Function C also happens with Function E. Figure 18 shows the L₂ error for displacements $u_y$ calculated at the isotropic problem. The mesh with 1280 elements starts diverging after using 1088 internal points and meshes of 1620 and 2000 show the loss of accuracy in results caused by saturation. The results of Functions E for problems 2 and 3 are given in figures 19 and 20. It can be seen at these cases that the saturation only appears at the most refined meshes, showing that Function C is more sensible to it than Function E.

Figure 18
L₂ error for displacements $u_y$ for problem 1 using Function E.

Figure 19
Relative diference for displacements $u_x$ for problem 2 using Function D.

Figure 20
Relative diference for displacements $u_z$ for problem 3 along line $\left(50,\mathrm{}0,\mathrm{}z\right)$ using Function D.

4.6 Comments on the results

One point to be observed is the necessity of internal points at the anisotropic problem for good results. The internal points are used to interpolate a body force given by Equation 10, composed of an elastic force and a thermal effect. Along all problems in this paper, no elastic force has been introduced, so only the dilation is responsible for the results. Figure 21 presents the isothermals which are a fraction of a tenth of the range of temperatures of each problem, with minimum temperature at blue colors, and maximum at red colors.

Figure 21
Temperature distribution and isothermals for each problem.

For the first two problems, the temperature distribution is well regular, but for the third problem, there is almost two-thirds of the figure with only 20% of the range of temperature, with the other third with 80% of the range. This is a very irregular distribution, and the use of boundary nodes only was insufficient to make a good representation of the domain temperature field by the interpolation functions. In this situation, the internal points were essential to represent the temperature field and its thermal effect for the thermoelasticity.

Another important matter is the saturation presented by Functions C and E at refined meshes or increasing the number of internal points. Since Functions A and B presented the same results, the second one will be used instead at this analysis. Figure 22 presents the graph of the last four RBF of Table 2, and it can be seen that Functions C and E (blue and orange lines respectively) keeps their value $\varphi \left(r\right)$ close to 1 until $r$ near 0.4, while Functions B and D (red and green lines) start to move away from 1 at values lower than $r=0.2$. As the mesh refinement or the number of internal points increases, the distance between them reduces and the values of the functions come to these regions. So, Functions C and E generate lines at the interpolation matrix similar to each other, affecting the evaluation of the interpolation coefficients. This also means that Functions B and D can also present stagnation, but the meshes used in this paper was not enough to generate them.

Figure 22
RBF functions at the interval $0<r<1$.

5 CONCLUSIONS

In this paper, the uncoupled thermoelasticity formulation for Boundary Element using Dual Reciprocity Method was presented and interpolation functions were used to approximate the gradient of the thermal field. Singularity needs to be avoided in these approximations and for this reason, a modification in polyharmonic spline using natural logarithm generally used at the Dual Reciprocity method was proposed. It was demonstrated that the choice of temperature interpolation function is very important to guarantee a reliable solution.

The multiquadric RBF was tested with a great range of values for the shape parameter but no good result was found. An optimization procedure for this parameter should be implemented to verify the performance of this function since the manual search was not capable of finding a good one.

The polyharmonic RBF of the form $r^m$ with $m$ even didn’t work out, showing unreal displacements, which agrees with limitations on the use of these functions presented by other authors (see, e.g., Bayona 2019). For odd coefficients, the singularity of $m=1$ and sensibility to saturation of $m=5$ proved that these coefficients are not a good choice. The only coefficient that produced good results was $m=3$, the cubic function.

For the polyharmonic RBF with logarithmic expression, the singularity when $r$ tends to $0$ in derivatives made them unable to be used. The proposed function in Equation 29 overcomes the singularity but only produced good results with $m=2$. For $m=3$, the saturation with mesh refinement makes the solution diverge and produce not representative results.

It was shown that high order functions are more susceptible to saturation than low order functions with the reduction of the distance between the interpolation points (mesh refinement or increasing the number of internal points), which makes similar for a range of close points. So, high order functions should be avoided.

For problems with regular temperature field, the results showed that the use of internal points is not necessary and refinement of the mesh is a better way to improve the solution. For irregular temperature fields, the use of internal points was mandatory, demanding mesh refinement and a great number of internal points to produce good results. The increasing number of internal points improved the solution until saturation, but this phenomenon was harder to achieve with the mesh refinement.

The benchmark function used by Kögl and Gaul (2003) (Function A), the cubic function (Functions B), and the modified proposed function (Function D) produced similar results, with the last one obtaining better results for problems without needing internal points, and better stability over saturation conditions.

Finally, Function B is favored over Function A because the presence of the term $r^2$ in Function B signifies only computational cost without accompanying benefits. Given the superior performance of the proposed modified PHS Function D, which operates without relying on internal points and demonstrates stability under saturation, it is advisable to choose it over the benchmark function.

Acknowledgments

The authors thank the University of Campinas (UNICAMP) for the facilities and structures provided for this work, and the financial support received by the following organizations: Capes (23038.005057/2015-68), Air Force Office of scientific Research (AFOSR) of United States (FA9550-18-1-0113 and FA9550-20-1-0133), University of Brasilia (UnB) and Instituto Federal do Espírito Santo (IFES).

References

• Bayona, V. (2019). An insight into RBF-FD approximations augmented with polynomials. Computers and Mathematics with Applications, 77:2337-2353.
• Benedetti, I. (2023). An integral framework for computational thermos-elastic homogenization of polycrystalline materials. Computer Methods in Applied Mechanics and Engineering.
• Brebbia, C.A. and Dominguez, J. (1992). Boundary Elements: An Introductory Course, WIT Press.
• Brebbia, C.A. and Nardini, D. (1983). Dynamic analysis in solid mechanics by an alternative boundary element procedure. Soil Dynamics and Earthquake Engineering 65:147-164.
• Carlson, R.E. and Foley, T.A. (1991). The parameter R2 in multiquadric interpolation. Computers & Mathematics Applications 21:29-42.
• Flyer, N., Barnett, G.A. and Wicker, L.J. (2016b). Enhancing finite differences with radial basis functions: Experiments on the Navier-Stokes equations. Journal of Computational Physics 316:39-62.
• Flyer, N., Fornberg, B., Bayona, V., Barnett, G.A. (2016a). On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy. Journal of Computational Physics 321:21-38.
• Franke, R. (1982). Scattered data interpolation: tests of some methods. Mathematics of Computations 48:181-200.
• Galvis, A.F., Prada, D.M., Moura, L.S., Zavaglia, C., Foster, J.M., Sollero, P., Wrobel, L.C. (2021). BESLE: Boundary element software for 3D linear elasticity. Computer Physics Communications 265:108009.
• Gao, X.W. (2003). Boundary element analysis in thermoelasticity with and without internal cells. International Journal for Numerical Methods in Engineering, 57:975-990.
• Goldberg, M.A., Chen, C.S. and Bowman, H. (1999). Some recent results and proposals for the use of radial basis functions in the BEM. Engineering Analysis with Boundary Elements, 23:285-296.
• Goldberg, M.A., Chen, C.S. and Karur, S.R. (1996). Improved multiquadric approximation for partial differential equations. Engineering Analysis with Boundary Elements 18:9-17.
• Goldberg, M.A., Chen, C.S., Bowman, H and Power, H. (1998). Some comments on the use of Radial Basis Functions in the Dual Reciprocity Method. Computational Mechanics, 21:141-148.
• Kansa, E.J. (1990). Multiquadrics – a scattered data approximation scheme with applications to computational fluid-dynamics – I. Computers & Mathematics Applications, 19:127-145.
• Karur, S.R. and Ramachandran, P.A. (1994). Radial Basis Function Approximation in the Dual Reciprocity Method. Mathematical and Computer Modelling, 20:59-70.
• Kögl, M. and Gaul, L. (2003). A boundary element method for anisotropic coupled thermoelasticity. Archive of Applied Mechanics, 73:377-398.
• Mai-Duy, N. and Tran-Cong, T. (2003). Approximation of function and its derivatives using radial basis function networks. Applied Mathematical Modelling, 27:197-220.
• Partridge, P.W., Brebbia, C.A. and Wrobel, L.C. (1992). The Dual Reciprocity Boundary Element Method, Computational Mechanics Publications.
• Pooladi, F. and Larsson, E. (2023). Stabilized interpolation using radial basis functions augmented with selected radial polynomials. Journal of Computational and Applied Mathematics 437:115482.
• Rizzo, F.J. and Shippy, D.J. (1977). An advanced boundary integral equation method for three-dimensional thermoelasticity. International Journal of Numerical Methods, 11:1753-1768.
• Shiah, Y.C. and Tan, C.L. (2004). BEM treatment of three-dimensional anisotropic field problems by direct domain mapping. Engineering Analysis with Boundary Elements 28:43-52.
• Tan, C.L., Shiah, Y.C. and Wang, C.Y. (2013). Boundary element elastic stress analysis of 3D generally anisotropic solids using fundamental solutions based on Fourier series. International Journal of Solids and Structures, 50:2701-2711.
• Tolstykh, A.I. and Shirobokov, D.A. (2003). On using radial basis functions in a “finite difference mode” with applications to elasticity problems. Computational Mechanics 33:68-79.

Edited by

Editor: Rogério José Marczak

Publication Dates

• Publication in this collection
  20 May 2024
• Date of issue
  2024

History

• Received
  04 Oct 2023
• Reviewed
  25 Nov 2023
• Accepted
  03 Mar 2024
• Published
  19 Mar 2024