Meshless Local Petrov-Galerkin Method for 3D Steady-State Heat Conduction Problems

Abstract

The Meshless Local Petrov-Galerkin (MLPG) method is applied for solving the three-dimensional steady state heat conduction problems. This method is a truly meshless approach; also neither the nodal connectivity nor the background mesh is required for solving the initial boundary-value problems. The penalty method is adopted to enforce the essential boundary conditions. The moving least squares (MLS) approximation is used for interpolation schemes and the Heviside step function is chosen for representing the test function. The numerical results are compared with the exact solutions of the problem and Finite Difference Method (FDM). This comparison illustrates the accuracy as well as the capability of this method.

1. Introduction

In the recent years, the meshless methods have evoked considerable interest among researches and engineers for finding solutions of initial- and boundary-value problems. The MLPG method is one of these meshless schemes. The main advantage of this method over other meshless methods is that this method is a “truly meshless” approach. In other words, no background mesh is required for the interpolation of the solution variables or for the evaluation of the various integrals appearing in the local weak formulation of the problem. The MLPG method has been demonstrated to be quite successful in solving various partial differential equations. The concept of this approach was first introduced by Atluri and Zhu [1]. They solved elastostatic problems in 2D domains. Lin and Atluri [2] introduced an upwinding scheme for analyzing steady-state convection-diffusion problems. Wu et al. [3] coupled the MLPG method with either the finite-element or the boundary-element method to enhance the efficiency of the MLPG method. Ching and Batra [4] augmented the polynomial basis functions with singular fields to determine deformations and stress fields near the crack tip for generally 2D mixed-mode problems. Gu and Liu [5] and Batra and Ching [6] applied the Newmark family of methods to analyze 2D transient elastodynamic problems. The bending of a thin plate has been studied by Gu and Liu [7] andLong and Atluri [8].

However, these works and several other published works in the literature were dealing with 2D boundary-value problems, and fewer researches were conducted on 3D problems. Han and Atluri [9] used the MLPG approach for 3D problems arising in elastostatics. They also applied the MLPG method for solving 3D elastic fracture [10] and 3D elastodynamic problems [11]. The main objective of this investigation is to develop the MLPG method for 3D steady-state heat conduction problems. First, the governing equations are given; then, the weak formulations of MLPG method and the MLS approximation are briefly introduced. Comparisons are made through some illustrative examples to show the accuracy of the present approach.

2. MLPG Formulation

In Figure 1, we illustrated the MLPG method. The problem domain is denoted by Ω, which is sounded by boundaries including essential boundary Γ₁, natural boundary Γ₂, and Robin's boundary Γ₃. In the MLPG method, the problem domain is represented by a set of arbitrarily distributed nodes, as shown in the figure. The weighted residual method is used to create the discrete system of equations.

Figure 1:

Subdomains and their boundaries.

The major idea of MLPG is that the implementation of the integral form of the weighted residual method is confined to a very small local subdomain of a node. This means that the weak form is satisfied at each node of the problem domain in a local integral sense. Therefore, the weak form is integrated over a local quadrature domain independent of other domains of other nodes. This becomes possible using the Petrov-Galerkin formulation, in which one has the freedom to choose the test and trial functions independently.

The heat conduction Poisson's equation with its corresponding boundary conditions can be written as

$\lambda \left(\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}+\frac{{\partial }^2\theta }{\partial z^2}\right)=\dot{\psi }\hspace{1em}\text{in}\hspace{0.17em}\Omega ,$

(1)

the essential (Dirichlet) boundary condition:

$\theta ={\theta }_1\hspace{1em}\mathrm{on}\hspace{0.17em}{\Gamma }_1,$

(2)

the natural (Neumann's) boundary condition:

$-\lambda \left(\frac{\partial \theta }{\partial x}{n}_x+\frac{\partial \theta }{\partial y}{n}_y+\frac{\partial \theta }{\partial z}{n}_z\right)=q\hspace{1em}\mathrm{on}\hspace{0.17em}{\Gamma }_2,$

(3)

and the Robin boundary condition:

$\lambda \left(\frac{\partial \theta }{\partial x}{n}_x+\frac{\partial \theta }{\partial y}{n}_y+\frac{\partial \theta }{\partial z}{n}_z\right)=h\left({\theta }_f-\theta \right)\hspace{1em}\mathrm{on}\hspace{0.17em}{\Gamma }_3,$

(4)

where θ represents temperature, λ the thermal conductivity, n _x , n _y , and n _z are the components of the outward unit vector to Γ, q the given heat flux, h the convection heat transfer coefficient, θ _f is the environmental temperature, $\dot{\psi }$ the heat source per unit mass, and Γ₁, Γ₂, and Γ₃ the boundaries at which the Dirichlet, Neumann, and Robin conditions apply, respectively.

In the Ω _s , the weighted integral form of (1) is given by

${{\displaystyle \int }}_{{\Omega }_s}^{}\left[\lambda \left(\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}+\frac{{\partial }^2\theta }{\partial z^2}\right)-\dot{\psi }\right]\upsilon d{\Omega }_s=0.$

(5)

Here, v is the test function, and Ω _s is the sub domain. To reduce this high-order differentiability requirement on θ, we can integrate (5) by parts. Using Gauss's theorem, we can obtain the following local weak formulation:

$\begin{array}{cc}\hfill & {{\displaystyle \int }}_{{\Omega }_s}^{}\left(\lambda \frac{\partial \theta }{\partial x}\frac{\partial \upsilon }{\partial x}+\lambda \frac{\partial \theta }{\partial y}\frac{\partial \upsilon }{\partial y}+\lambda \frac{\partial \theta }{\partial z}\frac{\partial \upsilon }{\partial z}+\dot{\psi }\upsilon \right)d{\Omega }_s\hfill \\ \hfill & \hspace{1em}\hspace{1em}-{{\displaystyle \int }}_{{\Gamma }_1}^{}\lambda \left(\frac{\partial \theta }{\partial x}{n}_x+\frac{\partial \theta }{\partial y}{n}_y+\frac{\partial \theta }{\partial z}{n}_z\right)\upsilon d{\Gamma }_1\hfill \\ \hfill & \hspace{1em}\hspace{1em}-{{\displaystyle \int }}_{L_s}^{}\lambda \left(\frac{\partial \theta }{\partial x}{n}_x+\frac{\partial \theta }{\partial y}{n}_y+\frac{\partial \theta }{\partial z}{n}_z\right)\upsilon d{L}_s=0.\hfill \end{array}$

(6)

In this equation, L _s is the other part of the local boundary which is inside the solution domain (Figure 1).

Substituting (3) and (4) into (6) leads to

$\begin{array}{cc}\hfill & {{\displaystyle \int }}_{{\Omega }_s}^{}\left(\lambda \frac{\partial \theta }{\partial x}\frac{\partial \upsilon }{\partial x}+\lambda \frac{\partial \theta }{\partial y}\frac{\partial \upsilon }{\partial y}+\lambda \frac{\partial \theta }{\partial z}\frac{\partial \upsilon }{\partial z}+\dot{\psi }\upsilon \right)d{\Omega }_s\hfill \\ \hfill & \hspace{1em}\hspace{1em}+{{\displaystyle \int }}_{{\Gamma }_2}^{}q\upsilon d{\Gamma }_2-{{\displaystyle \int }}_{{\Gamma }_3}^{}h\left({\theta }_f-\theta \right)\upsilon d{\Gamma }_3\hfill \\ \hfill & \hspace{1em}\hspace{1em}-{{\displaystyle \int }}_{{\Gamma }_1}^{}\lambda \left(\frac{\partial \theta }{\partial x}{n}_x+\frac{\partial \theta }{\partial y}{n}_y+\frac{\partial \theta }{\partial z}{n}_z\right)\upsilon d{\Gamma }_1\hfill \\ \hfill & \hspace{1em}\hspace{1em}-{{\displaystyle \int }}_{L_s}^{}\lambda \left(\frac{\partial \theta }{\partial x}{n}_x+\frac{\partial \theta }{\partial y}{n}_y+\frac{\partial \theta }{\partial z}{n}_z\right)\upsilon d{L}_s=0.\hfill \end{array}$

(7)

The MLS approximation function is given by

$\begin{array}{c}\theta ={\displaystyle \sum }_{I=1}^N{\varphi }_I\stackrel{\wedge }{{\theta }_I}\end{array}.$

(8)

Here, $\stackrel{\wedge }{\theta }$ is the unknown fictitious nodal values, and Ø _I is usually called as the shape function of MLS approximation corresponding to nodal point x _I (Figure 2).

Figure 2:

Schematics of the MLS approximation.

Substituting (8) into (7) for all the nodes and applying the penalty factor method for imposing the essential boundary conditions yield the following linear equation:

$\begin{array}{cc}\hfill & {\displaystyle \sum }_{J=1}^M{{\displaystyle \int }}_{{\Omega }_s}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}\frac{\partial {\upsilon }_I}{\partial x}+\frac{\partial {\phi }^J}{\partial y}\frac{\partial {\upsilon }_I}{\partial y}+\frac{\partial {\phi }^J}{\partial z}\frac{\partial {\upsilon }_I}{\partial z}\right)\stackrel{\wedge }{{\theta }^J}d{\Omega }_s\hfill \\ \hfill & \hspace{1em}\hspace{1em}+{\displaystyle \sum }_{J=1}^M{{\displaystyle \int }}_{{\Gamma }_3}^{}h{\phi }^J\stackrel{\wedge }{{\theta }^J}{\upsilon }_{I}d{\Gamma }_3\hfill \\ \hfill & \hspace{1em}\hspace{1em}-\hspace{0.17em}{\displaystyle \sum }_{J=1}^M{{\displaystyle \int }}_{{\Gamma }_1}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}{n}_x+\frac{\partial {\phi }^J}{\partial y}{n}_y+\frac{\partial {\phi }^J}{\partial z}{n}_z\right)\stackrel{\wedge }{{\theta }^J}{\upsilon }_{I}d{\Gamma }_1\hfill \\ \hfill & \hspace{1em}\hspace{1em}-\hspace{0.17em}{\displaystyle \sum }_{J=1}^M{{\displaystyle \int }}_{L_s}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}{n}_x+\frac{\partial {\phi }^J}{\partial y}{n}_y+\frac{\partial {\phi }^J}{\partial z}{n}_z\right)\stackrel{\wedge }{{\theta }^J}{\upsilon }_{I}d{L}_s\hfill \\ \hfill & \hspace{1em}\hspace{1em}+\alpha {\displaystyle \sum }_{J=1}^M{{\displaystyle \int }}_{{\Gamma }_1}^{}{\phi }^J\stackrel{\wedge }{{\theta }^J}{\upsilon }_{I}d{\Gamma }_1=-{{\displaystyle \int }}_{{\Omega }_s}^{}\dot{\psi }{\upsilon }_{I}d{\Omega }_s\hfill \\ \hfill & \hspace{1em}\hspace{1em}-{{\displaystyle \int }}_{{\Gamma }_2}^{}q{\upsilon }_{I}d{\Gamma }_2+{{\displaystyle \int }}_{{\Gamma }_3}^{}h{\theta }_f{\upsilon }_{I}d{\Gamma }_3+\alpha {{\displaystyle \int }}_{{\Gamma }_1}^{}{\theta }_1{\upsilon }_{I}d{\Gamma }_1.\hfill \end{array}$

(9)

In (9), M is the total number of nodes in the entire domain Ω, $\stackrel{\wedge }{{\theta }^J}$ is the unknown fictitious nodal value, and α is the penalty parameter.

The matrix form of (9) is as follows:

$\mathit{K}\cdot \stackrel{\wedge }{\mathit{\theta }}=\mathit{F}.$

(10)

K and F are the global stiffness matrix and the global vector, respectively, defined as

$\begin{array}{cc}\hfill K_{IJ}& ={{\displaystyle \int }}_{{\Omega }_s}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}\frac{\partial {\upsilon }_I}{\partial x}+\frac{\partial {\phi }^J}{\partial y}\frac{\partial {\upsilon }_I}{\partial y}+\frac{\partial {\phi }^J}{\partial z}\frac{\partial {\upsilon }_I}{\partial z}\right)d{\Omega }_s\hfill \\ \hfill & \hspace{1em}+{{\displaystyle \int }}_{{\Gamma }_3}^{}h{\phi }^J{\upsilon }_{I}d{\Gamma }_3\hfill \\ \hfill & \hspace{1em}-\hspace{0.17em}{{\displaystyle \int }}_{{\Gamma }_1}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}{n}_x+\frac{\partial {\phi }^J}{\partial y}{n}_y+\frac{\partial {\phi }^J}{\partial z}{n}_z\right){\upsilon }_{I}d{\Gamma }_1\hfill \\ \hfill & \hspace{1em}-\hspace{0.17em}{{\displaystyle \int }}_{L_s}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}{n}_x+\frac{\partial {\phi }^J}{\partial y}{n}_y+\frac{\partial {\phi }^J}{\partial z}{n}_z\right){\upsilon }_{I}d{L}_s\hfill \\ \hfill & \hspace{1em}+\alpha {{\displaystyle \int }}_{{\Gamma }_1}^{}{\phi }^J{\upsilon }_{I}d{\Gamma }_1,\hfill \\ \hfill F_I& =-{{\displaystyle \int }}_{{\Omega }_s}^{}\dot{\psi }{\upsilon }_{I}d{\Omega }_s-{{\displaystyle \int }}_{{\Gamma }_2}^{}q{\upsilon }_{I}d{\Gamma }_2\hfill \\ \hfill & \hspace{1em}+{{\displaystyle \int }}_{{\Gamma }_3}^{}h{\theta }_f{\upsilon }_{I}d{\Gamma }_3+\alpha {{\displaystyle \int }}_{{\Gamma }_1}^{}{\theta }_1{\upsilon }_{I}d{\Gamma }_1.\hfill \end{array}$

(11)

If a Heaviside step function is used as the test function for the nodes on the natural boundary or inside the domain, (11) may simplify as

$\begin{array}{cc}\hfill K_{IJ}& ={{\displaystyle \int }}_{{\Gamma }_3}^{}h{\phi }^{J}d{\Gamma }_3-{{\displaystyle \int }}_{{\Gamma }_1}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}{n}_x+\frac{\partial {\phi }^J}{\partial y}{n}_y+\frac{\partial {\phi }^J}{\partial z}{n}_z\right)d{\Gamma }_1\hfill \\ \hfill & \hspace{1em}-{{\displaystyle \int }}_{L_s}^{}\lambda \left(\frac{\partial {\phi }^J}{\partial x}{n}_x+\frac{\partial {\phi }^J}{\partial y}{n}_y+\frac{\partial {\phi }^J}{\partial z}{n}_z\right)d{L}_s+\alpha {{\displaystyle \int }}_{{\Gamma }_1}^{}{\phi }^{J}d{\Gamma }_1,\hfill \\ \hfill F_I& =-{{\displaystyle \int }}_{{\Omega }_s}^{}\dot{\psi }d{\Omega }_s-{{\displaystyle \int }}_{{\Gamma }_2}^{}qd{\Gamma }_2+{{\displaystyle \int }}_{{\Gamma }_3}^{}h{\theta }_{f}d{\Gamma }_3+\alpha {{\displaystyle \int }}_{{\Gamma }_1}^{}{\theta }_{1}d{\Gamma }_1.\hfill \end{array}$

(12)

3. The MLS Approximation Scheme

The MLS, originated by mathematicians for data fitting and surface construction, can be categorized as a method of finite series representation of functions. The MLS method is now a widely used alternative for constructing mesh-free shape functions for approximation.

The MLS approximation has two major features making it popular: (a) the approximated field function is continuous and smooth in the entire problem domain; (b) it is capable of producing an approximation with the desired order of consistency. The MLS approximation is described in the following.

Consider a domain of definition of point x or support domain Ω _x , which is located within the problem domain Ω (see Figure 2) and has a number of randomly located nodes x _I (I=1,…,N). The moving least squares approximate θ ^h (x) of θ(x) by the following definition:

${\theta }^h\left(x\right)={\displaystyle \sum }_{i=1}^m{p}_i\left(x\right)a_i\left(x\right)={\mathit{p}}^T\mathit{a}\left(x\right),$

(13)

where p ^T (x) = [p₁(x),p₂(x),…,p _m (x)] is a complete monomial basis, m is the number of terms in the basis, and a (x) = [a₁(x),a₂(x),…,a _m (x)] is the corresponding coefficient. For example, for a 3D problem, the basis can be chosen as

$\left(m=4\right):{\mathit{p}}^{\mathit{T}}\left(x\right)=\left[1,x,y,z\right],$

(14)

while the quadratic one is

$\left(m=10\right):{\mathit{p}}^{\mathit{T}}\left(x\right)=\left[1,x,y,z,x^2,y^2,z^2,xy,yz,xz\right].$

(15)

The coefficient vector a (x) is determined by minimizing the difference between the local approximation and the function, and it is defined as

$\begin{array}{cc}\hfill J\left(\mathit{a}\left(x\right)\right)& ={\displaystyle \sum }_{I=1}^N{w}_I\left(x\right){\left[{\mathit{P}}^T\left(x_I\right)\mathit{a}\left(x\right)-{\widehat{\theta }}^I\right]}^{\hspace{0.17em}2}\hfill \\ \hfill & ={\left[\mathit{P}\mathit{a}\left(x\right)-\widehat{\mathit{\theta }}\right]}^T\mathit{W}\left[\mathit{P}\mathit{a}\left(x\right)-\widehat{\mathit{\theta }}\right],\hfill \end{array}$

(16)

where x _I denotes the position vector of node I;w _I (x) is the weight function associated with the node I. N is the number of nodes in Ω for which the weight functions w _I (x)>0. The matrix P and the diagonal matrix W are defined as

$\begin{array}{c}\mathit{P}={\left[\begin{array}{cccc}{p}_1\left(x_1\right)& p_2\left(x_1\right)& \cdots & p_m\left(x_1\right)\\ p_1\left(x_2\right)& p_2\left(x_2\right)& \cdots & p_m\left(x_1\right)\\ & \cdot & & \\ & \cdot & & \\ & \cdot & & \\ p_1\left(x_N\right)& p_2\left(x_N\right)& \cdots & p_m\left(x_N\right)\end{array}\right]}_{N\times m},\end{array}$

(17)

$\begin{array}{c}\mathit{W}={\left[\begin{array}{ccccc}{w}_1\left(x\right)& \cdot & \cdot & \cdot & 0\\ \cdot & \cdot & & & \cdot \\ \cdot & & \cdot & & \cdot \\ \cdot & & & \cdot & \cdot \\ 0& \cdot & \cdot & \cdot & w_N\left(x\right)\end{array}\right]}_{N\times N},\end{array}$

(18)

$\begin{array}{c}{\widehat{\theta }}^T={\lfloor {\widehat{\theta }}^1,{\widehat{\theta }}^2,\dots ,{\widehat{\theta }}^N\rfloor }_{1\times N}.\end{array}$

(19)

In (19), ${\stackrel{\wedge }{\theta }}^I$ is the fictitious nodal value. It is not the nodal value of trial functions denoted by θ ^h (x).

In order to obtain the coefficient a (x), we make J( a (x)) in (24) stationary with respect to a (x) as follows:

$\begin{array}{c}\frac{\partial J\left(\mathit{a}\left(x\right)\right)}{\partial \left(\mathit{a}\left(x\right)\right)}=2{\displaystyle \sum }_{I=1}^N{w}_I\left(x\right)\left[{\displaystyle \sum }_{i=1}^m{p}_i\left(x_I\right)a\left(x\right)-{\widehat{\theta }}^I\right]p_i\left(x_I\right)=0\end{array}.$

(20)

This leads to the following set of linear relations:

$\begin{array}{c}\mathit{A}\left(x\right)\mathit{a}\left(x\right)=\mathit{B}\left(x\right)\widehat{\mathit{\theta }},\\ \left(m\times m\right)\left(m\times 1\right)=\left(m\times N\right)\left(N\times 1\right),\end{array}$

(21)

where the matrixes A (x) and B (x) are defined by

$\begin{array}{c}\mathit{A}\left(x\right)={\mathit{P}}^T\mathit{W}\mathit{P}=\mathit{B}\left(x\right)\mathit{P}={\displaystyle \sum }_{I=1}^N{w}_I\left(x\right)\mathit{p}\left(x_I\right){\mathit{p}}^T\left(x_I\right),\\ \mathit{B}\left(x\right)={\mathit{P}}^T\mathit{W}=\left[w_1\left(x\right)\mathit{p}\left(x_1\right),w_2\left(x\right)\mathit{p}\left(x_2\right),\dots ,w_N\left(x\right)\mathit{p}\left(x_N\right)\right].\end{array}$

(22)

Solving a (x) from (21) and substituting it into (13), we can obtain the final form of the MLS approximation as

$\begin{array}{c}{\theta }^h\left(x\right)={\mathbf{\Phi }}^T\left(x\right)\widehat{\mathit{\theta }}={\displaystyle \sum }_{I=1}^N{\varphi }^I\left(x\right){\widehat{\theta }}^I,\\ {\theta }^h\left(x_I\right)={\theta }^I\ne {\widehat{\theta }}^I,\hspace{1em}x∊{\Omega }_x,\end{array}$

(23)

where Φ ^T (x) = p ^T (x) A ^{− 1}(x) B (x) is the shape function, and its partial derivative is

${\varphi }_{,K}^I={\displaystyle \sum }_{j=1}^m\left[p_{j,k}{\left({\mathit{A}}^{-1}\mathit{B}\right)}_{jI}+p_j{\left({\mathit{A}}^{-1}{\mathit{B}}_{,K}+{\mathit{A}}_{,K}^{-1}\mathit{B}\right)}_{jI}\right].$

(24)

In practical applications, the weight function w _I (x) is generally nonzero over the small neighborhood of node x _I , and this neighborhood is called the domain of influence or the domain of definition (Figure 2). Typically, the shape of the domain in the two-dimensional space can be circular, elliptic, rectangular, or any other convenient regular closed line and in the three-dimensional space can be sphere, ellipsoid, cube, or any other convenient regular closed surface. In the present analysis, a spherical support domain has been selected. The choice of weight function w _I (x) affects the resulting approximation θ ^h (x); therefore, the selection of the appropriate weight function is crucially important. The quadratic spline weight function is used in this analysis as shown to be effective in numerical practices [1, 2]. Thus we have

$w_I\left(x\right)=\{\begin{array}{cc}1-6{\left(\frac{d_I}{r_I}\right)}^2+8{\left(\frac{d_I}{r_I}\right)}^3-3{\left(\frac{d_I}{r_I}\right)}^4,\hfill & 0\le d_I\le r_I,\hfill \\ 0,\hfill & d_I\ge r_I,\hfill \end{array}$

(25)

where d _I is the distance between points x and nod x _I , and r _I is the size of the support domain (Figure 2). It can be seen that the quadratic spline weight function is C¹ continuous over the entire domain.

4. Results of Numerical Examples

To study the accuracy and efficiency of the MLPG method, several examples are given in this section. Numerical results obtained from the present approach are compared with the exact solutions, and excellent agreement is achieved.

4.1. Illustrative Example 1

Consider the problem of the 3D heat conduction (1) with the following boundary conditions:

$\begin{array}{c}\theta =y+z+0.5\hspace{1em}\text{at}\hspace{0.17em}x=0,\\ \theta =a+y+z+0.5\hspace{1em}\text{at}\hspace{0.17em}x=a,\\ \theta =x+z+0.5\hspace{1em}\text{at}\hspace{0.17em}y=0,\\ q=-\lambda \hspace{1em}\text{at}\hspace{0.17em}y=b,\\ \theta =-\frac{c}{2}+x+y+0.5\hspace{1em}\text{at}\hspace{0.17em}z=-\frac{c}{2},\\ \theta =\frac{c}{2}+x+y+0.5\hspace{1em}\text{at}\hspace{0.17em}z=\frac{c}{2}.\end{array}$

(26)

The exact solutions for this problem are as follows:

$\begin{array}{c}\theta \left(x,y,z\right)=x+y+z+0.5,\\ q_x\left(x,y,z\right)=-\lambda ,\\ q_y\left(x,y,z\right)=-\lambda ,\\ q_z\left(x,y,z\right)=-\lambda .\end{array}$

(27)

The geometry for this example is shown in Figure 3. The node distribution with the grid of 64 nodes is presented in Figure 4 for the case of a=b=c=1. The temperature distribution diagrams along two lines are exhibited in Figures 5 and 6.

Figure 3:

Geometry used for example 1 ( $\dot{\psi }=0\hspace{0.17em}\text{and}\hspace{0.17em}\lambda =50$ ).

Figure 4:

Regular node distribution for example 1.

Figure 5:

Temperature distribution along $x=1/2$ and $y=1/2$ .

Figure 6:

Temperature distribution along $y=1/2$ and $z=1/4$ .

As it is observed from these figures, the MLPG results agree well with the exact solutions of the problem.

4.2. Illustrative Example 2

Consider the problem of the 3D heat conduction (1) with the following boundary conditions:

$\begin{array}{c}\theta =0\hspace{1em}\text{at}\hspace{0.17em}x=0,\\ \theta =\sinh \left(\frac{\pi a}{2}\right)\sin \left(\frac{\pi y}{2\sqrt{2}}\right)\sin \left(\frac{\pi \left(z+\left(c/2\right)\right)}{2\sqrt{2}}\right)\hspace{1em}\text{at}\hspace{0.17em}x=a,\\ \theta =0\hspace{1em}\text{at}\hspace{0.17em}y=0,\\ \theta =\sinh \left(\frac{\pi x}{2}\right)\sin \left(\frac{\pi b}{2\sqrt{2}}\right)\sin \left(\frac{\pi \left(z+\left(c/2\right)\right)}{2\sqrt{2}}\right)\hspace{1em}\text{at}\hspace{0.17em}y=b,\\ \theta =0\hspace{1em}\text{at}\hspace{0.17em}z=-\frac{c}{2},\\ \theta =\sinh \left(\frac{\pi x}{2}\right)\sin \left(\frac{\pi y}{2\sqrt{2}}\right)\sin \left(\frac{\pi c}{2\sqrt{2}}\right)\hspace{1em}\text{at}\hspace{0.17em}z=\frac{c}{2},\end{array}$

(28)

with the same geometry and node distribution as shown in Figures 3 and 4.

The exact solutions for this problem are

$\begin{array}{cc}\hfill \theta \left(x,y,z\right)& =\sinh \left(\frac{\pi x}{2}\right)\sin \left(\frac{\pi y}{2\sqrt{2}}\right)\sin \left(\frac{\pi \left(z+0.5\right)}{2\sqrt{2}}\right),\hfill \\ \hfill q_x\left(x,y,z\right)& =-\lambda \frac{\pi }{2}\cosh \left(\frac{\pi x}{2}\right)\sin \left(\frac{\pi y}{2\sqrt{2}}\right)\sin \left(\frac{\pi \left(z+0.5\right)}{2\sqrt{2}}\right),\hfill \\ \hfill q_y\left(x,y,z\right)& =-\lambda \frac{\pi }{2\sqrt{2}}\sinh \left(\frac{\pi x}{2}\right)\cos \left(\frac{\pi y}{2\sqrt{2}}\right)\sin \left(\frac{\pi \left(z+0.5\right)}{2\sqrt{2}}\right),\hfill \\ \hfill q_z\left(x,y,z\right)& =-\lambda \frac{\pi }{2\sqrt{2}}\sinh \left(\frac{\pi x}{2}\right)\sin \left(\frac{\pi y}{2\sqrt{2}}\right)\cos \left(\frac{\pi \left(z+0.5\right)}{2\sqrt{2}}\right).\hfill \end{array}$

(29)

The temperature distributions and the heat flux distributions are presented in Figures 7, 8, 9, 10, 11, 12, and 13. As shown in these figures, the MLPG results admit a remarkable accuracy.

Figure 7:

Temperature distribution along $x=1/2$ and $y=1/2$ .

Figure 8:

Heat flux distribution along $x=1/2$ and $y=1/2$ .

Figure 9:

Heat flux distribution along $x=1/2$ and $\hspace{0.17em}y=1/2$ .

Figure 10:

Temperature distribution along $y=1/2$ and $z=1/4$ .

Figure 11:

Heat flux distribution along $y=1/2$ and $z=1/4$ .

Figure 12:

Heat flux distribution along $y=1/2$ and $z=1/4$ .

Figure 13:

Heat flux distribution along $y=1/2$ and $z=1/4$ .

4.3. Illustrative Example 3

Consider the problem of the 3D heat conduction ((1)) with the following boundary conditions:

$\begin{array}{c}\theta =-y^2+z^2-1\hspace{1em}\text{at}\hspace{0.17em}x=0,\\ \theta =a^2-y^2+z^2-1\hspace{1em}\hspace{0.17em}\text{at}\hspace{0.17em}x=a,\\ \theta =x^2+z^2-1\hspace{1em}\text{at}\hspace{0.17em}y=0,\\ \theta =x^2-b^2+z^2-1\hspace{1em}\text{at}\hspace{0.17em}y=b,\\ \theta =x^2-y^2+\frac{c^2}{4}-1\hspace{1em}\text{at}\hspace{0.17em}z=-\frac{c}{2},\\ \theta =x^2-y^2+\frac{c^2}{4}-1\hspace{1em}\text{at}\hspace{0.17em}z=\frac{c}{2}.\end{array}$

(30)

The geometry is the same as in Figure 3 for the case of $\dot{\psi }=30\hspace{0.17em}\text{and}\hspace{0.17em}\lambda =15$ . The node distribution with 27 nodes is exhibited in Figure 14 for the case of a=b=c=1.

Figure 14:

Regular node distribution for example 3.

The exact solutions for this case are as follows:

$\begin{array}{c}\theta \left(x,y,z\right)=x^2-y^2+z^2-1,\\ q_x\left(x,y,z\right)=-2\lambda x,\\ q_y\left(x,y,z\right)=2\lambda y,\\ q_z\left(x,y,z\right)=-2\lambda z.\end{array}$

(31)

The temperature and heat flux distributions are presented in Figures 15, 16, 17, 18, and 19. These figures show that the MLPG results agree with the exact solutions.

Figure 15:

Temperature distribution along $x=1/4$ and $y=3/4$ .

Figure 16:

Heat flux distribution along $x=1/4$ and $y=3/4$ .

Figure 17:

Temperature distribution along $y=1/2$ and $z=1/4$ .

Figure 18:

Heat flux distribution along $y=1/2$ and $z=1/4$ .

Figure 19:

Heat flux distribution along $y=1/2$ and $z=1/4$ .

4.4. Illustrative Example 4

Consider the problem of the 3D heat conduction (1) with the following boundary conditions

$\begin{array}{c}q=0\hspace{1em}\text{at}\hspace{0.17em}x=0,\\ \theta =\cosh \left(\frac{a}{\sqrt{2}}\right)\cosh \left(\frac{y}{\sqrt{2}}\right)\sin \left(z\right)+z^2\hspace{1em}\text{at}\hspace{0.17em}x=a,\\ q=0\hspace{1em}\text{at}\hspace{0.17em}y=0,\\ \theta =\cosh \left(\frac{x}{\sqrt{2}}\right)\cosh \left(\frac{b}{\sqrt{2}}\right)\sin \left(z\right)+z^2\hspace{1em}\text{at}\hspace{0.17em}y=b,\\ \theta =-\cosh \left(\frac{x}{\sqrt{2}}\right)\cosh \left(\frac{y}{\sqrt{2}}\right)\sin \left(\frac{c}{2}\right)+\frac{c^2}{4}\hspace{1em}\text{at}\hspace{0.17em}z=-\frac{c}{2},\\ \theta =\cosh \left(\frac{x}{\sqrt{2}}\right)\cosh \left(\frac{y}{\sqrt{2}}\right)\sin \left(\frac{c}{2}\right)+\frac{c^2}{4}\hspace{1em}\text{at}\hspace{0.17em}z=\frac{c}{2}.\end{array}$

(32)

Here, the geometry is the same as in Figure 3 for the case of $\dot{\psi }=2\pi \hspace{0.17em}\text{and}\hspace{0.17em}\lambda =\pi$ . The values of a,b, and c are set equal to a=b=c=π, and the node distribution with the grid of 64 nodes is considered as shown in Figure 20. The exact solutions for this case are as follows:

$\begin{array}{cc}\hfill \theta \left(x,y,z\right)& =\cosh \left(\frac{x}{\sqrt{2}}\right)\cosh \left(\frac{y}{\sqrt{2}}\right)\sin \left(z\right)+z^2,\hfill \\ \hfill q_x\left(x,y,z\right)& =-\lambda \frac{1}{\sqrt{2}}\sinh \left(\frac{x}{\sqrt{2}}\right)\cosh \left(\frac{y}{\sqrt{2}}\right)\sin \left(z\right),\hfill \\ \hfill q_y\left(x,y,z\right)& =-\lambda \frac{1}{\sqrt{2}}\cosh \left(\frac{x}{\sqrt{2}}\right)\sinh \left(\frac{y}{\sqrt{2}}\right)\sin \left(z\right),\hfill \\ \hfill q_z\left(x,y,z\right)& =-\lambda \cosh \left(\frac{x}{\sqrt{2}}\right)\cosh \left(\frac{y}{\sqrt{2}}\right)\cos \left(z\right)+2z.\hfill \end{array}$

(33)

The temperature and heat flux distributions for this case are also illustrated in Figures 21, 22, 23, and 24. Once again, an excellent agreement is achieved with the exact solutions.

Figure 20:

Regular node distribution for example 4.

Figure 21:

Temperature distribution along $x=\pi /4$ and $y=3\pi /4$ .

Figure 22:

Heat flux distribution along $x=\pi /4$ and $y=3\pi /4$ .

Figure 23:

Heat flux distribution along $x=\pi /4$ and $y=3\pi /4$ .

Figure 24:

Heat flux distribution along $y=\pi /2$ and $z=\pi /4$ .

5. Results Comparison of MLPG and FDM

Also for comparison, the obtained results from the proposed method are illustrated with results of FDM in Tables 1, 2, 3, and 4. In these tables, errors are the error average of non-boundary nodes temperature.

Table 1:

Results’ comparison of MLPG and FDM for example 1.

Method	Error
MLPG	6.1×10− 3
FDM	4.8861×10− 4

Table 2:

Results’ comparison of MLPG and FDM for example 2.

Method	Error
MLPG	6×10− 3
FDM	3.7×10− 3

Table 3:

Results’ comparison of MLPG and FDM for example 3.

Method	Error
MLPG	4.4785×10− 4
FDM	8.33×10− 2

Table 4:

Results’ comparison of MLPG and FDM for example 4.

Method	Error
MLPG	0.1989
FDM	0.5711

It can be seen that, for example 3 and 4, the obtained results from MLPG are better than FDM; also, MLPG shows an acceptable accuracy in example 1 and 2. Nevertheless, MLPG can calculate the problem values in each desired position; however, FDM could calculate these values just in the node points.

6. Conclusion

The meshless local Petrov-Galerkin method is applied to analyze the three-dimensional steady-state heat conduction problems. Using a Heaviside test function, the domain of the integral in the weak form is simplified. This substantially reduces the computational efforts needed for constructing the stiffness matrix. Hence, the present approach is computationally efficient compared to the conventional MLPG method. The penalty approach is used to impose the essential boundary conditions, and the moving least squares approximation is used for interpolation schemes. Having established the accuracy of obtained numerical results, the efficiency as well as the capability of the present technique is clarified. Also, one of the most impotent benefits of MLPG rather than the other methods such as FDM is that MLPG can calculate the problem values in each desired position; however, FDM could calculate these values just in the node points.